diff --git a/abgabe4nochmal/src/abgabe4nochmal.cpp b/abgabe4nochmal/src/abgabe4nochmal.cpp index c78e48b..5492ac2 100644 --- a/abgabe4nochmal/src/abgabe4nochmal.cpp +++ b/abgabe4nochmal/src/abgabe4nochmal.cpp @@ -203,23 +203,36 @@ double lanczosArnoldi(gauss_leg_param *glp, } MatrixXd A_new(mx, mx); + MatrixXcd psi(dim, mx); + EigenSolver es(A_new); + MatrixXcd EV_new = es.eigenvectors(); for (int i = 0; i < mx; i++) { for (int j = 0; j < mx; j++) { A_new(i, j) = eigenvektoren[i].dot(A * eigenvektoren[j]); //printf("%f\n",A_new(i,j)); } + } - //VectorXcd psi[mx](dim); - MatrixXcd psi(dim, mx); - EigenSolver es(A_new); - MatrixXcd EV_new = es.eigenvectors(); + double lambda[mx]; + double L = 0; for (int j = 0; j < mx; j++) { for (int k = 0; k < dim; k++) { psi(k, j) = 0; } + for (int i = 0; i < mx; i++) { + psi.col(i) += EV_new(j, i) * eigenvektoren[j]; + psi.col(i) = psi.col(i) / (psi.col(i).norm()); + lambda[i] = psi.col(i).dot(A * psi.col(i)).real(); + + if (lambda[i] > L||i==0) { + L = lambda[i]; + } + + } } + /* for (int j = 0; j < mx; j++) { for (int i = 0; i < mx; i++) { @@ -227,7 +240,8 @@ double lanczosArnoldi(gauss_leg_param *glp, //printf("%f\n",psi(i,j).real()); } - } + }*/ + /* double lambda[mx]; double L = 0; psi.col(0) = psi.col(0) / (psi.col(0).norm()); @@ -241,7 +255,7 @@ double lanczosArnoldi(gauss_leg_param *glp, if (lambda[i] > L) { L = lambda[i]; } - } + }*/ return L; } @@ -312,7 +326,7 @@ void variationH(gauss_leg_param *glp, double e1, double e2, int main() { gauss_leg_param glp; glp.mu = 1; - glp.n = 2000; + glp.n = 1000; glp.v0 = 400; glp.x1 = 0; glp.x2 = 10; @@ -320,7 +334,7 @@ int main() { glp.p2 = 0.5; glp.E = -145; glp.pmax = 200; - glp.h = 0.1; + glp.h = 0.5; glp.a = 1; getStuetzstellen(&glp); //printf("%f\t", gauss_Int(&glp, potential)); @@ -328,10 +342,10 @@ int main() { //printf("%f", abweichung(&glp, potential, potExact)); //lanczosArnoldi(&glp, potential); - //printf("%f",secant(&glp,-200,-100,potential)); + printf("%f",secant(&glp,-200,-100,potential)); glp.x1 = 0.4; - printf("%f", secant(&glp, -200, -100, potential2)); + //printf("%f", secant(&glp, -200, -100, potential2)); //printf("%15.6e\t", gauss_Int(&glp, potential2)); glp.a = 2; //printf("%f", secant(&glp, -200, -100, potExact));