Strongly Implicit Procedure (SIP) 代碼

此代碼爲SIP方法求解五對角矩陣，在結構上是對稱的，可是A能夠不是對稱矩陣

int SIPsolver(double *ALL, double *AL, double *AC, double *AR, double *ARR, double *x,double* r, int m,int n)
{
	//求解以下形式矩陣：
	//AC0	AR0			ARR0						
	//AL1	AC1	AR1			ARR1					
	//	    AL2	AC2	AR2			ARR2				
	//		    AL3	AC3	AR3			ARR3			
	//ALLm			AL4	AC4	AR4			ARR4		
	//	ALL5			AL5	AC5	AR5			ARRn-m	
	//		ALL6			AL6	AC6	AR6			
	//			ALL7			AL7	AC7	AR7		
	//				ALL8			AL8	AC8	AR8	
	//					ALLn			ALn	ACn	
	double a=0.6;
	double P1=0;
	double P2=0;
	double *LLL,*LL,*LC;
	double *UR,*URR;
	double *RES,*y;
	double RESN=0;
	int MAXITER=1000;
	int ITER;
	LLL=new double [n+1];
	LL=new double [n+1];
	LC=new double [n+1];
	UR=new double [n+1];
	URR=new double [n+1];
	RES=new double [n+1];
	for(int i=0;i<=n;i++){
		LLL[i]=0;
		LL[i]=0;
		LC[i]=0;
		UR[i]=0;
		URR[i]=0;
		RES[i]=0;
	}
	for(int i=0;i<=n;i++)
	{
		if(i>=m){
			LLL[i]=ALL[i]/(1+a*UR[i-m]);
		}
		if(i>=1){
			LL[i]=AL[i]/(1+a*URR[i-1]);
		}
		if(i>=m){
			P1=a*LLL[i]*UR[i-m];
		}else{
			P1=0;
		}
		if(i>=1){
			P2=a*LL[i]*URR[i-1];
		}else{
			P2=0;
		}
		if(i<1){
			LC[i]=AC[i]+(P1+P2);
		}else if(i>=1&&i<m){
			LC[i]=AC[i]+(P1+P2)-LL[i]*UR[i-1];
		}else{
			LC[i]=AC[i]+(P1+P2)-LLL[i]*URR[i-m]-LL[i]*UR[i-1];
		}
		if(i<=n-1){
			UR[i]=(AR[i]-P1)/LC[i];
		}else{
			UR[i]=0;
		}
		if(i<=n-m){
			URR[i]=(ARR[i]-P2)/LC[i];
		}else{
			URR[i]=0;
		}
		
	}
	ITER=0;
	do{
		RESN=0;
		for(int i=0;i<=n;i++)
		{
			RES[i]=r[i]-(ALL[i]*(i-m>=0?x[i-m]:0)+AL[i]*(i-1>=0?x[i-1]:0)+AC[i]*x[i]+AR[i]*(i+1<=n?x[i+1]:0)+ARR[i]*(i+m<=n?x[i+m]:0));
			RESN+=abs(RES[i]);
		}

		for(int i=0;i<=n;i++)
		{
			RES[i]=(RES[i]-(i-1>=0?RES[i-1]:0)*LL[i]-(i-m>=0?RES[i-m]:0)*LLL[i])/LC[i];
		}
		
		for(int i=n;i>=0;i--)
		{
			x[i]=x[i]+(RES[i]-(i+1<=n?RES[i+1]:0)*UR[i]-(i+m<=n?RES[i+m]:0)*URR[i]);
		}
		ITER++;
	}while(abs(RESN)>1E-6&&ITER<MAXITER);
	
	delete [] LLL;
	delete [] LL;
	delete [] LC;
	delete [] UR;
	delete [] URR;
	delete [] RES;
	return ITER>MAXITER?-1:1;
}

根據論文

Iterative solution of implicit approximations of multidimensional partial differential equations
HL Stone - SIAM Journal on Numerical Analysis, 1968 - SIAM

編寫