crout.cpp

用CROUT分解法求解三对角方程.对于对角占优的三对角方程可以求出数值稳定解.

//author:张鹏飞 ID :00512080  NO.3
#include"stdio.h"
#include"math.h"
#include"string.h"
#define N 100   //N表示程序可以求解的方程组的阶数
int n;

bool inputerror=false; //标识数据输入有没有错误，若取true，则表示有错误

bool methodfailed=false;   //标识Crout分解能否对此方程组求解 
double a[N],b[N],c[N],f[N]; 
char ch;

double chartodouble(char* c)//此函数用于将一个字符数组转换成数字，如果含有非法字符，则返回0，同时向inputerror报错
{
	int i,flag=0,negetive;
	double num=0;
	if(c[0]=='-'){negetive=-1;i=1;}
	else if(c[0]=='+'){negetive=1;i=1;}
	else if(c[0]>='0'&&c[0]<='9'){negetive=1;i=0;}
	else {inputerror=true;return 0;}
	while(c[i]!='\0'&&flag==0)
	{
		if(c[i]=='.'){flag=1;i++;break;}
		else if(c[i]<'0'||c[i]>'9')
		{inputerror=true;return 0;}
		else {num=(double)(c[i]-'0')+num*10.0;}	
		i++;
	}	
	double d=0.1;
	while(c[i]!='\0'&&flag==1)
	{
	
		if(c[i]<'0'||c[i]>'9'){inputerror=true;return 0;}
		else
		{
			num=d*(double)(c[i]-'0')+num;
			d=0.1*d;
		}
		i++;
		
	}
	return num*negetive;
}
	
bool getdata(FILE* fp)
{	                  //此函数用于从fp指向的文件中读取方程的系数，并传递给相应的数组，
                      //同时兼容了输入文件格式的不规范可能引起的一些错误
	char cha[20];
	int i=0,j;
	while(ch!=' '&&ch!='\n')
	{
		cha[i]=ch;
		i++;
		ch=fgetc(fp);
	}
	cha[i]='\0';
	n=(int)chartodouble(cha);
  if(n==-1)return false;
  else if(n<-1) inputerror=true;
  else
  {
	while(ch==' '||ch=='\n')
	{
		ch=fgetc(fp);
	}
	//input a[n]
	for(j=1;j<n;j++)
	{
		i=0;
		while(ch!=' '&&ch!='\n')
		{
		   cha[i]=ch;
		   i++;
		   ch=fgetc(fp);
		}
        cha[i]='\0';
	    a[j]=chartodouble(cha);
		if(ch=='\n'){ch=fgetc(fp);break;}
		ch=fgetc(fp);
		while(ch==' ')
		{
		  ch=fgetc(fp);
		}
		if(ch=='\n'){ch=fgetc(fp);break;}
	}
	if(j<n-1)inputerror=1;

	//input b[n]
	for(j=0;j<n;j++)
	{
		i=0;
		while(ch!=' '&&ch!='\n')
		{
		   cha[i]=ch;
		   i++;
		   ch=fgetc(fp);
		}
        cha[i]='\0';
	    b[j]=chartodouble(cha);
		if(ch=='\n'){ch=fgetc(fp);break;}
		ch=fgetc(fp);
		while(ch==' ')
		{
		  ch=fgetc(fp);
		}
		if(ch=='\n'){ch=fgetc(fp);break;}
	}
	if(j<n-1)inputerror=true;

	//input c[n]
	for(j=0;j<n-1;j++)
	{
		i=0;
		while(ch!=' '&&ch!='\n')
		{
		   cha[i]=ch;
		   i++;
		   ch=fgetc(fp);
		}
        cha[i]='\0';
	    c[j]=chartodouble(cha);
		if(ch=='\n'){ch=fgetc(fp);break;}
		ch=fgetc(fp);
		while(ch==' ')
		{
		  ch=fgetc(fp);
		}
		if(ch=='\n'){ch=fgetc(fp);break;}
	}
	if(j<n-2)inputerror=true;

	//input f[n]
	for(j=0;j<n;j++)
	{
		i=0;
		while(ch!=' '&&ch!='\n')
		{
		   cha[i]=ch;
		   i++;
		   ch=fgetc(fp);
		}
        cha[i]='\0';
	    f[j]=chartodouble(cha);
		if(ch=='\n'){ch=fgetc(fp);break;}
		ch=fgetc(fp);
		while(ch==' ')
		{
		  ch=fgetc(fp);
		}
		if(ch=='\n'){ch=fgetc(fp);break;}
	}
	if(j<n-1)inputerror=true;
	return true;
  }
}



void crout(double a[],double b[],double c[],double f[],int n)//用Crout分解求解三对角方程，对于不能用Crout分解法求解的方程，向 methodfailed函数报错
{
	 c[0]=c[0]/b[0];
	 for(int i=1;i<n-1;i++)
	 {
		 b[i]=-c[i-1]*a[i]+b[i];
		 if(b[i]==0){methodfailed=true;return;}
		 c[i]=c[i]/b[i];
	 }
		 b[n-1]=-c[n-2]*a[n-1]+b[n-1];
		 if(b[n-1]==0){methodfailed=true;return;}
	 f[0]=f[0]/b[0];
	 for(i=1;i<n;i++)
	 {
		 f[i]=(f[i]-a[i]*f[i-1])/b[i];
	 }
	 for(i=n-2;i>=0;i--)
	 {
		 f[i]=f[i]-c[i]*f[i+1];
	 }
}
void main()
{
	FILE *fp1,*fp2;
	int i;
	if((fp1=fopen("input.txt","r"))==NULL)
	{
		printf("cannot open input.txt!\n");
		return;
	}
	if((fp2=fopen("output.txt","w"))==NULL)
	{
		printf("cannot open output.txt\n");
	}
	ch=fgetc(fp1);
	while(getdata(fp1))
	{
	    crout(a,b,c,f,n);
        if(inputerror==true)
		{
			fprintf(fp2,"inputdata error!\n");
		}
	    else if(methodfailed==true)
		{
			fprintf(fp2,"The Crout method failed.\n");
		}
        else
		{	
		    for(i=0;i<n;i++)
			{ 
	    	   fprintf(fp2,"%16.12f\n",f[i]);
			}
	   
		}
		fprintf(fp2,"\n");
	}
	fclose(fp1);
	fclose(fp2);
}