北航数值分析大作业三.doc

数值分析大作业三一、 算法设计方案：1、使用牛顿迭代法(简单迭代法不收敛)，对原题中给出的，()的11*21组分别求出原题中方程组的一组解，于是得到一组和对应的。2、对于已求出的，使用分片二次代数插值法对原题中关于的数表进行插值得到。于是产生了z=f(x,y)的11*21个数值解。3、从k=1开始逐渐增大k的值，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，得到每次的。当时结束计算，输出拟合结果。4、计算的值并输出结果，以观察逼近的效果。其中。二、源程序：#include #include #define N 6 /矩阵的阶；#define M 4 /方程组未知数个数； #define EPSL 1.0e-12 /迭代的精度水平；#define K 10 /最小二乘法拟合时的次数上限；#define SIGSUM 1.0e-7#define L 500 /迭代最大次数；#define X 11#define Y 21double mat_uN=0,0.4,0.8,1.2,1.6,2.0, mat_tN=0,0.2,0.4,0.6,0.8,1.0, mat_ztuNN= -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, , mat_val_zXY=0, init_tuvw4=1,2,1,2, mat_crsKK=0;FILE *p;int main()int i,j,x,y,kmin,tmp=0;double tmpval;int getval_z(double,double,int,int);int leasquare();void result_out(int); p=fopen(Result.txt,w+); fprintf(p,1、差值得到的数表为：n); fclose(p); for(i=0;iX;i+) for(j=0;jY;j+) tmp+=getval_z(0.08*i,0.5+0.05*j,i,j); if(!tmp) printf(迭代求解 z=f(x,y) 完毕。n); else printf(迭代超过最大次数，计算结果可能不准确！n); if(kmin=leasquare() printf(近似表达式已求出！n); for(i=0;i8;i+) for(j=0;j5;j+) tmp+=getval_z(0.1*(i+1),0.5+0.2*(j+1),i,j); if(!tmp) printf(迭代求解 z=f(x*,y*) 完毕。n); for(x=0;x8;x+) for(y=0;y5;y+) tmpval=0; for(i=0;ikmin;i+) for(j=0;jkmin;j+) tmpval=tmpval+mat_crsij*pow(0.1*(x+1),i)*pow(0.5+0.2*(y+1),j); p=fopen(Result.txt, at+); fprintf(p,x*%d=%f, y*%d=%f: f(x*%d,y*%d)=%.12le, p(x*%d,y*%d)=%.12len,x+1,0.1*(x+1),y+1,0.5+0.2*(y+1),x+1,y+1,mat_val_zxy,x+1,y+1,tmpval); fclose(p); printf(结果见Result.txt！); getchar(); else printf(近似表达式未求出，增大K值后重试n); getchar(); /牛顿迭代法求方程的解int getval_z(double x,double y,int i,int j) double vec_tuvwM,vec_deltaM,fdaoMM,fM,mo_k,mo_delta_k,shang,t,u; int n=0,k,flag_max,flag_stat; void solve(double fdaoMM,double vec_deltaM,double fM); double interpolation(double,double); flag_stat = 1; for(k=0;kM;k+) vec_tuvwk=init_tuvwk; do f0=-1.0*(0.5*cos(vec_tuvw0)+vec_tuvw1+vec_tuvw2+vec_tuvw3-x-2.67); f1=-1.0*(vec_tuvw0+0.5*sin(vec_tuvw1)+vec_tuvw2+vec_tuvw3-y-1.07); f2=-1.0*(0.5*vec_tuvw0+vec_tuvw1+cos(vec_tuvw2)+vec_tuvw3-x-3.74); f3 =-1.0*(vec_tuvw0+0.5*vec_tuvw1+vec_tuvw2+sin(vec_tuvw3)-y-0.79); fdao00=-0.5*sin(vec_tuvw0); fdao01=1; fdao02=1; fdao03=1; fdao10=1; fdao11=0.5*cos(vec_tuvw1); fdao12=1; fdao13=1; fdao20=0.5; fdao21=1; fdao22=-1*sin(vec_tuvw2); fdao23=1; fdao30=1; fdao31=0.5; fdao32=1; fdao33=cos(vec_tuvw3); solve(fdao, vec_delta, f); flag_max = 0; for(k=1;kfabs(vec_deltaflag_max) flag_max=k; mo_delta_k=vec_deltaflag_max; flag_max=0; for(k=1;kfabs(vec_tuvwflag_max) flag_max=k; mo_k=vec_tuvwflag_max; shang=fabs(mo_delta_k/mo_k); if(shangEPSL) t=vec_tuvw0+vec_delta0,u=vec_tuvw1+vec_delta1; flag_stat=0; break; else for(k=0;kM;k+) vec_tuvwk+=vec_deltak; while(n+ = L); if(!flag_stat) p=fopen(Result.txt, at+); fprintf(p,x%d=%f, y%d=%f: f(x%d,y%d)=%.12len,i,0.08*i,j,0.5+0.05*j,i,j,mat_val_zij=interpolation(u,t); fclose(p); return flag_stat;/消元法求方程的解void solve(double fdaoMM,double vec_deltaM,double fM) int i,j,k,ik;double tmp,mik;for(k=0;kM-1;k+) ik = k; for(i=k;iM;i+) if(fabs(fdaoikk)fabs(fdaoik) ik=i; for(j=k;jM;j+) tmp=fdaokj;fdaokj=fdaoikj;fdaoikj=tmp; tmp=fk;fk=fik;fik=tmp; for(i=k+1;iM;i+) mik=fdaoik/fdaokk; for(j=k;j=0;k-) tmp=0; for(j=k+1;jM;j+) tmp+=fdaokj*vec_deltaj; vec_deltak=(fk-tmp)/fdaokk; /二次插值double interpolation(double u, double t) int r,i,j,k; double coe_u3,coe_t3,z; z=0; r=int(fabs(t/0.2)+0.5); k=int(fabs(u/0.4)+0.5); if(r=0) r=1; if(r=5) r=4; if(k=0) k=1; if(k=5) k=4; coe_u0=(u-mat_uk)*(u-mat_uk+1)/(mat_uk-1-mat_uk)/(mat_uk-1-mat_uk+1); coe_u1=(u-mat_uk-1)*(u-mat_uk+1)/(mat_uk-mat_uk-1)/(mat_uk-mat_uk+1); coe_u2=(u-mat_uk-1)*(u-mat_uk)/(mat_uk+1-mat_uk-1)/(mat_uk+1-mat_uk); coe_t0=(t-mat_tr)*(t-mat_tr+1)/(mat_tr-1-mat_tr)/(mat_tr-1-mat_tr+1); coe_t1=(t-mat_tr-1)*(t-mat_tr+1)/(mat_tr-mat_tr-1)/(mat_tr-mat_tr+1); coe_t2=(t-mat_tr-1)*(t-mat_tr)/(mat_tr+1-mat_tr-1)/(mat_tr+1-mat_tr); for(i=r-1;i=r+1;i+) for(j=k-1;j=k+1;j+) z+=coe_ti-r+1*coe_uj-k+1*mat_ztuij; return z;/曲面拟合int leasquare() int x,y,i,j,k,k_max,k_now; double vec_xX,vec_yY,mat_btbKK=0,mat_gtgKK=0,mat_btbbtKX,mat_btbbtuKY,tmp,sigma; void inv(double matrixKK,int); k_now=1; for(i=0;iX;i+) vec_xi=0.08*i; for(j=0;jY;j+) vec_yj=0.5+0.05*j; do for(i=0;ik_now;i+) for(j=0;jk_now;j+) tmp=0; for(k=0;kX;k+) tmp+=pow(vec_xk,i)*pow(vec_xk,j); mat_btbij=tmp; inv(mat_btb,k_now); for(i=0;ik_now;i+) for(j=0;jk_now;j+) tmp=0; for(k=0;kY;k+) tmp+=pow(vec_yk,i)*pow(vec_yk,j); mat_gtgij=tmp; inv(mat_gtg, k_now); for(i=0;ik_now;i+) for(j=0;jX;j+) tmp=0; for(k=0;kk_now;k+) tmp+=mat_btbik*pow(vec_xj,k); mat_btbbtij=tmp; for(i=0;ik_now;i+) for(j=0;jY;j+) tmp=0; for(k=0;kX;k+) tmp+=mat_btbbtik*mat_val_zkj; mat_btbbtuij=tmp; for(i=0;ik_now;i+) for(j=0;jk_now;j+) tmp=0; for(k=0;kY;k+) tmp+=mat_btbbtuik*pow(vec_yk,j); mat_btbij=tmp; for(i=0;ik_now;i+) for(j=0;jk_now;j+) tmp=0; for(k=0;kk_now;k+) tmp+=mat_btbik*mat_gtgkj; mat_crsij=tmp; sigma=0; for(x=0;xX;x+) for(y=0;yY;y+) tmp=0; for(i=0;ik_now;i+) for(j=0;jk_now;j+) tmp+=mat_crsij*pow(0.08*x,i)*pow(0.5+0.05*y,j); sigma+=(tmp-mat_val_zxy)*(tmp-mat_val_zxy); p=fopen(Result.txt, at+); fprintf(p,k=%d =%.12lenn,k_now-1,sigma); fclose(p); if(sigmaSIGSUM) p=fopen(Result.txt, at+); fprintf(p,Crs%d%d=nn,k_now,k_now); for(i=0;ik_now;i+) fprintf(p,); for(j=0;jk_now-1;j+) fprintf(p,%.12le, , mat_crsij); fprintf(p,%.12le, mat_crsik_now-1); fprintf(p,n); fprintf(p,n); fclose(p); return k_now; while(k_now+K);return 0; /dolittle分解void inv(double matrixKK, int rank) int i,j,k,t,n; double mat_tmpKK=0,matrix_bakKK,vec_tmpK,vec_xK=0,vec_dxK=0,vec_rK=0,tmp,max_x,max_dx; void preprocess(double KK,double K,int); void LUDecomposition(double KK,int); void LUSolve(double KK,double K,double K,int); n=0; for(i=0;irank;i+) for(j=0;jrank;j+) matrix_bakij=mat_tmpij=matrixij; vec_tmpi=0; LUDecomposition(mat_tmp,rank); for(i=0;irank;i+) if(i=0) vec_tmpi=1; else vec_tmpi-1=0; vec_tmpi=1; LUSolve(mat_tmp,vec_x,vec_tmp,rank); while(n+=3) for(j=0;jrank;j+) tmp=0; for(t=0;trank;t+) tmp+=matrix_bakjt*vec_xt; vec_rj=vec_tmpj-tmp; LUSolve(mat_tmp,vec_dx,vec_r,rank); max_x=vec_x0,max_dx=vec_dx0; for(j=0;jmax_x) max_x=vec_xj; if(vec_dxjmax_dx) max_dx=vec_dxj; for(j=0;jrank;j+) vec_xj+=vec_dxj; for(j=0;jrank;j+) matrixji=vec_xj; /改变条件数void preprocess(double matrixKK,double vec_tmpK,int rank) int i,k; double tmp; for(i=0;irank;i+) tmp=matrixi0; for(k=0;ktmp) tmp=matrixik; for(k=0;krank;k+) matrixik/=tmp; vec_tmpi/=tmp; /LU分解void LUDecomposition(double matrixKK,int rank) int i,j,k,t; double tmp; for(k=0;krank;k+) for(i=k;irank;i+) tmp=0; for(t=0;tk;t+) tmp+=matrixit*matrixtk; matrixik-=tmp; if(krank-1) for(j=k+1;jrank;j+) tmp=0; for(t=0;tk;t+) tmp+=matrixkt*matrixtj; matrixkj=(matrixkj-tmp)/matrixkk; else /求解矩阵的逆void LUSolve(double matrixKK,double vec_xK,double vec_tmpK,int rank) int i,j,t; double tmp; vec_x0=vec_tmp0/matrix00; for(i=1;irank;i+) tmp=0; for(t=0;t=0;i-) tmp=0; for(t=i+1;trank;t+) tmp+=matrixit*vec_xt; vec_xi-=tmp; 三、 运算结果：1、差值得到的数表为：x0=0.000000, y0=0.500000: f(x0,y0)=4.465040184799e-001x0=0.000000, y1=0.550000: f(x0,y1)=3.246832629274e-001x0=0.000000, y2=0.600000: f(x0,y2)=2.101596866825e-001x0=0.000000, y3=0.650000: f(x0,y3)=1.030436083159e-001x0=0.000000, y4=0.700000: f(x0,y4)=3.401895562658e-003x0=0.000000, y5=0.750000: f(x0,y5)=-8.873581363801e-002x0=0.000000, y6=0.800000: f(x0,y6)=-1.733716327497e-001x0=0.000000, y7=0.850000: f(x0,y7)=-2.505346114666e-001x0=0.000000, y8=0.900000: f(x0,y8)=-3.202765063876e-001x0=0.000000, y9=0.950000: f(x0,y9)=-3.826680661097e-001x0=0.000000, y10=1.000000: f(x0,y10)=-4.377957667384e-001x0=0.000000, y11=1.050000: f(x0,y11)=-4.857589414438e-001x0=0.000000, y12=1.100000: f(x0,y12)=-5.266672548835e-001x0=0.000000, y13=1.150000: f(x0,y13)=-5.606384797965e-001x0=0.000000, y14=1.200000: f(x0,y14)=-5.877965387677e-001x0=0.000000, y15=1.250000: f(x0,y15)=-6.082697790899e-001x0=0.000000, y16=1.300000: f(x0,y16)=-6.221894528764e-001x0=0.000000, y17=1.350000: f(x0,y17)=-6.296883781856e-001x0=0.000000, y18=1.400000: f(x0,y18)=-6.308997600028e-001x0=0.000000, y19=1.450000: f(x0,y19)=-6.259561525454e-001x0=0.000000, y20=1.500000: f(x0,y20)=-6.149885466094e-001x1=0.080000, y0=0.500000: f(x1,y0)=6.380152265102e-001x1=0.080000, y1=0.550000: f(x1,y1)=5.066117551462e-001x1=0.080000, y2=0.600000: f(x1,y2)=3.821763692772e-001x1=0.080000, y3=0.650000: f(x1,y3)=2.648634911536e-001x1=0.080000, y4=0.700000: f(x1,y4)=1.547802002848e-001x1=0.080000, y5=0.750000: f(x1,y5)=5.199268349093e-002x1=0.080000, y6=0.800000: f(x1,y6)=-4.346804020491e-002x1=0.080000, y7=0.850000: f(x1,y7)=-1.316010567885e-001x1=0.080000, y8=0.900000: f(x1,y8)=-2.124310883088e-001x1=0.080000, y9=0.950000: f(x1,y9)=-2.860045510580e-001x1=0.080000, y10=1.000000: f(x1,y10)=-3.523860789794e-001x1=0.080000, y11=1.050000: f(x1,y11)=-4.116554565222e-001x1=0.080000, y12=1.100000: f(x1,y12)=-4.639049115188e-001x1=0.080000, y13=1.150000: f(x1,y13)=-5.092367247005e-001x1=0.080000, y14=1.200000: f(x1,y14)=-5.477611179623e-001x1=0.080000, y15=1.250000: f(x1,y15)=-5.795943883391e-001x1=0.080000, y16=1.300000: f(x1,y16)=-6.048572588895e-001x1=0.080000, y17=1.350000: f(x1,y17)=-6.236734213318e-001x1=0.080000, y18=1.400000: f(x1,y18)=-6.361682484133e-001x1=0.080000, y19=1.450000: f(x1,y19)=-6.424676566900e-001x1=0.080000, y20=1.500000: f(x1,y20)=-6.426971026996e-001x2=0.160000, y0=0.500000: f(x2,y0)=8.400813957651e-001x2=0.160000, y1=0.550000: f(x2,y1)=6.997641656726e-001x2=0.160000, y2=0.600000: f(x2,y2)=5.660614423514e-001x2=0.160000, y3=0.650000: f(x2,y3)=4.391716081175e-001x2=0.160000, y4=0.700000: f(x2,y4)=3.192421380407e-001x2=0.160000, y5=0.750000: f(x2,y5)=2.063761923874e-001x2=0.160000, y6=0.800000: f(x2,y6)=1.006385238914e-001x2=0.160000, y7=0.850000: f(x2,y7)=2.060740067835e-003x2=0.160000, y8=0.900000: f(x2,y8)=-8.935402476698e-002x2=0.160000, y9=0.950000: f(x2,y9)=-1.736269688648e-001x2=0.160000, y10=1.000000: f(x2,y10)=-2.507999561599e-001x2=0.160000, y11=1.050000: f(x2,y11)=-3.209322694446e-001x2=0.160000, y12=1.100000: f(x2,y12)=-3.840977350046e-001x2=0.160000, y13=1.150000: f(x2,y13)=-4.403821754175e-001x2=0.160000, y14=1.200000: f(x2,y14)=-4.898811523126e-001x2=0.160000, y15=1.250000: f(x2,y15)=-5.326979655338e-001x2=0.160000, y16=1.300000: f(x2,y16)=-5.689418792921e-001x2=0.160000, y17=1.350000: f(x2,y17)=-5.987265495151e-001x2=0.160000, y18=1.400000: f(x2,y18)=-6.221686297503e-001x2=0.160000, y19=1.450000: f(x2,y19)=-6.393865356972e-001x2=0.160000, y20=1.500000: f(x2,y20)=-6.504993507878e-001x3=0.240000, y0=0.500000: f(x3,y0)=1.051515091801e+000x3=0.240000, y1=0.550000: f(x3,y1)=9.029274308302e-001x3=0.240000, y2=0.600000: f(x3,y2)=7.605802668593e-001x3=0.240000, y3=0.650000: f(x3,y3)=6.247151981455e-001x3=0.240000, y4=0.700000: f(x3,y4)=4.955197560009e-001x3=0.240000, y5=0.750000: f(x3,y5)=3.731340427746e-001x3=0.240000, y6=0.800000: f(x3,y6)=2.576567488723e-001x3=0.240000, y7=0.850000: f(x3,y7)=1.491505594102e-001x3=0.240000, y8=0.900000: f(x3,y8)=4.764698677337e-002x3=0.240000, y9=0.950000: f(x3,y9)=-4.684932320146e-002x3=0.240000, y10=1.000000: f(x3,y10)=-1.343567603849e-001x3=0.240000, y11=1.050000: f(x3,y11)=-2.149133449274e-001x3=0.240000, y12=1.100000: f(x3,y12)=-2.885737006348e-001x3=0.240000, y13=1.150000: f(x3,y13)=-3.554063647857e-001x3=0.240000, y14=1.200000: f(x3,y14)=-4.154913964886e-001x3=0.240000, y15=1.250000: f(x3,y15)=-4.689182499695e-001x3=0.240000, y16=1.300000: f(x3,y16)=-5.157838831247e-001x3=0.240000, y17=1.350000: f(x3,y17)=-5.561910752001e-001x3=0.240000, y18=1.400000: f(x3,y18)=-5.902469305629e-001x3=0.240000, y19=1.450000: f(x3,y19)=-6.180615482412e-001x3=0.240000, y20=1.500000: f(x3,y20)=-6.397468392579e-001x4=0.320000, y0=0.500000: f(x4,y0)=1.271246751481e+000x4=0.320000, y1=0.550000: f(x4,y1)=1.115002018146e+000x4=0.320000, y2=0.600000: f(x4,y2)=9.646077272154e-001x4=0.320000, y3=0.650000: f(x4,y3)=8.203473694749e-001x4=0.320000, y4=0.700000: f(x4,y4)=6.824476781794e-001x4=0.320000, y5=0.750000: f(x4,y5)=5.510852085975e-001x4=0.320000, y6=0.800000: f(x4,y6)=4.263923859018e-001x4=0.320000, y7=0.850000: f(x4,y7)=3.084629956332e-001x4=0.320000, y8=0.900000: f(x4,y8)=1.973571296919e-001x4=0.320000, y9=0.950000: f(x4,y9)=9.310562085941e-002x4=0.320000, y10=1.000000: f(x4,y10)=-4.285992234034e-003x4=0.320000, y11=1.050000: f(x4,y11)=-9.483392529689e-002x4=0.320000, y12=1.100000: f(x4,y12)=-1.785729903640e-001x4=0.320000, y13=1.150000: f(x4,y13)=-2.555537790546e-001x4=0.320000, y14=1.200000: f(x4,y14)=-3.258401501575e-001x4=0.320000, y15=1.250000: f(x4,y15)=-3.895069883634e-001x4=0.320000, y16=1.300000: f(x4,y16)=-4.466382045995e-001x4=0.320000, y17=1.350000: f(x4,y17)=-4.973249517677e-001x4=0.320000, y18=1.400000: f(x4,y18)=-5.416640326994e-001x4=0.320000, y19=1.450000: f(x4,y19)=-5.797564797951e-001x4=0.320000, y20=1.500000: f(x4,y20)=-6.117062881476e-001x5=0.400000, y0=0.500000: f(x5,y0)=1.498321052481e+000x5=0.400000, y1=0.550000: f(x5,y1)=1.334998632066e+000x5=0.400000, y2=0.600000: f(x5,y2)=1.177125123739e+000x5=0.400000, y3=0.650000: f(x5,y3)=1.025024055020e+000x5=0.400000, y4=0.700000: f(x5,y4)=8.789600231743e-001x5=0.400000, y5=0.750000: f(x5,y5)=7.391451087035e-001x5=0.400000, y6=0.800000: f(x5,y6)=6.057448714871e-001x5=0.400000, y7=0.850000: f(x5,y7)=4.788838610666e-001x5=0.400000, y8=0.900000: f(x5,y8)=3.586506258818e-001x5=0.400000, y9=0.950000: f(x5,y9)=2.451022361964e-001x5=0.400000, y10=1.000000: f(x5,y10)=1.382683509285e-001x5=0.400000, y11=1.050000: f(x5,y11)=3.815486540699e-002x5=0.400000, y12=1.100