今天学懵了。。。写的很乱，程序没错但是运行崩溃额。。。求大神给我敲敲！！只要能正常运行就好。。。
#include<stdio.h>#include<stdlib.h>
#include<math.h>
#define M 4000
#define N 100
//解三对角矩阵的算法1
void sanduijiao1(double s[N+1][N+2],double u[M][N+1][N+1],double bs,int a,int b)
{
int i,j,k;
//追的过程
for(i=0;i<N;i++)
{
for(j=0;j<N+2;j++)
{
s[i][j]=s[i][j]/s[i][i];
}
for(k=0;k<N+2;k++)
{
s[i+1][k]=s[i+1][k]+bs*s[i][k];
}
}
//赶的过程
for(i=N;i>=0;i--)
{
if(i==N)u[a][i][b]=s[i][N+1];
else{
u[a][i][b]=s[i][N+1]-bs/(2*bs+1)*u[a][i+1][b];
}
}
}
//解三对角矩阵的算法2
void sanduijiao2(double s[N+1][N+2],double u[M][N+1][N+1],double bs,int a,int b)
{
int i,j,k;
//追的过程
for(i=0;i<N;i++)
{
for(j=0;j<N+2;j++)
{
s[i][j]=s[i][j]/s[i][i];
}
for(k=0;k<N+2;k++)
{
s[i+1][k]=s[i+1][k]+bs*s[i][k];
}
}
//赶的过程
for(i=N;i>=0;i--)
{
if(i==N)u[a][b][i]=s[i][N+1];
else{
u[a][b][i]=s[i][N+1]-bs/(2*bs+1)*u[a][b][i+1];
}
}
}
int main(void)
{
FILE *fp;
fp=fopen("u.txt","w");
int i,j,k;
double u[M][N+1][N+1],v[M][N+1][N+1],a[M][N+1][N+1],d[M][N+1][N+1],x[N+1],y[N+1];
double a1=0.1305,b=0.7695,du=0.05,dv=1.0,T=4.0,dt=T/M,h=1.0/N,bs=du*dt/(h*h),bs1=dv*dt/(h*h),kai=200;
//三对角矩阵
double s[N+1][N+2],s1[N+1][N+2];
//产生网格
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[0][i][j] = a1 + b + 0.001*exp(-100*(pow(x[i]-1.0/3, 2.0)+pow(y[j]-1.0/2,2.0)));
v[0][i][j] = b/((a1+b)*(a1+b));
}
//初始三对角矩阵
for(i=0;i<N+1;i++)
for(j=0;j<N+2;j++)
{
if(i==j)s[i][j]=2*bs+1;
else if(j==i+1&&j==i-1)s[i][j]=-bs;
else s[i][j]=0.0;
}
for(i=0;i<N+1;i++)
for(j=0;j<N+2;j++)
{
if(i==j)s1[i][j]=2*bs+1;
else if(j==i+1&&j==i-1)s1[i][j]=-bs;
else s1[i][j]=0.0;
}
//主算法程序部分
for(j=0;j<M;j++)
{
printf("%d",j);
if(j%100==0)fprintf(fp,"%lf",u[j][N+1][N+1]);
if(j%2==1)
{
for(k=0;k<N+1;k++)
{
for(i=0;i<N+1;i++)
{
a[j][i][k]=u[j-1][i][k]+kai*(a1-u[j-1][i][k]+u[j-1][i][k]*u[j-1][i][k]*v[j-1][i][k])+dt/(h*h)*(u[j-1][i][k+1]-2*u[j-1][i][k]+u[j-1][i][k-1]);
d[j][i][k]=v[j-1][i][k]+kai*(b-u[j-1][i][k]*u[j-1][i][k]*v[j-1][i][k])+dt/(h*h)*(v[j-1][i][k+1]-2*v[j-1][i][k]+v[j-1][i][k-1]);
s[i][N+1]=a[j][i][k];
s1[i][N+1]=d[j][i][k];
}
sanduijiao1(s,u,bs,j,k);
sanduijiao1(s1,v,bs1,j,k);
}
}
else
{
for(i=0;i<N+1;i++)
{
for(k=0;k<N+1;k++)
{
a[j][i][k]=u[j-1][i][k]+kai*(a1-u[j-1][i][k]+u[j-1][i][k]*u[j-1][i][k]*v[j-1][i][k])+dt/(h*h)*(u[j-1][i][k+1]-2*u[j-1][i][k]+u[j-1][i][k-1]);
d[j][i][k]=v[j-1][i][k]+kai*(b-u[j-1][i][k]*u[j-1][i][k]*v[j-1][i][k])+dt/(h*h)*(v[j-1][i][k+1]-2*v[j-1][i][k]+v[j-1][i][k-1]);
s[k][N+1]=a[j][i][k];
s1[k][N+1]=a[j][i][k];
}
sanduijiao2(s,u,bs,j,i);
sanduijiao2(s1,v,bs1,j,i);
}
}
}
return 0;
}
