数值分析上机作业3——求解非线性方程组以及二元函数的插值拟合1. 算法设计对于全部的插值节点(,),0,1,...,10,0,1,...,20i j x y i j ==，带入非线性方程组中，用Newton 迭代法解非线性方程组，得到(,),0,1,...,10,0,1,...,20i j t u i j ==。

对(,)i j t u ，在二维数表中进行插值，采用分片双二次插值法。

插值过程中，先选择分片区域的中心节点，在数表中的列记为(0:5)tt ，行记为(0:5)uu ，中心节点记为(,)a b ，生成向量_(0:2)t temp ，_(0)(())((1))/(((1)())((1)(1)))i i t temp t tt a t tt a tt a tt a tt a tt a =--+----+， _(1)((1))((1))/((()(1))(()(1)))i i t temp t tt a t tt a tt a tt a tt a tt a =---+---+， _(2)((1))(())/(((1)(1))((1)()))i i t temp t tt a t tt a tt a tt a tt a tt a =---+--+-，同理，生成向量_(0:2)u temp ，_(0)(())((1))/(((1)())((1)(1)))_(1)((1))((1))/((()(1))(()(1)))_(2)((1))(())/(((1)(1))((1)())j j j j j j u temp u uu a u uu a uu a uu a uu a uu a u temp u uu a u uu a uu a uu a uu a uu a u temp u uu a u uu a uu a uu a uu a uu a =--+----+=---+---+=---+--+-)记数表中以分片区域中心节点为中心的3×3的矩阵为T , 对于(,)i j t u 插值结果为(_)()(_)T t temp T u temp 。

在拟合,0(,)kr s rsr s p x y cx y ==∑时，需要计算11()()T T T C B B B UG G G --=，令,T T W B B M G G ==，计算,W M 时，根据对称性只需要计算对角线元素和对角线以上元素即可，节省运算时间。

于是TWCM B UG E ==，用选主元的LU 分解法求解WF E =，再计算TTMC F =，这里,C F 只需要按行取元素进行运算即可，故不需要进行转置运算。

k 从1到9依次增加，计算σ的值，当710σ-≤时，得到达到精度的最小的k 。

