亲，问题又来了。。。求帮，程序中出现除零错误额，请教如何避免。。。个人觉得是（追赶法的地方有问题，求助）
#include<stdio.h>#include<stdlib.h>
#include<math.h>
#include "debug.h"
#include "rd_gs.h"
#define N 100
#define M 400
int thomas(double a[N], double b[N+1], double c[N], double d[N+1], double x[N+1][N+1], int n,int j)
/*
* 使用追赶法求解三对角系统
* 输入
* a,b,c: 三对角矩阵的三个对角
* n : 矩阵的维数
* d ：右边项
*
* 输出
* x ：解向量
*/
{
int k;
double tmp;
for(k=0;k<n-1;k++)
{
tmp = a[k+1]/b[k];
b[k+1]=b[k+1]-tmp*c[k];
d[k+1]=d[k+1]-tmp*d[k];
}
x[n-1][j]=d[n-1]/b[n-1];
for(k=n-2;k>=0;k--)
x[k][j]=(d[k]-c[k]*x[k+1][j])/b[k];
return 1;
}
int thomas1(double a[N], double b[N+1], double c[N], double d[N+1], double x[N+1][N+1], int n,int j)
/*
* 使用追赶法求解三对角系统
* 输入
* a,b,c: 三对角矩阵的三个对角
* n : 矩阵的维数
* d ：右边项
*
* 输出
* x ：解向量
*/
{
int k;
double tmp;
for(k=0;k<n-1;k++)
{
tmp = a[k+1]/b[k];
b[k+1]=b[k+1]-tmp*c[k];
d[k+1]=d[k+1]-tmp*d[k];
}
x[j][n-1]=d[n-1]/b[n-1];
for(k=n-2;k>=0;k--)
x[j][k]=(d[k]-c[k]*x[j][k+1])/b[k];
return 1;
}
//封装一下运行的算法函数。
/*int thomas(double a[N], double b[N+1], double c[N], double d[N+1], double x[N+1][N+1], int n,int j)
/*
* 使用追赶法求解三对角系统
* 输入
* a,b,c: 三对角矩阵的三个对角
* n : 矩阵的维数
* d ：右边项
*
* 输出
* x ：解向量
*/
/*{
int k;
double tmp;
double *r,*y;
r=(double *)malloc(sizeof(int)*n);
y=(double *)malloc(sizeof(int)*n);
r[0] = c[0]/b[0];
y[0] = d[0]/b[0];
for(k=1;k<n;k++)
{
r[k]=c[k]/(b[k]-r[k-1]*a[k]);
y[k]=(d[k]-y[k-1]*a[k])/(b[k]-r[k-1]*a[k]);
}
x[n-1][j]=y[n-1];
for(k=n-2;k>=0;k--)
x[k][j]=y[k]-r[k]*x[k+1][j];
return 1;
}
int thomas1(double a[N], double b[N+1], double c[N], double d[N], double x[N+1][N+1], int n,int j)
/*
* 使用追赶法求解三对角系统
* 输入
* a,b,c: 三对角矩阵的三个对角
* n : 矩阵的维数
* d ：右边项
*
* 输出
* x ：解向量
*/
/*{
int k;
double tmp;
double *r,*y;
r=(double *)malloc(sizeof(int)*n);
y=(double *)malloc(sizeof(int)*n);
r[0] = c[0]/b[0];
y[0] = d[0]/b[0];
for(k=0;k<n-1;k++)
{
r[k]=c[k]/(b[k]-r[k-1]*a[k]);
y[k]=(d[k]-y[k-1]*a[k])/(b[k]-r[k-1]*a[k]);
}
x[j][n-1]=d[n-1]/b[n-1];
for(k=n-2;k>=0;k--)
x[j][k]=y[k]-r[k]*x[j][k+1];
return 1;
}*/
int main(void)
{
double T = 4.0,t;
double a = 0.1305;
double b = 0.7695;
double kai = 200.0;
double du = 0.05;
double dv = 1.0;
double dt = T/M;
double h = 1.0/N;
double bs=du*dt/(h*h),bs1=dv*dt/(h*h);
double u[N+1][N+1], v[N+1][N+1],d1[N+1],d2[N+1],a1[N+1][N+1],d[N+1][N+1],s1[N],s2[N+1],s3[N],c1[N],c2[N+1],c3[N];
double x[N+1], y[N+1];
int k,i, j;
char fname[100];
//FILE *fp;
//fp=fopen("u.txt","w");
//产生网格
for(i=0;i<N+1;i++)
x[i] = i*h;
for(j=0;j<N+1;j++)
y[j] = j*h;
//初值条件
for(i=0;i<N+1;i++)
for(j=0;j<N+1;j++){
u[i][j] = a + b + 0.001*exp(-100*(pow(x[i]-1.0/3, 2.0)+pow(y[j]-1.0/2,2.0)));
v[i][j] = b/((a+b)*(a+b));
}
//初始化三对角矩阵
for(i=0;i<N;i++){s1[i]=-bs;c1[i]=-bs1;s3[i]=-bs;c3[i]=-bs1;}
for(i=0;i<N+1;i++){s2[i]=(2*bs+1);c2[i]=(2*bs1+1);}
//算法描述
for(j=1;j<M;j++)
{
printf("%d\n",j);
if(k%1==0){
t = k*dt;
sprintf(fname, "u%.4lf.txt", t);
output_data(u, fname);
sprintf(fname, "v%.4lf.txt", t);
output_data(v, fname);
}
//初始化三对角矩阵2
for(i=0;i<N;i++){s1[i]=-bs;c1[i]=-bs1;s3[i]=-bs;c3[i]=-bs1;}
for(i=0;i<N+1;i++){s2[i]=2*bs+1;c2[i]=2*bs1+1;}
//2n+1时间层
if(j%2==1)
{
for(k=0;k<N+1;k++)
{
if(k==0)
{
for(i=0;i<N+1;i++)
{
a1[i][k]=u[i][k]+dt*kai*(a-u[i][k]+u[i][k]*u[i][k]*v[i][k])+bs*(2*u[i][1]-2*u[i][0]);
d[i][k]=v[i][k]+dt*kai*(b-u[i][k]*u[i][k]*v[i][k])+bs1*(2*v[i][1]-2*v[i][0]);
d1[i]=a1[i][k];
d2[i]=d[i][k];
}
}
else if(k==N)
{
for(i=0;i<N+1;i++)
{
a1[i][k]=u[i][k]+dt*kai*(a-u[i][k]+u[i][k]*u[i][k]*v[i][k])+bs*(u[i][N-1]-2*u[i][N]);
d[i][k]=v[i][k]+dt*kai*(b-u[i][k]*u[i][k]*v[i][k])+bs1*(2*v[i][N-1]-2*v[i][N]);
d1[i]=a1[i][k];
d2[i]=d[i][k];
}
}
else
{
for(i=0;i<N+1;i++)
{
a1[i][k]=u[i][k]+dt*kai*(a-u[i][k]+u[i][k]*u[i][k]*v[i][k])+bs*(u[i][k+1]-2*u[i][k]+u[i][k-1]);
d[i][k]=v[i][k]+dt*kai*(b-u[i][k]*u[i][k]*v[i][k])+bs1*(v[i][k+1]-2*v[i][k]+v[i][k-1]);
d1[i]=a1[i][k];
d2[i]=d[i][k];
}
}
thomas(s1,s2,s3,d1,u,N+1,k);
thomas(c1,c2,c3,d2,v,N+1,k);
}
}
//2n时间层
else
{
for(i=0;i<N+1;i++)
{
if(i==0)
{
for(k=0;k<N+1;k++)
{
a1[i][k]=u[i][k]+dt*kai*(a-u[i][k]+u[i][k]*u[i][k]*v[i][k])+bs*(2*u[1][k]-2*u[0][k]);
d[i][k]=v[i][k]+dt*kai*(b-u[i][k]*u[i][k]*v[i][k])+bs1*(2*v[1][k]-2*v[0][k]);
d1[k]=a1[i][k];
d2[k]=d[i][k];
}
}
else if(i==N)
{
for(k=0;k<N+1;k++)
{
a1[i][k]=u[i][k]+dt*kai*(a-u[i][k]+u[i][k]*u[i][k]*v[i][k])+bs*(2*u[N-1][k]-2*u[N][k]);
d[i][k]=v[i][k]+dt*kai*(b-u[i][k]*u[i][k]*v[i][k])+bs1*(2*v[N-1][k]-2*v[N][k]);
d1[k]=a1[i][k];
d2[k]=d[i][k];
}
}
else
{
for(k=0;k<N+1;k++)
{
a1[i][k]=u[i][k]+dt*kai*(a-u[i][k]+u[i][k]*u[i][k]*v[i][k])+bs*(u[i+1][k]-2*u[i][k]+u[i-1][k]);
d[i][k]=v[i][k]+dt*kai*(b-u[i][k]*u[i][k]*v[i][k])+bs1*(v[i+1][k]-2*v[i][k]+v[i-1][k]);
d1[k]=a1[i][k];
d2[k]=d[i][k];
}
}
}
//thomas1(s1,s2,s3,d1,u,N+1,i);
//thomas1(c1,c2,c3,d2,v,N+1,i);
}
}
//output data
output_data(u, "u.txt");
output_data(v, "v.txt");
return 0;
}
void output_data(double u[N+1][N+1], char *fname)
{
int i, j;
FILE *fd;
fd = OPENfile("w", fname);
for(i=0;i<N+1;i++){
for(j=0;j<N+1;j++)
fprintf(fd, "%.6lf ", u[i][j]);
if(i<N) fprintf(fd, "\n");
}
fclose(fd);
}
FILE *OPENfile( char *mode, char *fname )
{
FILE *fd;
if((fd = fopen(fname, mode)) == (FILE *)NULL)
{
fprintf(stderr,"ERROR: Cannot open file %s in (%s) mode.\n", fname, mode);
exit(-1);
}
return(fd);
}
[ 本帖最后由 zxd675816777 于 2012-7-14 12:21 编辑 ]
好厉害，这这这。。。。。完全看不懂了。。。。。。。