本人做了一个二维对称翼型数值解的算法，但方程组的解太糟糕了，求知道原因以及修正方法。
程序如下：
//-------------二维对称翼型数值解----------------//
#include<stdio.h>

class Airfoil
{
public:
    //----------构造函数----------//
    void setAirfoil(int dimension,double step1,double V)
    {
        n=dimension;
        step=step1;
        Velocity=V;
    }//end

     //-----------辅助函数-------------//
    //the transfer function
    void tansf(double *a,double *b)
    {
       double c;
       c=*a;
       *a=*b;
       *b=c;
    }//end the transf function

    void printfx(double x[])
    {
        int i;
        printf("\n\n");
        for(i=0;i<n;i++)
            printf("%0.3lf   ",x[i]);
        printf("\n\n");
    }//end

    void print(double *a)
    {
        printf("\n\n");
        int i,j;
        for(i=0;i<n;i++)
        {   for(j=0;j<n;j++)
             printf("%0.3lf  ",a[i*n+j]);
             printf("\n");
        }
        printf("\n\n");
    }


    //the max function
    int max(double *a,int k)
    {
       int i,j;
       double m;
       j=k;
       m=a[k*n+k];
       for(i=k+1;i<n;i++)
       {
          if(m<a[i*n+k])
          {  m=a[i*n+k];
             j=i;
          }
       }
       return(j);
    }//end the max function
    ////////////////////////////////////////////////////////

     //----------------主要函数--------------------------//
     //Column main element gauss elimination method
    void gauss(double *a,double *b,double *x)
    {
        int i,j,k,p;
        double m;

        //System of linear equations gaussian elimination process
        for(k=0;k<n-1;k++)
    {   

            p=max(a,k);

            //---------行变换------------//
            if(p!=k)
            {
                for(i=k;i<n;i++)
                   tansf(&a[k*n+i],&a[p*n+i]);
                   tansf(&b[k],&b[p]);
            }
            //---------常数项交换----------//

            
            //---------消去过程-----------//
            for(i=k+1;i<n;i++)
            {
                m=a[i*n+k]/a[k*n+k];
                for(j=k;j<n;j++)
                    a[i*n+j]=a[i*n+j]-m*a[k*n+j];
                b[i]=b[i]-m*b[k];
            }
            //--------------------------------//
        
        }

        //-----------回代过程-------------//
        x[n-1]=b[n-1]/a[n*n-1];
        for(k=n-2;k>=0;k--)
        {
            m=0.0;
            for(i=k+1;i<n;i++)
                m=m+a[k*n+i]*x[i];
            x[k]=(b[k]-m)/a[k*n+k];
        }
        //----------------------------------//

    }//end the gauss function

    //-------------生成方程组系数------------//
    void setCB(double x[],double y[],double *a)
    {
        int i,j;
        for(i=0;i<n;i++)
        {   
            for(j=0;j<n;j++)
            {   a[i*n+j]=y[i]*step/((x[i]-step*(j+0.5))*(x[i]-step*(j+0.5))+y[i]*y[i]);
                a[i*n+j]=a[i*n+j]/Velocity;
            }
        }
    }//end

    //---------------求解偶极子强度------------//
    void solve(double *a,double y[],double m[])
    {
        gauss(a,y,m);
    }//end

    //--------------流函数---------------//
    double stream(double x2,double y2,double m[])
    {
        int i;
        double z=0;

        for(i=0;i<n;i++)
            z=z+m[i]*step*y2/((x2-step*(i+0.5))*(x2-step*(i+0.5))+y2*y2);
        return(Velocity*y2-z);
    }//end

public:
    int n;
    double step;
    double Velocity;
};

void main()
{
    static double x[25]={0.500,0.750,1.250,2.500,5.000,7.500,10.000,15.000,20.000,25.000,
    30.000,35.000,40.000,45.000,50.000,55.000,60.000,65.000,70.000,75.000,80.000,85.000,90.000,95.000,97.50};
    static double y[25]={0.250,0.350,0.535,0.930,1.580,2.120,2.585,3.365,3.980,4.475,4.860,
    5.150,5.355,5.475,5.515,5.475,5.355,5.150,4.860,4.475,3.980,3.365,2.585,1.580,0.780};

    //static double x[5]={10.0,35.0,90.0};
    //static double y[5]={2.585,5.150,2.585};
    double a[25][25];
    double m[25];

    Airfoil myAirfoil;

    myAirfoil.setAirfoil(25,4.0,10.0);

    //--------测试系统值----------//
    printf("%d  \n\n",myAirfoil.n);
    printf("%0.6lf  \n\n",myAirfoil.step);
    printf("%0.6lf  \n\n",myAirfoil.Velocity);

    myAirfoil.setCB(x,y,*a);
    //myAirfoil.print(*a);

    myAirfoil.solve(*a,y,m);



    myAirfoil.printfx(m);
    printf("%0.6lf  \n\n",myAirfoil.stream(35.0,35.0,m));
}