2.打印输出结果f(x,y):4.46504018480e-001 3.24683262927e-001 2.10159686683e-001 1.03043608316e-0013.40189556266e-003 -8.87358136384e-0026.38015226510e-001 5.0661*******e-001 3.82176369277e-001 2.64863491154e-0011.54780200285e-001 5.199********e-0028.40081395765e-001 6.99764165673e-001 5.66061442351e-001 4.39171608118e-0013.19242138041e-001 2.06376192387e-0011.0515*******e+000 9.029********e-001 7.60580266859e-001 6.24715198145e-0014.95519756001e-001 3.73134042774e-0011.27124675148e+000 1.11500201815e+000 9.64607727215e-001 8.20347369475e-0016.82447678179e-001 5.51085208597e-0011.49832105248e+000 1.33499863207e+000 1.177********e+000 1.025********e+0008.78960023174e-001 7.39145108703e-0011.73189274038e+000 1.56203457721e+000 1.39721691821e+000 1.23780100674e+0001.0840*******e+000 9.36322772315e-0011.97122178640e+000 1.79532959950e+000 1.62406711323e+000 1.45783058271e+0001.29695464975e+000 1.14171810545e+000k=1, sigma=3.22090897363e+000k=2, sigma=4.65996003320e-003k=3, sigma=1.72117537926e-004k=4, sigma=3.30953430190e-006k=5, sigma=2.54137773513e-008C_rs:2.021********e+000 -3.66842591519e+000 7.09246688452e-001 8.48607444215e-001-4.158********e-001 6.74322022261e-0023.19192646165e+000 -7.41209812962e-001 -2.69690653026e+000 1.63095275395e+000-4.84601977266e-001 6.05908514261e-0022.56706343343e-001 1.58096413206e+000 -4.65701259544e-001 -7.89186706497e-0021.00853116187e-001 -2.0768*******e-002-2.68608304872e-001 -7.33963449843e-001 1.08429601112e+000 -8.156********e-0013.0728*******e-001 -4.68489486618e-0022.16521800059e-001 -1.73026852965e-001 -8.41324310602e-002 2.55736987891e-001-1.47683427939e-001 2.77711894906e-002-5.54328606191e-002 1.40518220408e-001 -1.30388672239e-001 3.44966421960e-0026.95959872154e-003 -3.30022941240e-003f(x*,y*):0.194720 -0.183037 -0.445498 -0.597567 -0.646460 0.405979 -0.022516 -0.338221 -0.544438 -0.647361 0.634777 0.158801 -0.207366 -0.465358 -0.6202710.878960 0.358651 -0.055253 -0.362680 -0.5675651.136611 0.574980 0.115992 -0.238568 -0.491434 1.406042 0.805941 0.304429 -0.095016 -0.393902 1.685784 1.049881 0.508294 0.066149 -0.276834 1.974575 1.305335 0.725992 0.243254 -0.141950p(x*,y*):0.194730 -0.183042 -0.445500 -0.597559 -0.646446 0.405990 -0.022521 -0.338224 -0.544430 -0.647348 0.634787 0.158796 -0.207369 -0.465350 -0.6202570.878970 0.358646 -0.055255 -0.362671 -0.5675511.136620 0.574976 0.115989 -0.238560 -0.491421 1.406051 0.805937 0.304426 -0.095009 -0.393890 1.685791 1.049878 0.508291 0.066156 -0.276822 1.974581 1.305332 0.725989 0.243261 -0.141939f x y插值的数值情况：对(,)(,)(,)i i i i f x y p x y 的数值情况（z 坐标为10-7）：3. 源程序 1) 主函数#include<stdio.h> #include<math.h>#define length_x 11//x 向量长度 #define length_y 21//y 向量长度 #define x_length 8//x*向量长度 #define y_length 5//y*向量长度extern double interpolation(double tt,double uu,double **table); //输入一点，输出一点插值extern double row_multi_sum(double **B,int i,int j,int length); //矩阵列相乘int main() {double tt=0; double uu=0;double table[6][6]={0}; double ff=0;double xx[length_x]={0};double yy[length_y]={0};double sum=0.0;double U[length_x][length_y]={0};double B[length_x][10]={0};double G[length_y][10]={0};double W[10][10]={0};double T[10][10]={0};double M[10][10]={0};double C[10][10]={0};double F[10][length_y]={0};double E[10][10]={0};double XX[10][10]={0};//double ff[length_x][length_y]={0}; //记录f阵double pp[length_x][length_y]={0}; //记录p阵int r=0;int s=0;int ii,jj,kk;int k_max=0;FILE *fd;fd = fopen("tmptmp.txt", "w");ii=0;jj=0;kk=0;table[0][0]=-0.5;table[0][1]=-0.34;table[0][2]=0.14;table[0][3]=0.94;table[0][4]=2.06;table[0][5]=3.5;table[1][0]=-0.42;table[1][1]=-0.5;table[1][2]=-0.26;table[1][3]=0.3;table[1][4]=1.18;table[1][5]=2.38;table[2][0]=-0.18;table[2][1]=-0.5;table[2][2]=-0.5;table[2][3]=-0.18;table[2][4]=0.46;table[2][5]=1.42;table[3][0]=0.22;table[3][1]=-0.34;table[3][2]=-0.58;table[3][3]=-0.5;table[3][4]=-0.1;table[3][5]=0.62;table[4][0]=0.78;table[4][1]=-0.02;table[4][2]=-0.5;table[4][3]=-0.66;table[4][4]=-0.5;table[4][5]=-0.02;table[5][0]=1.5;table[5][1]=0.46;table[5][2]=-0.26;table[5][3]=-0.66;table[5][4]=-0.74;table[5][5]=-0.5;for(ii=0;ii<length_x;ii++)xx[ii]=0.08*ii;for(ii=0;ii<length_y;ii++)yy[ii]=0.5+0.05*ii;for(ii=0;ii<length_x;ii++) //生成U阵{for(jj=0;jj<length_y;jj++){ite_equations(&tt,&uu,xx[ii],yy[jj]); //牛顿迭代解方程，输入一点（x,y）输出一点(t,u)U[ii][jj]=(double)interpolation(tt,uu,(double **)table); //输入一点，输出一点插值}}fprintf(fd,"\nf(x,y):\n");for(ii=0;ii<x_length;ii++){for(jj=0;jj<=y_length;jj++){fprintf(fd,"%.11e ",U[ii][jj]);}fprintf(fd,"\n");}for(k_max=1;k_max<=6;k_max++){for(kk=0;kk<(k_max+1);kk++) //生成B阵{for(ii=0;ii<length_x;ii++){B[ii][kk]=pow(xx[ii],kk);}}for(kk=0;kk<(k_max+1);kk++){for(ii=0;ii<length_y;ii++) //生成G阵{G[ii][kk]=pow(yy[ii],kk);}}for(ii=0;ii<=k_max;ii++) //BT * B{for(jj=0;jj<=k_max;jj++){W[ii][jj]=row_multi_sum((double **)B,ii,jj,length_x);}}for(ii=0;ii<=k_max;ii++) //GT * G{for(jj=0;jj<=k_max;jj++){M[ii][jj]=row_multi_sum((double **)G,ii,jj,length_y);}}for(ii=0;ii<=k_max;ii++) //F=B'*Ufor(jj=0;jj<length_y;jj++){F[ii][jj]=0;for(kk=0;kk<length_x;kk++)F[ii][jj]=F[ii][jj]+B[kk][ii]*U[kk][jj];}for(ii=0;ii<=k_max;ii++) //E=F*Gfor(jj=0;jj<=k_max;jj++){E[ii][jj]=0;for(kk=0;kk<length_y;kk++)E[ii][jj]=E[ii][jj]+F[ii][kk]*G[kk][jj];}Matrix_LU((double **)W,(double **)E,(double **)M,(double**)C,k_max,0);//(double **)XX,for(ii=0;ii<length_x;ii++)for(jj=0;jj<length_y;jj++){sum=0;for(r=0;r<=k_max;r++)for(s=0;s<=k_max;s++)sum=sum+C[r][s]*pow(xx[ii],r)*pow(yy[jj],s);pp[ii][jj]=sum;}sum=0;for(ii=0;ii<length_x;ii++)for(jj=0;jj<length_y;jj++){sum=sum+(pp[ii][jj]-U[ii][jj])*(pp[ii][jj]-U[ii][jj]);}fprintf(fd,"\nk=%d, sigma=%.11e\n",k_max,sum);if(sum<=1e-7){fprintf(fd,"\nC_rs:\n");for(r=0;r<=k_max;r++){for(s=0;s<=k_max;s++){fprintf(fd,"%.11e ",C[r][s]);}fprintf(fd,"\n");}break;}}for(ii=0;ii<x_length;ii++)xx[ii]=0.1*(ii+1);for(ii=0;ii<y_length;ii++)yy[ii]=0.5+0.2*(ii+1);for(ii=0;ii<x_length;ii++) //生成f*阵{for(jj=0;jj<y_length;jj++){ite_equations(&tt,&uu,xx[ii],yy[jj]); //牛顿迭代解方程，输入一点（x,y）输出一点(t,u)U[ii][jj]=(double)interpolation(tt,uu,(double **)table); //输入一点，输出一点插值}}fprintf(fd,"\nf(x*,y*):\n");for(ii=0;ii<x_length;ii++){for(jj=0;jj<y_length;jj++){fprintf(fd,"%-10f ",U[ii][jj]);}fprintf(fd,"\n");}for(ii=0;ii<x_length;ii++) //生成p*阵{for(jj=0;jj<y_length;jj++){sum=0;for(r=0;r<=k_max;r++)for(s=0;s<=k_max;s++)sum=sum+C[r][s]*pow(xx[ii],r)*pow(yy[jj],s);pp[ii][jj]=sum;}}fprintf(fd,"\np(x*,y*):\n");for(ii=0;ii<x_length;ii++){for(jj=0;jj<y_length;jj++){fprintf(fd,"%-10f ",pp[ii][jj]);}fprintf(fd,"\n");}fclose(fd);}2)矩阵2列相乘#include<stdio.h>#include<math.h>#define size 10#define TTT(a,b) *((double *)TTT+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组double row_multi_sum(double **B,int i,int j,int length) //矩阵列相乘{double sum=0;double **TTT;int kk=0;TTT=B;for(kk=0;kk<length;kk++){sum=sum+(TTT(kk,i))*(TTT(kk,j));}return sum;}3)9点分片双二次插值#include<stdio.h>#include<math.h>#define size 6#define table(a,b) *((double *)table+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组double interpolation(double tt,double uu,double **table) //输入一点，输出一点插值{int ii=0;int jj=0;int section_t=0;int section_u=0;double answer=0;double temp=0;double temp_u[3]={0};double temp_t[3]={0};double my_t[6]={0};double my_u[6]={0};for(ii=1;ii<6;ii++){my_u[ii]=my_u[ii-1]+0.4;my_t[ii]=my_t[ii-1]+0.2;}//判断输入点在哪个区域//判断u用哪段插值if ((uu>=0)&&(uu<0.6))section_u=1;else if((uu>=0.6)&&(uu<1.0))section_u=2;else if((uu>=1.0)&&(uu<1.4))section_u=3;else if((uu>=1.4)&&(uu<=2))section_u=4;else;//判断t用哪段插值if ((tt>=0)&&(tt<0.3))section_t=1;else if((tt>=0.3)&&(tt<0.5))section_t=2;else if((tt>=0.5)&&(tt<0.7))section_t=3;else if((tt>=0.7)&&(tt<=1.0))section_t=4;else;temp_u[0]=(uu-my_u[section_u])*(uu-my_u[section_u+1])/((my_u[section _u-1]-my_u[section_u])*(my_u[section_u-1]-my_u[section_u+1]));temp_u[1]=(uu-my_u[section_u-1])*(uu-my_u[section_u+1])/((my_u[secti on_u]-my_u[section_u-1])*(my_u[section_u]-my_u[section_u+1]));temp_u[2]=(uu-my_u[section_u-1])*(uu-my_u[section_u])/((my_u[section _u+1]-my_u[section_u-1])*(my_u[section_u+1]-my_u[section_u]));temp_t[0]=(tt-my_t[section_t])*(tt-my_t[section_t+1])/((my_t[section _t-1]-my_t[section_t])*(my_t[section_t-1]-my_t[section_t+1]));temp_t[1]=(tt-my_t[section_t-1])*(tt-my_t[section_t+1])/((my_t[secti on_t]-my_t[section_t-1])*(my_t[section_t]-my_t[section_t+1]));temp_t[2]=(tt-my_t[section_t-1])*(tt-my_t[section_t])/((my_t[section _t+1]-my_t[section_t-1])*(my_t[section_t+1]-my_t[section_t]));answer=0;for(ii=0;ii<3;ii++){answer=answer+temp_t[ii]*(temp_u[0]*table((section_t-1+ii),section_u -1)+temp_u[1]*table((section_t-1+ii),section_u)+temp_u[2]*table((sectio n_t-1+ii),section_u+1));}return answer;}4)Newton迭代求解非线性方程组#include<stdio.h>#include<math.h>void ite_equations(double *tt,double *uu,double xx,double yy) //牛顿迭代解非线性方程组{double t=1;double u=1;double v=1;double w=1;double d[4]={0}; //增量double x=xx;double y=yy;double F_Jacob[4][4]={0};double F[4]={0};double L[9][9]={0};double U[9][9]={0};double temp1=0;double temp2=0;int ii=0;int jj=0;int kk=0;for(ii=0;ii<=1e2;ii++){F_Jacob[0][0]=-0.5*sin(t);F_Jacob[0][1]=1;F_Jacob[0][2]=1;F_Jacob[0][3]=1;F_Jacob[1][0]=1;F_Jacob[1][1]=0.5*cos(u);F_Jacob[1][2]=1;F_Jacob[1][3]=1;F_Jacob[2][0]=0.5;F_Jacob[2][1]=1;F_Jacob[2][2]=-sin(v);F_Jacob[2][3]=1;F_Jacob[3][0]=1;F_Jacob[3][1]=0.5;F_Jacob[3][2]=1;F_Jacob[3][3]=cos(w);F[0]=0.5*cos(t)+u+v+w-x-2.67;F[1]=t+0.5*sin(u)+v+w-y-1.07;F[2]=0.5*t+u+cos(v)+w-x-3.74;F[3]=t+0.5*u+v+sin(w)-y-0.79;LU((double **)L,(double **)U,(double **)F_Jacob,(double **)F);huidai((double **)L,(double **)U,(double **)F,(double **)d);t=t+d[0];u=u+d[1];v=v+d[2];w=w+d[3];temp1=sqrt(t*t+u*u+v*v+w*w);temp2=sqrt(d[0]*d[0]+d[1]*d[1]+d[2]*d[2]+d[3]*d[3]);if((temp2/temp1)<1e-12){*tt=t;*uu=u;break;}}5)解线性方程组用到LU分解#include<stdio.h>#include<math.h>#define size 4//#define halfBW 3#define L(a,b) *((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组#define U(a,b) *((double *)U+(a)*(size)+b)#define A(a,b) *((double *)A+(a)*(size)+b) //A 的宏定义，以方便操作二维数组//#define B(a,b) *((double *)B+(a)*(size)+b) //B 的宏定义，以方便操作二维数组//#define X(a,b) *((double *)X+(a)*(size)+b) //A 的宏定义，以方便操作二维数组#define b(ii) *((double *)b+ii)void LU(double **LL,double **UU,double **AA,double **bb) // //int a //需要选主元，{double **L;double **U;double **A;double A_temp[size][size]={0};double **b;double y[size]={0};double x[size]={0};int kk,ii,jj,pp,qq;double s,s_temp=0;double sum_answer=0;double temp=0;int s_count=0;int n=size;////////////////L=LL;U=UU;A=AA;b=bb;/////////////kk=0;ii=0;jj=0;pp=0;qq=0;//////////////////////////////////for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)A_temp[ii][jj]=A(ii,jj);for(ii=0;ii<n;ii++)L(ii,ii)=1.0;for (kk=0;kk<n;kk++) //总循环求出L U{ sum_answer=0;s_count=kk;for(pp=kk;pp<n;pp++)//doolitte 选主元{if (kk==0)sum_answer=0;elsefor(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(pp,qq)*U(qq,pp);s=A_temp[pp][kk]-sum_answer;if (s>s_temp){s_temp=s;s_count=pp;}else if (s<(-1.0)*s_temp){s_temp=(-1.0)*s;s_count=pp;}else;}for(pp=0;pp<n;pp++) //换位A_temp{temp=A_temp[s_count][pp];A_temp[s_count][pp]=A_temp[kk][pp];A_temp[kk][pp]=temp;}for(pp=0;pp<kk;pp++) //换位L{temp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;}temp=b(s_count);b(s_count)=b(kk);b(kk)=temp;for (ii=kk;ii<n;ii++) //先求一行U{if(kk==0)U(kk,ii)=A_temp[kk][ii];else{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_temp[kk][ii]-sum_answer;}}for(ii=kk+1;ii<n;ii++) //求出一列L{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_temp[ii][kk]-sum_answer)/U(kk,kk);}}}6)回代求解方程组#include<stdio.h>#include<math.h>#define size 4#define L(a,b) *((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组#define U(a,b) *((double *)U+(a)*(size)+b)#define x(a) *((double *)x+a) //L U 的宏定义，以方便操作二维数组#define y(a) *((double *)y+a)#define b(a) *((double *)b+a)//#define xx(a) xx[a]//A*uu=yy => L*U*uu=yy =>L*xx=yy && U*uu=xxvoid huidai(double **LL,double **UU,double **bb,double **xx) //回代法求xx, LUx=b{double **L;double **U;double **b;double **x;double y[size]={0};double sum_answer=0;int ii=0;int jj=0;L=LL;U=UU;b=bb;x=xx;//迭代Ly=by[0]=-b(0);for(ii=1;ii<size;ii++){sum_answer=0;for(jj=1;jj<ii;jj++){sum_answer=sum_answer+L(ii,jj)*y[jj];}y[ii]=-b(ii)-sum_answer;}//迭代Ux=yx(size-1)=y[size-1]/U((size-1),(size-1));for(ii=size-2;ii>=0;ii--){sum_answer=0;for(jj=ii+1;jj<size;jj++){sum_answer=sum_answer+U(ii,jj)*x(jj);}x(ii)=(y[ii]-sum_answer)/U(ii,ii);}}7)求解C=（B’*B）-1*(B’*U*G)*(G’*G)#include<stdio.h>#include<math.h>#define size 10//#define halfBW 3#define L(a,b) L[a][b]//*((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组#define U(a,b) U[a][b]//*((double *)U+(a)*(size)+b)#define A(a,b) *((double *)A+(a)*(size)+b) //A 的宏定义，以方便操作二维数组#define B(a,b) *((double *)B+(a)*(size)+b) //B 的宏定义，以方便操作二维数组#define M(a,b) *((double *)M+(a)*(size)+b) //B 的宏定义，以方便操作二维数组#define C(a,b) *((double *)C+(a)*(size)+b) //B 的宏定义，以方便操作二维数组//#define X(a,b) *((double *)X+(a)*(size)+b) //A 的宏定义，以方便操作二维数组#define X(a,b) X[a][b]//#define b(ii) *((double *)b+ii)void Matrix_LU(double **AA,double **BB,double **MM,double **CC,intk_max,int mode){double **A;double **B;double **M;double **C;double X[size][size]={0};int n=k_max+1;double A_temp[size][size]={0}; double L[size][size]={0};double U[size][size]={0};double Y[size][size]={0};double **b;int kk,ii,jj,pp,qq;double s,s_temp=0;double sum_answer=0;double temp=0;int s_count=0;//////////////////L=LL;//U=UU;A=AA;B=BB;M=MM;C=CC;//X=XX;/////////////kk=0;ii=0;jj=0;pp=0;qq=0;////////////////////////////////// for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)A_temp[ii][jj]=A(ii,jj);for(ii=0;ii<n;ii++)L(ii,ii)=1.0;for (kk=0;kk<n;kk++) //总循环求出L U{ sum_answer=0;s_count=kk;for(pp=kk;pp<n;pp++)//doolitte 选主元{if (kk==0)sum_answer=0;elsefor(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(pp,qq)*U(qq,pp);s=A_temp[pp][kk]-sum_answer;if (s>s_temp){s_temp=s;s_count=pp;}else if (s<(-1.0)*s_temp){s_temp=(-1.0)*s;s_count=pp;}else;}for(pp=0;pp<n;pp++) //换位A_temp{temp=A_temp[s_count][pp];A_temp[s_count][pp]=A_temp[kk][pp];A_temp[kk][pp]=temp;}for(pp=0;pp<kk;pp++) //换位L{temp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;}for(pp=0;pp<n;pp++) //换位B /////////////{temp=B(s_count,pp);B(s_count,pp)=B(kk,pp);B(kk,pp)=temp;}for (ii=kk;ii<n;ii++) //先求一行U{if(kk==0)U(kk,ii)=A_temp[kk][ii];else{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_temp[kk][ii]-sum_answer;}}for(ii=kk+1;ii<n;ii++) //求出一列L{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_temp[ii][kk]-sum_answer)/U(kk,kk);}}//迭代Ly=bfor(kk=0;kk<n;kk++){Y[0][kk]=B(0,kk);for(ii=1;ii<n;ii++){sum_answer=0;for(jj=0;jj<ii;jj++){sum_answer=sum_answer+L(ii,jj)*Y[jj][kk];}Y[ii][kk]=B(ii,kk)-sum_answer;}//迭代Ux=yX((n-1),kk)=Y[n-1][kk]/U((n-1),(n-1));for(ii=n-2;ii>=0;ii--){sum_answer=0;for(jj=ii+1;jj<n;jj++){sum_answer=sum_answer+U(ii,jj)*X(jj,kk);}X(ii,kk)=(Y[ii][kk]-sum_answer)/U(ii,ii);}}////////////////////////////////////////////////////////////// // 再操作，（G'G）C'=X'//////////////////////////////////for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)A_temp[ii][jj]=M(ii,jj);for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)B(ii,jj)=X(jj,ii);for(ii=0;ii<n;ii++)L(ii,ii)=1.0;for (kk=0;kk<n;kk++) //总循环求出L U{ sum_answer=0;s_count=kk;for(pp=kk;pp<n;pp++)//doolitte 选主元{if (kk==0)sum_answer=0;elsefor(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(pp,qq)*U(qq,pp);s=A_temp[pp][kk]-sum_answer;if (s>s_temp){s_temp=s;s_count=pp;}else if (s<(-1.0)*s_temp){s_temp=(-1.0)*s;s_count=pp;}else;}for(pp=0;pp<n;pp++) //换位A_temp{temp=A_temp[s_count][pp];A_temp[s_count][pp]=A_temp[kk][pp];A_temp[kk][pp]=temp;}for(pp=0;pp<kk;pp++) //换位L{temp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;}for(pp=0;pp<n;pp++) //换位B /////////////{temp=B(s_count,pp);B(s_count,pp)=B(kk,pp);B(kk,pp)=temp;}for (ii=kk;ii<n;ii++) //先求一行U{if(kk==0)U(kk,ii)=A_temp[kk][ii];else{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_temp[kk][ii]-sum_answer;}}for(ii=kk+1;ii<n;ii++) //求出一列L{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_temp[ii][kk]-sum_answer)/U(kk,kk);}}//迭代Ly=bfor(kk=0;kk<n;kk++){Y[0][kk]=B(0,kk);for(ii=1;ii<n;ii++){sum_answer=0;for(jj=0;jj<ii;jj++){sum_answer=sum_answer+L(ii,jj)*Y[jj][kk];}Y[ii][kk]=B(ii,kk)-sum_answer;}//迭代Ux=yX((n-1),kk)=Y[n-1][kk]/U((n-1),(n-1));for(ii=n-2;ii>=0;ii--){sum_answer=0;for(jj=ii+1;jj<n;jj++){sum_answer=sum_answer+U(ii,jj)*X(jj,kk);}X(ii,kk)=(Y[ii][kk]-sum_answer)/U(ii,ii);}}for(ii=0;ii<n;ii++) // 给C赋值，直接转置for(jj=0;jj<n;jj++)C(ii,jj)=X(jj,ii);}。