北航研究生数值分析上机作业四(报告+所有程序大全).doc

BY1102130 刘晚春 电子信息工程学院 数值分析上机作业4数值方法求解常微分方程及数值积分1. 算法设计i. 由，插值求解ii. 将，带入Newton迭代法求解非线性方程组iii. 由得到iv. 计算同理以上过程再进行3次得到，由经典R-K方法公式：可得。再由复化Simpson积分公式：可得积分值。2. 打印输出结果5位有效数字格式：t(i) y(i) x(i) z(i) p(i) q(i)0.0 2.0000e+0009.1445e-0011.6333e+0002.8181e+0008.6000e+0000.2 1.9533e+0008.9816e-0011.6605e+0002.8477e+0008.5396e+0000.4 1.9108e+0008.8411e-0011.6844e+0002.8765e+0008.4951e+0000.6 1.8722e+0008.7209e-0011.7048e+0002.9043e+0008.4660e+0000.8 1.8373e+0008.6194e-0011.7216e+0002.9307e+0008.4516e+0001.0 1.8056e+0008.5353e-0011.7347e+0002.9558e+0008.4516e+0001.2 1.7770e+0008.4675e-0011.7440e+0002.9792e+0008.4656e+0001.4 1.7513e+0008.4151e-0011.7495e+0003.0008e+0008.4932e+0001.6 1.7283e+0008.3772e-0011.7511e+0003.0205e+0008.5342e+0001.8 1.7078e+0008.3533e-0011.7489e+0003.0382e+0008.5883e+0002.0 1.6897e+0008.3432e-0011.7430e+0003.0537e+0008.6552e+0002.2 1.6739e+0008.3465e-0011.7334e+0003.0669e+0008.7349e+0002.4 1.6603e+0008.3632e-0011.7203e+0003.0779e+0008.8272e+0002.6 1.6489e+0008.3934e-0011.7039e+0003.0866e+0008.9320e+0002.8 1.6395e+0008.4374e-0011.6844e+0003.0929e+0009.0491e+0003.0 1.6320e+0008.6317e-0011.6204e+0003.0924e+0009.3945e+0003.2 1.6263e+0008.6990e-0011.5984e+0003.0947e+0009.5260e+0003.4 1.6226e+0008.7709e-0011.5769e+0003.0952e+0009.6532e+0003.6 1.6210e+0008.8477e-0011.5559e+0003.0940e+0009.7760e+0003.8 1.6216e+0008.9294e-0011.5354e+0003.0910e+0009.8945e+0004.0 1.6244e+0009.0160e-0011.5153e+0003.0865e+0001.0008e+0014.2 1.6294e+0009.1072e-0011.4958e+0003.0804e+0001.0117e+0014.4 1.6366e+0009.2027e-0011.4770e+0003.0728e+0001.0219e+0014.6 1.6462e+0009.3020e-0011.4590e+0003.0639e+0001.0315e+0014.8 1.6581e+0009.4045e-0011.4418e+0003.0537e+0001.0403e+0015.0 1.6724e+0009.0334e-0011.5304e+0003.0511e+0009.8131e+0005.2 1.6891e+0009.1466e-0011.5099e+0003.0369e+0009.9016e+0005.4 1.7082e+0009.2998e-0011.4823e+0003.0211e+0001.0034e+0015.6 1.7299e+0009.4964e-0011.4481e+0003.0042e+0001.0212e+0015.8 1.7544e+0009.7407e-0011.4082e+0002.9865e+0001.0438e+0016.0 1.7817e+0001.0038e+0001.3632e+0002.9686e+0001.0714e+0016.2 1.8122e+0001.0395e+0001.3140e+0002.9510e+0001.1043e+0016.4 1.8462e+0001.0817e+0001.2615e+0002.9342e+0001.1428e+0016.6 1.8840e+0001.1312e+0001.2068e+0002.9186e+0001.1871e+0016.8 1.9261e+0001.1889e+0001.1507e+0002.9049e+0001.2375e+0017.0 1.9732e+0001.3243e+0001.0353e+0002.9047e+0001.3693e+0017.2 2.0274e+0001.4003e+0009.8689e-0012.8963e+0001.4300e+0017.4 2.0875e+0001.4798e+0009.4271e-0012.8890e+0001.4907e+0017.6 2.1536e+0001.5626e+0009.0237e-0012.8826e+0001.5514e+0017.8 2.2262e+0001.6485e+0008.6549e-0012.8770e+0001.6118e+0018.0 2.3055e+0001.7373e+0008.3174e-0012.8720e+0001.6717e+0018.2 2.3917e+0001.8286e+0008.0082e-0012.8674e+0001.7309e+0018.4 2.4850e+0001.9224e+0007.7242e-0012.8631e+0001.7892e+0018.6 2.5853e+0002.0103e+0007.4900e-0012.8571e+0001.8391e+0018.8 2.6930e+0002.1088e+0007.2455e-0012.8532e+0001.8959e+0019.0 2.8084e+0002.2125e+0007.0107e-0012.8501e+0001.9546e+0019.2 2.9316e+0002.3222e+0006.7835e-0012.8480e+0002.0157e+0019.4 3.0628e+0002.4387e+0006.5623e-0012.8470e+0002.0798e+0019.6 3.2024e+0002.5627e+0006.3465e-0012.8472e+0002.1474e+0019.8 3.3505e+0002.6945e+0006.1362e-0012.8486e+0002.2188e+00110.0 3.5072e+0000.0000e+0000.0000e+0000.0000e+0000.0000e+000I = 1.99597053717e+00112位有效数字格式：t(i) y(i) x(i) z(i) p(i) q(i)0.0 2.00000000000e+0009-0011.63331102473e+0002.81808518237e+0008.60000000000e+0000.2 1.95327865082e+0008.98164241433e-0011.66048303187e+0002.84771740736e+0008.53960839959e+0000.4 1.91081324161e+0008.84105976349e-0011.68436757832e+0002.87651131035e+0008.49513396605e+0000.6 1.87224799520e+0008.72085264796e-0011.70479197341e+0002.90425238594e+0008.46597925076e+0000.8 1.83727150662e+0008.61940369583e-0011.72161140444e+0002.93073678797e+0008.45162816553e+0001.0 1.80561130796e+0008.53534537640e-0011.73471260234e+0002.95577332733e+0008.45163452750e+0001.2 1.77702931229e+0008.46753884324e-0011.74401667275e+0002.97918535843e+0008.46561262996e+0001.4 1.75131800336e+0008.41505645592e-0011.74948103464e+0003.00081255077e+0008.49322947826e+0001.6 1.72829726288e+0008.37716725063e-0011.75110045266e+0003.02051255798e+0008.53419840577e+0001.8 1.70781174751e+0008.35332477826e-0011.74890718752e+0003.03816260319e+0008.58827384401e+0002.0 1.68972874326e+0008.34315682686e-0011.74297031891e+0003.05366099808e+0008.65524706583e+0002.2 1.67393643665e+0008.34645660590e-0011.73339431614e+0003.06692860151e+0008.73494275605e+0002.4 1.66034255185e+0008.36317500229e-0011.72031694322e+0003.07791020405e+0008.82721629080e+0002.6 1.64887330889e+0008.39341353488e-0011.70390658776e+0003.08657579961e+0008.93195162776e+0002.8 1.63947266351e+0008.43741764569e-0011.68435909936e+0003.09292167559e+0009.04905972495e+0003.0 1.63203030387e+0008.63170285577e-0011.62038769698e+0003.09241637899e+0009.39453217083e+0003.2 1.62627371156e+0008.69897566347e-0011.59842096759e+0003.09472449777e+0009.52602920877e+0003.4 1.62260292300e+0008.77088688884e-0011.57694286086e+0003.09521665197e+0009.65317695129e+0003.6 1.62104289459e+0008.84765633406e-0011.55592302009e+0003.09395468087e+0009.77604354317e+0003.8 1.62162386648e+0008.92936706465e-0011.53536864498e+0003.09100776284e+0009.89450509805e+0004.0 1.62438155891e+0009.01595493705e-0011.51532001552e+0003.08645100005e+0001.00082431038e+0014.2 1.62935714722e+0009.10720246398e-0011.49584523459e+0003.08036324592e+0001.01167489922e+0014.4 1.63659702493e+0009.20273847590e-0011.47703429201e+0003.07282434712e+0001.02193368205e+0014.6 1.64615237898e+0009.30204536018e-0011.45899253832e+0003.06391202853e+0001.03151653186e+0014.8 1.65807861700e+0009.40447589839e-0011.44183362478e+0003.05369869065e+0001.04032708170e+0015.0 1.67241984210e+0009.03344994719e-0011.53038303379e+0003.05108284433e+0009.81308179470e+0005.2 1.68911306061e+0009-0011.50988213520e+0003.03691455618e+0009.90159651142e+0005.4 1.70824801577e+0009.29979382362e-0011.48227329348e+0003.02114506447e+0001.00341761251e+0015.6 1.72994237125e+0009.49635997655e-0011.44813478693e+0003.00419290127e+0001.02124401879e+0015.8 1.75435424888e+0009.74070357442e-0011.40816281857e+0002.98651746980e+0001.04384069653e+0016.0 1.78168659441e+0001.00380720990e+0001.36315240323e+0002.96861094212e+0001.07144487926e+0016.2 1.81219159789e+0001.03945293119e+0001.31397346313e+0002.95098687119e+0001.10432339885e+0016.4 1.84617465882e+0001.08168460110e+0001.26154381758e+0002.93416468883e+0001+0016.6 1.88399724950e+0001+0001.20680157287e+0002.91864983040e+0001+0016.8 1.92607788055e+0001+0001+0002.90490976875e+0001.23750358067e+0017.0 1.97323470415e+0001.32427622759e+0001.03534718100e+0002.90474041413e+0001.36931375097e+0017.2 2.02744672072e+0001.40025078086e+0009.86892009355e-0012.89626492151e+0001.42997226019e+0017.4 2.08747295883e+0001.47976212259e+0009.42712066644e-0012.88895660163e+0001.49074113864e+0017.6 2+0001.56260007661e+0009.02370127868e-0012.88261062733e+0001.55143109186e+0017.8 2.22621725695e+0001.64852115579e+0008.65491766933e-0012.87702888558e+0001.61182507517e+0018.0 2.30548684995e+0001.73726912188e+0008.31744928215e-0012.87202657228e+0001.67170138115e+0018.2 2.39167062745e+0001.82860921084e+0008.00821966465e-0012.86744151570e+0001.73086729372e+0018.4 2.48496772940e+0001.92237786577e+0007.72423666478e-0012.86314665927e+0001.78920431992e+0018.6 2.58531620491e+0002.01026889185e+0007.49002231821e-0012.85706715680e+0001.83905000696e+0018.8 2.69303002280e+0002.10878488517e+0007.24549531544e-0012.85319063934e+0001.89592151585e+0019.0 2.80836877473e+0002.21248106116e+0007.01072328195e-0012.85012072436e+0001.95456065112e+0019.2 2.93155998567e+0002.32220824682e+0006.78347551560e-0012.84801694276e+0002.01565392666e+0019.4 3.06283265420e+0002.43874511600e+0006.56231027963e-0012.84701318119e+0002.07979441462e+0019.6 3.20240826926e+0002.56272268664e+0006.34654094338e-0012.84719637152e+0002+0019.8 3.35048866424e+0002.69453247700e+0006-0012.84858471102e+0002.21876220069e+00110.0 3.50724107778e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+000I = 1.99597053717e+001的数值情况：3. 源程序1. 主函数#include#include#define length_t 51/xx向量长度#define size 51/xx向量长度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 uu=0;double table67=0;double ff=0;double ttlength_t=0;double yylength_t=0;double xxlength_t=0;double zzlength_t=0;double pplength_t=0;double qqlength_t=0;double sum=0.0;double q=0;double p=0;double x=0;double y=0;double z=0;double t=0;double h=0.2;double kk4=0;int r=0; int s=0;int ii,jj;int k_max=0;int m=0;double sum1=0;double sum2=0;double answer=0;FILE *fd;fd = fopen(tmptmp.txt, w);ii=0;jj=0;table00=7.3;table01=8.6;table02=10.1;table03=11.8;table04=13.7;table05=15.8;table05=18.1;table10=7.78;table11=9.08;table12=10.58;table13=12.28;table14=14.18;table15=16.28;table16=18.58;table20=9.22;table21=10.52;table22=12.02;table23=13.72;table24=15.62;table25=17.72;table26=20.02;table30=11.62;table31=10.92;table32=14.42;table33=16.12;table34=18.02;table35=20.12;table36=22.42;table40=14.98;table41=16.28;table42=17.78;table43=19.48;table44=21.38;table45=23.38;table46=25.78;table50=19.3;table51=20.6;table52=22.1;table53=23.8;table54=25.7;table55=27.8;table56=30.1;for(ii=0;iilength_t;ii+)ttii=h*ii;/yy0=2;for(ii=0;iisize-1;ii+)t=ttii;y=yyii;q=(double)interpolation(t,y,(double *)table); /输入一点，输出一点插值ite_equations(y,q,&x,&p,&z); /牛顿迭代解方程，输入一点（x,y）输出一点(t,u)xxii=x;zzii=z;ppii=p;qqii=q;kk0=(t*x*z-y*p)/(5*z*p);t=ttii+0.5*h;y=yyii+0.5*h*kk0;q=(double)interpolation(t,y,(double *)table); /输入一点，输出一点插值ite_equations(y,q,&x,&p,&z); /牛顿迭代解方程，输入一点（x,y）输出一点(t,u)kk1=(t*x*z-y*p)/(5*z*p);t=ttii+0.5*h;y=yyii+0.5*h*kk1;q=(double)interpolation(t,y,(double *)table); /输入一点，输出一点插值ite_equations(y,q,&x,&p,&z); /牛顿迭代解方程，输入一点（x,y）输出一点(t,u)kk2=(t*x*z-y*p)/(5*z*p);t=ttii+h;y=yyii+h*kk2;q=(double)interpolation(t,y,(double *)table); /输入一点，输出一点插值ite_equations(y,q,&x,&p,&z); /牛顿迭代解方程，输入一点（x,y）输出一点(t,u)kk3=(t*x*z-y*p)/(5*z*p);yyii+1=yyii+h/6*(kk0+2*kk1+2*kk2+kk3);m=25; sum1=0;for(ii=1;ii=m;ii+) sum1=sum1+yy2*ii-1;sum2=0;for(ii=1;iim;ii+) sum2=sum2+yy2*ii;answer=h/3*(yy0+yysize-1+4*sum1+2*sum2);fprintf(fd,t(i)y(i) x(i) z(i) p(i) q(i)n);for(ii=0;iisize;ii+)fprintf(fd,%-.1f %-.11e%-.11e%-.11e%-.11e%-.11en,ttii,yyii,xxii,zzii,ppii,qqii);fprintf(fd,nI = %.11e,answer);fclose(fd);2. Newton迭代求解非线性方程组#include#includevoid ite_equations(double y,double q,double *xx,double *pp,double *zz) /牛顿迭代解非线性方程组double x=1;double p=1;double z=1;double d3=0; /增量double F_Jacob33=0;double F3=0;double L33=0;double U33=0;double temp1=0;double temp2=0;int ii=0;int jj=0;int kk=0;for(ii=0;ii=1e2;ii+) F_Jacob00=8;F_Jacob01=-y;F_Jacob02=q*pow(2.718281828459,-z);F_Jacob10=-z;F_Jacob11=-q*pow(2.718281828459,-p);F_Jacob12=-x+0.6;F_Jacob20=-y*pow(2.718281828459,-x);F_Jacob21=4.7;F_Jacob22=-q;F0=8*x-q*pow(2.718281828459,-z)-p*y;F1=-x*z+0.6*z+q*pow(2.718281828459,-p);F2=y*pow(2.718281828459,-x)-q*z+4.7*p;LU(double *)L,(double *)U,(double *)F_Jacob,(double *)F);huidai(double *)L,(double *)U,(double *)F,(double *)d);x=x+d0;p=p+d1;z=z+d2; temp1=sqrt(x*x+p*p+z*z);temp2=sqrt(d0*d0+d1*d1+d2*d2);if(temp2/temp1)1e-12)*xx=x;*pp=p;*zz=z;break;3. LU分解求解线性方程组#include#include#define size 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_tempsizesize=0;double *b;double ysize=0;double xsize=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;iin;ii+) /对A_temp赋值for(jj=0;jjn;jj+)A_tempiijj=A(ii,jj);for(ii=0;iin;ii+)L(ii,ii)=1.0;for (kk=0;kkn;kk+) /总循环求出L U sum_answer=0;s_count=kk;for(pp=kk;ppn;pp+)/doolitte 选主元if (kk=0)sum_answer=0;elsefor(qq=0;qqs_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;ppn;pp+) /换位A_temptemp=A_temps_countpp; A_temps_countpp=A_tempkkpp;A_tempkkpp=temp;for(pp=0;ppkk;pp+) /换位Ltemp=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;iin;ii+) /先求一行Uif(kk=0)U(kk,ii)=A_tempkkii;elsesum_answer=0;for(qq=0;qqkk;qq+)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_tempkkii-sum_answer;for(ii=kk+1;iin;ii+) /求出一列Lsum_answer=0;for(qq=0;qqkk;qq+)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_tempiikk-sum_answer)/U(kk,kk);4. 2次回代解方程组（LU的基础上）#include#include#define size 3#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)void huidai(double *LL,double *UU,double *bb,double *xx) /回代法求xx, LUx=b double *L; double *U; double *b; double *x; double ysize=0; double sum_answer=0; int ii=0; int jj=0; L=LL; U=UU; b=bb; x=xx;/迭代Ly=by0=-b(0);for(ii=1;iisize;ii+