C++ MatrixXf::colPivHouseholderQr方法代码示例

void TestLeastSquares() {
	MatrixXf A = MatrixXf::Random(10, 2);
	VectorXf b = VectorXf::Random(10);
	Vector2f x;

	std::cout << "=============================" << std::endl;
	std::cout << "Testing least squares solvers" << std::endl;
	std::cout << "=============================" << std::endl;

	x = A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b);
	std::cout << "Solution using Jacobi SVD = " << x.transpose() << std::endl;

	x = A.colPivHouseholderQr().solve(b);
	std::cout << "Solution using column pivoting Householder QR = " << x.transpose() << std::endl;

	// If the matrix A is ill-conditioned, then this is not a good method
	x = (A.transpose() * A).ldlt().solve(A.transpose() * b);
	std::cout << "Solution using normal equation = " << x.transpose() << std::endl;
}

int main(int argc, char **argv){

   double t = atof(argv[1]);
   float time = 0;
   float x_values [x_spaces+1];
   //float u_values [x_spaces+1];
   float u_old[x_spaces+1];
	float correction_factor;

	for (int i=0; i<=x_spaces;i++){
	   x_values[i] = x_lower + (x_higher-x_lower)/x_spaces * i;
	   u_old[i] = func(x_values[i]);
	}


   MatrixXf A = MatrixXf::Zero(x_spaces+1,x_spaces+1);
   VectorXf b(x_spaces+1,1);
   VectorXf u_values(x_spaces+1,1);

// The next section deals with the euler explicit method

   while (time <= t) {

      for (int i=0;i<=x_spaces;i++){
         if (i-1 >= 0){
            A(i,i-1) = -delta_t/(4*delta_x)*(u_old[i-1]);
         }

         if (i+1 <= x_spaces){
            A(i,i+1) = delta_t/(4*delta_x)*(u_old[i+1]); 
         }
         A(i,i) = 1;

			if (i-2<0){
				correction_factor = -omega/8*(u_old[i+2] -4*u_old[i+1]+6*u_old[i]);
			}else if (i-1<0){
				correction_factor = -omega/8*(u_old[i+2] -4*u_old[i+1]+6*u_old[i]-4*u_old[i-1]);
			}else if (i+1>x_spaces){
				correction_factor = -omega/8*(-4*u_old[i+1]+6*u_old[i]-4*u_old[i-1]+u_old[i-2]);
			}else if (i+2>x_spaces){
				correction_factor = -omega/8*(6*u_old[i]-4*u_old[i-1]+u_old[i-2]);
			}else{
				correction_factor = -omega/8*(u_old[i+2] -4*u_old[i+1]+6*u_old[i]-4*u_old[i-1]+u_old[i-2]);
			}


         if (i-1 < 0){
         	b(i) = -delta_t/(2*delta_x)*(pow(u_old[i+1],2)/2-pow(ustart,2)/2)-delta_t/(4*delta_x)*pow(ustart,2) + u_old[i]+delta_t/(4*delta_x)*pow(u_old[i+1],2)+correction_factor;
			}else if (i+1>x_spaces){
				b(i) = -delta_t/(2*delta_x)*(pow(uend,2)/2-pow(u_old[i-1],2)/2)-delta_t/(4*delta_x)*pow(u_old[i-1],2) + u_old[i]+delta_t/(4*delta_x)*pow(uend,2)+correction_factor;
			}else{
         	b(i) = -delta_t/(2*delta_x)*(pow(u_old[i+1],2)/2-pow(u_old[i-1],2)/2)-delta_t/(4*delta_x)*pow(u_old[i-1],2) + u_old[i]+delta_t/(4*delta_x)*pow(u_old[i+1],2)+correction_factor;
			}
     
      }
      //b(x_spaces) = u_old[x_spaces] - v_Neumann /2*u_old[x_spaces]; 
      u_values = A.colPivHouseholderQr().solve(b);

		time = time + delta_t;

		if (time >= t){
			break;
		}
		for (int i=0;i<=x_spaces;i++){
		   u_old[i] = u_values(i);
		}
		u_old[0]=ustart;
		u_values[0] = ustart;
		u_old[x_spaces]=uend;
		u_old[x_spaces]=uend;

  }    

		for (int j=0; j<=x_spaces;j++){
		   if (j==0){
		      cout << "X Implicit Value Start" << "," << x_values[j] << ",";
		   }else if (j==x_spaces){
		      cout << x_values[j] << "," << "X Implicit Value End" << ",";
		   }else{
		      cout << x_values[j] << ",";
		   }
		
		} 

		for (int k=0; k<=x_spaces;k++){
		   if (k==0){
		      cout << "Y Implicit Value Start" << "," << u_old[k] << ",";
		   }else if (k==x_spaces){
		      cout << u_old[k] << "," << "Y Implicit Value End"<<","<<endl;
		   }else{
		      cout << u_old[k] << ",";
		   }
	  
}

 
 



//.........这里部分代码省略.........

int main(int argc, char **argv){
   int n = atof(argv[1]);
	string method = argv[2];
   double h = 1.0/(n+1.0);
	double max_error = pow(10,-6);
	double error =1;
   time_t t_start,t_end;
	double  dif_iterative, dif_matrix;
   VectorXf  x (n,1);
   VectorXf  xk (n,1);
   VectorXf  xkplus1 (n,1);
   VectorXf  f (n,1);
   MatrixXf A = MatrixXf::Zero(n,n);
   MatrixXf I = MatrixXf::Zero(n,n);
   VectorXf  b (n,1);
   VectorXf xc (n,1);

   for (int i=0; i<n;i++){
      x(i) = h*(i+1);
		f(i) = func(x(i));
   }

	for (int i=0;i<n;i++){
		A(i,i) = 2;
		if (i!=n-1){
			A(i,i+1) = -1;
			A(i+1,i) = -1;
		} 
	}
	b = f*pow(h,2);

	//Solving the Equation
   t_start = time(0);
   xc = A.colPivHouseholderQr().solve(b);
  	t_end=time(0);
  	dif_matrix = difftime (t_end,t_start);

/*
   cout << "Here is the matrix A:\n" << A << endl;
   cout << "Here is the vector x:\n" << x << endl;
   cout << "Here is the vector b:\n" << b << endl;
   cout << "Here is the vector uc:\n" << xc << endl;
*/
	//Jacobin part
	MatrixXf P = MatrixXf::Zero(n,n);
	MatrixXf Pinverse = MatrixXf::Zero(n,n);

	if (method == "Jacobian"){
		for (int i=0; i<n;i++){
      	P(i,i) = A(i,i);
   	}
	} else if (method=="Gauss-Seidel"){
		for (int i=0; i<n;i++){
      	P(i,i) = A(i,i);
			for (int j=0; j<i;j++){
				P(i,j) = A(i,j); 
			}
   	}
	} else if (method=="SOR"){
		double omega = atof(argv[3]);
		for (int i=0; i<n;i++){
      	P(i,i) = A(i,i);
			for (int j=0; j<i;j++){
				P(i,j) = omega*A(i,j); 
			}
   	}

	} else {
		return 1;
	}

	Pinverse = P.inverse();

   for (int i=0; i<n;i++){
      xk(i) = 0;
		I(i,i) = 1;
   }


   t_start = time(0);
	double repeat_error=100*n;
	int loop=1;
	while(error >=max_error){
		xkplus1 = (I-Pinverse*A)*xk+Pinverse*b;
		xk = xkplus1;
		error = maximum_array(xc,xkplus1,n);
		//cout << loop<<endl;
		//cout << error<<endl;

		if (repeat_error == error){
			break;
		}
		repeat_error = error;
		loop +=1;
	}
  t_end=time(0);
	dif_iterative = difftime (t_end,t_start);
	cout << dif_iterative <<"," <<loop<<"," << dif_matrix <<","<<maximum_array(xc,xkplus1,n);

   return 0; 
//.........这里部分代码省略.........