微分方程差分三对角矩阵解法C语言程序

#include<stdio.h>
#include<math.h>
#define N 1000
//三对角方程组Thomas算法，计算两点边值问题
//-u"(x)+u(x)=exp(x)(sin(x)-cos(x),0<=x<=pi
//u(00=0, u(pi)=0.
void ThomasMethodsf(double a[],double b[],double c[],double d[],int n);
//Thomas算法函数
void Errorf(double x[],double d[],double e[],double *maxe_h,int n);
//计算结点处数值解的误差的绝对值和数值解的最大误差函数
void prtf(double a[],double b[],double c[],double d[],int n);
//打印三对角方程组系数矩阵函数
double fi(double xi);//微分方程右端项
//double qi(double xi);//微分方程u(x)系数
double ui(double xi);//微分方程精确解u（x)
void res(double d[],double e[],double maxe_h,int n);//输出结果到Result.txt文档
void main()
{
double a[N];//系数矩阵左下非零元素集合，顺序从左到右，从上到下
double b[N];//系数矩阵对角线元素集合，顺序同上
double c[N];//系数矩阵右上非零元素集合，顺序同上
double d[N];//方程组常数项集合
double x[N];//网格节点集合
double e[N];//误差
double maxe_h=0;//最大误差
double *p_1,*p_2,*p_3;//*p_1指向数组d[N],*p_2指向数组e[N],*p_3指向maxe_h double h;//网格步长
int m;//求解区间等分度
int n;//系数矩阵阶数
int i;
printf("请输入求解区间等分度:\nm=");
scanf("%d",&m);
h=3.1415926/m;
n=m-1;
for(i=1;i<m;i++)//系数矩阵和常数项赋值
{
x[i]=i*h;
a[i]=-1;
b[i]=2+h*h;
c[i]=-1;
d[i]=h*h*fi(x[i]);
}
p_1=d;
ThomasMethodsf(a,b,c,p_1,n);
p_2=e;
p_3=&maxe_h;
Errorf(x,d,p_2,p_3,n);
res(d,e,maxe_h,n);
printf("Reslut:\n");
for(i=(n+1)/5;i<=n;i=i+(n+1)/5)
printf("u%d=%f\n",i,d[i]);
printf("\n");
return;
}
double fi(double xi)
{
double fx;
fx=exp(xi)*(sin(xi)-2*cos(xi));
return fx;
}
/*double qi(double xi)
{
}*/
double ui(double xi)
{
double ux;
ux=exp(xi)*sin(xi);
return ux;
}
void ThomasMethodsf(double a[],double b[],double c[],double d[],int n) {
int j;
double p;
// prtf(a,b,c,d,n);
if(fabs(b[1])<1e-15)//判断系数矩阵是否为奇异矩阵函数
{
printf("Error: Matrix A is sigular.");
return;
}
c[1]=c[1]/b[1];//追
d[1]=d[1]/b[1];
for(j=2;j<n;j++)
{
p=b[j]-a[j]*c[j-1];
if(fabs(b[j])<1e-15)//判断系数矩阵是否为奇异矩阵函数
{
printf("Error: Matrix A is sigular.");
return;
}
c[j]=c[j]/p;
d[j]=(d[j]-a[j]*d[j-1])/p;
}
p=b[n]-a[n]*c[n-1];//赶
if(fabs(b[n])<1e-15)//判断系数矩阵是否为奇异矩阵函数
{
printf("Error: Matrix A is sigular.");
return;
}
d[n]=(d[n]-a[n]*d[n-1])/p;
for(j=n-1;j>=1;j--)
d[j]=d[j]-c[j]*d[j+1];
}
void res(double d[],double e[],double maxe_h,int n)
{
FILE *fp;
int j;
if((fp=fopen("Result.txt","a"))==NULL)//打开文件
{
printf("Connot open this file.");
return;
}
fprintf(fp,"m=%d时的数值解：\n",n+1);
for(j=1;j<=n;j=j++)
{
fprintf(fp,"u%d=%f ",j,d[j]);//写入文件
}
fprintf(fp,"\n对应数值解的误差的绝对值：\n");
for(j=1;j<=n;j=j++)
{
fprintf(fp,"%10.3e ",e[j]);//写入文件
}
fprintf(fp,"\n数值解的最大误差：%10.3e\n\n",maxe_h);
fclose(fp);//关闭文件
}
void Errorf(double x[],double d[],double e[],double *maxe_h,int n) {
int i;
for(i=1;i<=n;i++)
{
e[i]=fabs(ui(x[i])-d[i]);
if(e[i]>*maxe_h)
*maxe_h=e[i];
}
}
void prtf(double a[],double b[],double c[],double d[],int n)
{
int i,j;
printf("Coefficient Matrix:\n");
printf("%f %f \n",b[1],c[1]);
printf("%f %f %f\n",a[2],b[2],c[2]);
for(j=3;j<n;j++)
{
for(i=2;i<j;i++)
printf(" ");
printf("%f %f %f\n",a[j],b[j],c[j]);
}
for(i=2;i<n;i++)
printf(" ");
printf("%f %f\n",a[n],b[n]);
printf("Constant:\n");
for(j=1;j<=n;j++) printf("%f ",d[j]);
printf("\n");
}。