poisson.cpp

偏微分方程数值解- 这是矩阵的紧致存储及在紧致存储的矩阵上用共轭梯度(CG)方法解Poisson方程的实现, gcc mingw 3.4.4下通过

#include <cstdlib>
#include <iomanip.h>
#include "CG.h"

#define FUNCTION(x,y) (pi*pi*sin(pi*(x)))
#define BOUNDARY(x,y) 0
#define SOLUTION(x,y) (sin(pi*(x)))

const double pi = 3.141592653589793;

double value(int i, int j, int N){
	double x = double(i)/N, y = double(j)/N;
	return FUNCTION(x,y);
}

double boundary(int i, int j, int N){
	double x = double(i)/N, y = double(j)/N;
	if(i == 1 ){
		if(j==1)
			return BOUNDARY(0,y)+BOUNDARY(x,0);
		else if(j==N-1)
			return BOUNDARY(0,y)+BOUNDARY(x,1);
		else
			return BOUNDARY(0,y);
	} else if(i == N-1){
		if(j==1)
			return BOUNDARY(1,y)+BOUNDARY(1-y,0);
		else if(j==N-1)
			return BOUNDARY(1,y)+BOUNDARY(y,1);
		else
			return BOUNDARY(1,y);
	} else{
		if(j==1)
			return BOUNDARY(x,0);
		else if(j==N-1)
			return BOUNDARY(x,1);
		else
			return 0;
	}
}

void exactSolution(Vector& u,int N){
	u.rescale((N-1)*(N-1));
	for(int i=1;i<N;i++)
		for(int j=1;j<N;j++)
			u((i-1)*(N-1)+j) = SOLUTION(double(i)/N,double(j)/N);
}



bool initialize(Matrix& B, Vector& u, Vector& r, int N){
	double h = 1.0/N;
	int l = (N-1)*(N-1);
	B.rescale(l,N);
	u.rescale(l);
	r.rescale(l);
	for(int i=1;i<N;i++){
		for(int j=1;j<N;j++){
			//B
			if(j==1 || j==N-1) B((i-1)*(N-1)+j, N) = -3;
			else B((i-1)*(N-1)+j, N) = -4;
			if(j!=N-1){
				B((i-1)*(N-1)+j+1, N-1) = 1;
			}
			if((i-2)*(N-1)+j>0){
				B((i-1)*(N-1)+j, 1) = 1;
			}
			//r
			r((i-1)*(N-1)+j) = value(i,j,N) + boundary(i,j,N);
		}
	}
	r *= (h*h*(-1));
}

int main(void){
	int N;
	Matrix  B(1,1),result(1,1); //Always use "B" to denote compact mode of matrix A
	Vector u(1), r(1), exact(1);
	cout<<"Input N:";
	cin>>N;
	cout<<endl;
	initialize(B,u,r,N);
	cout<<"Initialization complete.>>>>>>>>>>>>>>>>>"<<endl;
	cout<<"Matrix B is "<<B.getrow()<<"*"<<B.getcol()<<endl;
	system("PAUSE");
//	int a = SORsolve(B,r,1.0,u);
//	printVector(u);
	int a = CGsolve(B,r,u);
	cout<<"iterations:"<<a<<endl;
//	printVector(u);
	exactSolution(exact, N);
//	printVector(exact);
	cout<<"error:"<<error(u,exact)<<endl<<endl;
	

	result.rescale(N-1,1);
	for(int i=1;i<N;i++)
		result(i,1) = u((i-1)*(N-1)+1);
	
	printMatrix(result);

    system("PAUSE");
    return EXIT_SUCCESS;
}