牛顿法求解非线性方程组，插值，北航数值分析第三次大作业

几年前写的代码，我现在把它贴出来，用C++编写的。

利用牛顿法求解非线性方程组

#include "stdio.h"#include "math.h"void newton(double x[11],double y[11]);//牛顿法求解非线性方程组void eryuanchazhi(double te[11][21],double ue[11][21]);//分片二次代数插值，求z=f(x,y)void qumianchazhi();    //曲面插值，求p(x,y)void compare();//比较f(x*,y*)与p(x*,y*)void inverse(double X[10][10],int N);//求逆矩阵double fanshu(double *p);//求向量的无穷范数double te[11][21]={0};//存储非线线方程组的解tdouble ue[11][21]={0};//存储非线线方程组的解udouble ze[11][21]={0};//存储分片二次插值的解zdouble x[11]={0};//xdouble y[21]={0};//ydouble inv[10][10]={0};//存储某矩阵对应的逆矩阵double  C[10][10]={0};//存储曲面插值的系数矩阵double rs=0;//存储k值，用于比较函数中使用void main(){int i,j;for(i=0;i<11;i++)for(j=0;j<21;j++){x[i]=0.08*i;y[j]=0.5+0.05*j;}newton(x,y);eryuanchazhi(te,ue);qumianchazhi();compare();}void newton(double x[11],double y[11]){double t=0;double u=0;double w=0;double v=0;double dF[5][5]={0};double F[5]={0};double d[5]={0};int i,j,k;int ik;int cnt,count;double temp=0;double m=0;double jie[4]={0};for(i=0;i<11;i++)for(j=0;j<21;j++){t=1;u=1;w=1;v=1;//初始化//===============求解非线性方程组=============//do{ dF[1][1]=-0.5*sin(t);  dF[1][2]=1;  dF[1][3]=1;  dF[1][4]=1;  dF[2][1]=1;  dF[2][2]=0.5*cos(u);  dF[2][3]=1;  dF[2][4]=1;  dF[3][1]=0.5;  dF[3][2]=1;  dF[3][3]=-sin(v);  dF[3][4]=1;  dF[4][1]=1;  dF[4][2]=0.5;  dF[4][3]=1;  dF[4][4]=cos(w);  F[1]=-(0.5*cos(t)+u+v+w-x[i]-2.67);  F[2]=-(t+0.5*sin(u)+v+w-y[j]-1.07);  F[3]=-(0.5*t+u+cos(v)+w-x[i]-3.74);  F[4]=-(t+0.5*u+v+sin(w)-y[j]-0.79); //对Newton法的系数矩阵进行赋值,注意这里是-F()//=========================列主元高斯消元法求解方程组====================//for (k=1;k<=3;k++){ik=k;for(cnt=k;cnt<=4;cnt++){if (fabs(dF[cnt][k])>fabs(dF[ik][k])){ik=cnt;}}temp=F[k];F[k]=F[ik];F[ik]=temp;for (cnt=k;cnt<=4;cnt++){temp=dF[k][cnt];dF[k][cnt]=dF[ik][cnt];dF[ik][cnt]=temp;}for (cnt=k+1;cnt<=4;cnt++){m=dF[cnt][k]/dF[k][k];for (count=k+1;count<=4;count++)dF[cnt][count]=dF[cnt][count]-m*dF[k][count];F[cnt]=F[cnt]-m*F[k];}}d[4]=F[4]/dF[4][4];for (k=4-1;k>=1;k--){temp=0;for (cnt=k+1;cnt<=4;cnt++)temp+=dF[k][cnt]*d[cnt];d[k]=(F[k]-temp)/dF[k][k];}  t=t+d[1];  u=u+d[2];  v=v+d[3];  w=w+d[4];  te[i][j]=t;  ue[i][j]=u;  jie[1]=t;  jie[2]=u;  jie[3]=v;  jie[4]=w;}while(fanshu(d)/fanshu(jie)>1e-12);printf("te[%d][%d]=%13.11e\t",i,j,te[i][j]);printf("ue[%d][%d]=%13.11e\n",i,j,ue[i][j]);}}void eryuanchazhi(double te[11][21],double ue[11][21]){int i,j,ii,jj,m;int r[11][21]={0};int k[11][21]={0};double temp=0;double u[6]={0,0.4,0.8,1.2,1.6,2};double t[6]={0,0.2,0.4,0.6,0.8,1.0};double 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}};for(i=0;i<=10;i++)for(j=0;j<=20;j++){if(te[i][j]<=0.3)r[i][j]=1;else if(te[i][j]>=0.7)r[i][j]=4;else{ for(m=2;m<=3;m++)if((te[i][j]>0.2*m-0.1)&&(te[i][j]<0.2*m+0.1))r[i][j]=m;}  //t的插值点if(ue[i][j]<=0.6)k[i][j]=1;else if(ue[i][j]>=1.4)k[i][j]=4;else{ for(m=2;m<=3;m++)if((ue[i][j]>0.4*m-0.2)&&(ue[i][j]<0.4*m+0.2))k[i][j]=m;}    //u的插值点ze[i][j]=0;for(ii=k[i][j]-1;ii<=k[i][j]+1;ii++)for(jj=r[i][j]-1;jj<=r[i][j]+1;jj++){temp=1;for(m=k[i][j]-1;m<=k[i][j]+1;m++)if(m!=ii)temp*=(te[i][j]-t[m])/(t[ii]-t[m]);for(m=r[i][j]-1;m<=r[i][j]+1;m++)if(m!=jj)temp*=(ue[i][j]-u[m])/(u[jj]-u[m]);ze[i][j]+=temp*z[ii][jj];}}FILE *fp=fopen("shubiao.txt","a");for(i=0;i<=10;i++)for(j=0;j<=20;j++)fprintf(fp,"x[%d]=%f\ty[%d]=%f\tf(x,y)=%13.11e\n",i,x[i],j,y[j],ze[i][j]);fclose(fp);}void qumianchazhi(){int i,j;int ii,jj;int k,t;double temp=0;double temp1=0;double sigma=0;doubleB[11+1][10]={0};doubleG[21+1][10]={0};double  BB[10][10]={0};double  GG[10][10]={0};double  BBN[10][10]={0};double  GGN[10][10]={0};double  BBB[10][12]={0};double  BBBU[10][22]={0};double  BBBUG[10][10]={0};double p=0;FILE *fp;fp=fopen("xishu.txt","a");k=0+1;//为了以后矩阵运算方便，这里先加上1do{for(i=1;i<=11;i++)for(j=1;j<=k;j++){B[i][j]=pow(x[i-1],j-1);}//B阵赋值for(i=1;i<=21;i++)for(j=1;j<=k;j++){G[i][j]=pow(y[i-1],j-1);}   //G阵赋值for(i=1;i<=k;i++)for(j=1;j<=k;j++){temp=0;for(t=1;t<=11;t++)temp+=B[t][i]*B[t][j];BB[i][j]=temp;}//Bt*Bfor(i=1;i<=k;i++)for(j=1;j<=k;j++){temp=0;for(t=1;t<=21;t++)temp+=G[t][i]*G[t][j]; GG[i][j]=temp;}//Gt*Ginverse(BB,k);for(i=1;i<=k;i++)for(j=1;j<=k;j++)BBN[i][j]=inv[i][j];//Bt*B的逆inverse(GG,k);for(i=1;i<=k;i++)for(j=1;j<=k;j++)GGN[i][j]=inv[i][j];//Gt*G的逆for(i=1;i<=k;i++)for(j=1;j<=11;j++){temp=0;for(t=1;t<=k;t++)temp+=BBN[i][t]*B[j][t];BBB[i][j]=temp;}//Bt*B的逆*Btfor(i=1;i<=k;i++)for(j=1;j<=21;j++){temp=0;for(t=1;t<=11;t++)temp+=BBB[i][t]*ze[t-1][j-1];BBBU[i][j]=temp;}//再乘以Ufor(i=1;i<=k;i++)for(j=1;j<=k;j++){temp=0;for(t=1;t<=21;t++)temp+=BBBU[i][t]*G[t][j];BBBUG[i][j]=temp;}//再乘以Gfor(i=1;i<=k;i++)for(j=1;j<=k;j++){temp=0;for(t=1;t<=k;t++)temp+=BBBUG[i][t]*GGN[t][j];C[i][j]=temp;}//最后乘以Gt*G的逆得到系数矩阵for(i=0;i<11;i++)for(j=0;j<21;j++){temp=0;for(ii=0;ii<k;ii++)for(jj=0;jj<k;jj++)temp+=C[ii+1][jj+1]*pow(x[i],ii)*pow(y[j],jj);//这里的下标，矩阵下标与数组下标差1p=temp;temp1+=(ze[i][j]-p)*(ze[i][j]-p);}sigma=temp1;temp1=0;fprintf(fp,"k=%d,sigma=%13.11e\n",k-1,sigma);k++;}while (sigma>=1e-7);k--;//减去最后一次判断循环条件前加的1rs=k;fprintf(fp,"k=%d\n",k-1);//减去一开始为了矩阵运算加的1for(i=1;i<=k;i++)for(j=1;j<=k;j++)fprintf(fp,"C[%d][%d]=%13.11e\n",i-1,j-1,C[i][j]);//这里为了匹配矩阵下标和书上所示的系数下标所以减1}void compare(){double pn[8+1][5+1]={0};//下标都加1，方便以后运算int i,j;double temp=0;int ii,jj;for(i=1;i<=8;i++)x[i]=0.1*i;for(j=1;j<=5;j++)y[j]=0.5+0.2*j;for(i=1;i<=8;i++)for(j=1;j<=5;j++){temp=0;for(ii=0;ii<=rs;ii++)for(jj=0;jj<=rs;jj++)temp+=C[ii+1][jj+1]*pow(x[i],ii)*pow(y[j],jj);pn[i][j]=temp;}newton(x,y);eryuanchazhi(te,ue);FILE *fp;fp=fopen("shubiao1.txt","a");for(i=1;i<=8;i++)for(j=1;j<=5;j++){fprintf(fp,"x[%d]=%f\t y[%d]=%f\t \n",i,x[i],j,y[j]); fprintf(fp,"z[%d][%d]=%13.11e\t pn[%d][%d]=%13.11e\t\n",i,j,ze[i][j],i,j,pn[i][j]);fprintf(fp,"两者差值sigma=%13.11e\n",ze[i][j]-pn[i][j]);}}double fanshu(double *p){ int i=0;double max=fabs(p[1]);for(i=1;i<=4;i++){ if(fabs(p[i])>max)max=fabs(p[i]);}return(max);}void inverse(double X[10][10],int N){double matrix[10][10];double I[10][10]={0};int i,j;int k,ik,cnt,count;double m;double temp;for(i=1;i<=N;i++)for(j=1;j<=N;j++)matrix[i][j]=X[i][j];for(i=1;i<=N;i++)I[i][i]=1;//单位阵//================列主元的高斯消元法求逆矩阵======================//for (k=1;k<=N-1;k++){ik=k;for(cnt=k;cnt<=N;cnt++){if (fabs(matrix[cnt][k])>fabs(matrix[ik][k])){ik=cnt;}}for (cnt=1;cnt<=N;cnt++){temp=matrix[k][cnt];matrix[k][cnt]=matrix[ik][cnt];matrix[ik][cnt]=temp;temp=I[k][cnt];I[k][cnt]=I[ik][cnt];I[ik][cnt]=temp;}for (cnt=k+1;cnt<=N;cnt++){m=matrix[cnt][k]/matrix[k][k];for (count=1;count<=N;count++){matrix[cnt][count]=matrix[cnt][count]-m*matrix[k][count];I[cnt][count]=I[cnt][count]-m*I[k][count];}}}for(i=1;i<=N;i++){inv[N][i]=I[N][i]/matrix[N][N];for (k=N-1;k>=1;k--){temp=0;for (cnt=k+1;cnt<=N;cnt++)temp+=matrix[k][cnt]*inv[cnt][i];inv[k][i]=(I[k][i]-temp)/matrix[k][k];}}}