app.c

数值分析中用Gauss消去法求非线性方程组的解

/*数值分析中求非线性方程组的解的*/
#include "stdio.h"
#include "math.h"

double max4();      /*求四个数中的最大值*/
double Gauss();     /*Gauss消去法求非线性方程组*/
double Interpolation();     /*作插值*/
double Inverse();       /*求逆矩阵函数*/
double Fitting();       /*拟和函数*/

void main(void)
{
    int i,j,l,k;
    double m,x,y,z=0,a[2],u[11][21],t[11][21];
    double B[11][21]={0},b[21][21],b8[21][21]={0},G[21][21]={0},g[21][21],g8[21][21]={0},c[21][21]={0};
    for(i=0;i<=10;i++){
        for(j=0;j<=20;j++){
            x=0.08*i,y=0.5+0.05*j;
            a[0]=x,a[1]=y;
            Gauss(a);
            u[i][j]=a[0],t[i][j]=a[1];
        }
    }
    for(i=0;i<=10;i++){
        for(j=0;j<=20;j++)
            u[i][j]=Interpolation(t[i][j],u[i][j]);
    }
    k=0;
    do{
        for(i=0;i<=20;i++){
            for(j=0;j<=20;j++){
                b[i][j]=0;
                g[i][j]=0;
                b8[i][j]=0;
                g8[i][j]=0;
                c[j][j]=0;
            }
        }
        for(i=0;i<=10;i++){
            B[i][0]=1,x=0.08*i;
            for(j=1;j<=k;j++)
                B[i][j]=pow(x,j);
        }
        for(j=0;j<=20;j++){
            G[j][0]=1,y=0.5+0.05*j;
            for(i=1;i<=k;i++)
                G[j][i]=pow(y,i);
        }
        for(i=0;i<=k;i++){
            for(j=0;j<=k;j++){
                  for(l=0;l<=10;l++)
                      b[i][j]=b[i][j]+B[l][i]*B[l][j];
                  for(l=0;l<=20;l++)
                     g[i][j]=g[i][j]+G[l][i]*G[l][j];
            }
        }
        Inverse(b,k);
        Inverse(g,k);
        for(i=0;i<=k;i++){
            for(j=0;j<=10;j++){
                  z=0;
                for(l=0;l<=k;l++)
                    z=z+b[i][l]*B[j][l];
                b8[i][j]=z;
            }
        }
        for(i=0;i<=20;i++){
            for(j=0;j<=k;j++){
                  z=0;
                  for(l=0;l<=k;l++)
                    z=z+G[i][l]*g[l][j];
                  g8[i][j]=z;
            }
         }
         for(i=0;i<=k;i++){
            for(j=0;j<=20;j++){
                   z=0;
                   for(l=0;l<=10;l++)
                     z=z+b8[i][l]*u[l][j];
                   b[i][j]=z;
            }
        }
        for(i=0;i<=k;i++){
            for(j=0;j<=k;j++){
                z=0;
                for(l=0;l<=20;l++)
                     z=z+b[i][l]*g8[l][j];
                c[i][j]=z;
            }
        }
        z=0;
        for(i=0;i<=10;i++){
            for(j=0;j<=20;j++)
                 z=z+(Fitting(c,k,B[i][1],G[j][1])-u[i][j])*(Fitting(c,k,B[i][1],G[j][1])-u[i][j]);
        }
        printf(" k=%d",k);
        printf(" %18.11e\n",z);
        k++;
    }while(z>1e-7);
    getch();
}

/*求四个值中的最大值*/
double max4(double x1,double x2,double x3,double x4)
{
    double z1,z2,z;
    if(x1>x2)z1=x1;
    else z1=x2;
    if(x3>x4)z2=x3;
    else z2=x4;
    if(z1>z2)z=z1;
    else z=z2;
    return z;
}

/*使用列主元的Gauss消去法求解非线性方程组*/
double Gauss(double x[2])
{
    int i,j,k,n;
    double t,u,v,w,f1,m,mi,T=1,U=1,V=1,W=1;
    double a[5][6],b[5][6],c[2];
    for(i=1;i<=4;i++){
        for(j=1;j<=5;j++)
            a[i][j]=1;
    }
    do{
        a[1][1]=-0.5*sin(T),a[1][5]=-(0.5*cos(T)+U+V+W-x[0]-2.67),
        a[2][2]=0.5*cos(U), a[2][5]=-(T+0.5*sin(U)+V+W-x[1]-1.07),
        a[3][1]=0.5, a[3][3]=-sin(V), a[3][5]=-(0.5*T+U+cos(V)+W-x[0]-3.74),
        a[4][2]=0.5, a[4][4]=cos(W),  a[4][5]=-(T+0.5*U+V+sin(W)-x[1]-0.79);
        for(i=1;i<=4;i++){
            for(j=1;j<=5;j++)
                b[i][j]=a[i][j];
        }
        for(k=1;k<=3;k++){
            n=k;
            f1=fabs(b[k][k]);
            for(i=k+1;i<4;i++){             /*选取主元*/
                m=fabs(b[i][k]);
                if(m<=f1)continue;
                n=i;
                f1=m;
            }
            for(j=k;j<=5;j++){
                f1=b[k][j];
                b[k][j]=b[n][j];
                b[n][j]=f1;
            }
            for(i=k+1;i<=4;i++){
                mi=b[i][k]/b[k][k];
                for(j=k+1;j<=5;j++)
                    b[i][j]=b[i][j]-mi*b[k][j];
            }
        }
        w=b[4][5]/b[4][4];
        v=(b[3][5]-b[3][4]*w)/b[3][3];
        u=(b[2][5]-b[2][3]*v-b[2][4]*w)/b[2][2];
        t=(b[1][5]-b[1][2]*u-b[1][3]*v-b[1][4]*w)/b[1][1];
        W=W+w;
        V=V+v;
        U=U+u;
        T=T+t;
        m=max4(fabs(w),fabs(v),fabs(u),fabs(t))/max4(fabs(W),fabs(V),fabs(U),fabs(T));
    }while (m>=1e-12);
    x[0]=U,x[1]=T;
    return 1;
}

