#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;
  }
  

很久不见这位兄弟了。只看了第一个函数，其中归一化处理这块就有问题。估计其它部分若有问题应该也是这个原因。

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];
       }
   }

带下划线的这句
当j < i时没有问题（其实j<i的部分没必要计算）
当j == i时这句的计算会使s[i][i] = 1
当j > i时因为这时s[i][i] == 1，所以这时的归一化处理已经出错了。

你说的运行崩溃怀疑与除0有关。你先修复上面的问题咱们再讨论。

。。。。。。
还有就是，你的算法看起来没用到三对角矩阵的特性，当作普通矩阵处理的。应该利用三对角阵的特性简化算法。



[ 本帖最后由 beyondyf 于 2012-7-10 13:18 编辑 ]

杨大哥啊，谢谢啦！！就是这个问题额，本来想三对角每一个对角都用一个一维数组解决的。。。最近要弄一个反应扩散方程，跟同学一赌气就想都没想直接敲了。。。后面debug的时候那叫一个痛苦额。。。。