hvcf

abgabe4nochmal/src/abgabe4nochmal.cpp

@@ -203,23 +203,36 @@ double lanczosArnoldi(gauss_leg_param *glp,
 
 	}
 	MatrixXd A_new(mx, mx);
+    MatrixXcd psi(dim, mx);
+	EigenSolver<MatrixXd> 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<MatrixXd> 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));