/*插值函数*/
double Interpolation(double x,double y)
{   
    int i,j,k,n,a,t;
    float m[6]={0,0.2,0.4,0.6,0.8,1.0},b[6]={0,0.4,0.8,1.2,1.6,2.0};
    float z[6][6]={ {-0.5,-0.34,0.14,0.94,2.06,3.5},{-0.42,-0.5,-0.26,0.3,1.18,2.38},
                {-0.18,-0.5,-0.5,-0.18,0.46,1.42},{0.22,-0.34,-0.58,-0.5,-0.1,0.62},
                {0.78,-0.02,-0.5,-0.66,-0.5,-0.02},{1.5,0.46,-0.26,-0.66,-0.74,-0.5}};
    double Lt[6]={0},Lu[6]={0},p;
    if (x<=0.3) n=1;
    else if(x>0.7) n=4;
    else{    
        for(i=2;i<=3;i++){
            if(x>(m[i]-0.1)&&x<=(m[i]+0.1))
                n=i;
            else continue;
        }
    }
    for(k=n-1;k<=n+1;k++){
        Lt[k]=1;
        for (t=n-1;t<=n+1;t++){
            if(k==t) continue;
            Lt[k]=Lt[k]*(x-m[t])/(m[k]-m[t]);
        }
    }

    if (y<=0.6) a=1;
    else if(y>1.4) a=4;
    else{
        for(i=2;i<=3;i++){
            if(y>(b[i]-0.2)&&y<=(b[i]+0.2))
                a=i;
            else continue;
        }
    }
    for(k=a-1;k<=a+1;k++){
        Lu[k]=1;
        for (t=a-1;t<=a+1;t++){
            if(k==t) continue;
            Lu[k]=Lu[k]*(y-b[t])/(b[k]-b[t]);
        }
    }
    p=0;
    for(i=n-1;i<=n+1;i++){
        for(j=a-1;j<=a+1;j++){
            p=p+Lt[i]*Lu[j]*z[i][j];
        }
    }
    return p;
}

/*Guass-Jordan法求矩阵的逆*/
double Inverse(double A[21][21],int n)
{
    int i,j,k,n1;
    double m,temp,f1,B[21][42];
    for(i=0;i<=20;i++){
        for(j=0;j<=41;j++)
            B[i][j]=0;
    }
    for(i=0;i<=n;i++){
        for(j=0;j<=n;j++)
            B[i][j]=A[i][j];
        B[i][i+n+1]=1;
    }
    for(k=0;k<=n;k++){
        n1=k;
        f1=fabs(B[k][k]);
        for(i=k+1;i<=n;i++){
            m=fabs(B[i][k]);
            if(m<=f1)continue;
            n1=i;
            f1=m;
        }
        for(j=0;j<=2*n+1;j++){
            f1=B[k][j];
            B[k][j]=B[n1][j];
            B[n1][j]=f1; 
        }
        temp=B[k][k];
        for (j=k;j<=2*n+1;j++) 
            B[k][j]/=temp;
        for (i=0;i<=n;i++){
            if(i==k) continue;
            temp=B[i][k];
            for (j=k;j<=2*n+1;j++) 
                B[i][j]-=temp*B[k][j];
        }
    }
    for(i=0;i<=n;i++){
        for(j=0;j<=n;j++)
            A[i][j]=B[i][j+n+1];
    }
    return 1;
}

/*求拟合值*/
double Fitting(double a[21][21],int k,double x,double y)
{
    int i,j;
    double z=0;
    for(i=0;i<=k;i++){
        if(i==0){
            for(j=0;j<=k;j++){
                if(j==0)
                    z=z+a[i][j]*1*1;
                else
                    z=z+a[i][j]*1*pow(y,j);
            }
        }
        else{
            for(j=0;j<=k;j++){
                if(j==0)
                    z=z+a[i][j]*pow(x,i)*1;
                else
                    z=z+a[i][j]*pow(x,i)*pow(y,j);
            }
        }
    }
    return z;
}