数 值 分 析 课程设计摸版

《数值分析课程设计》报告专业：x学号：x学生姓名：xxx指导教师：xx一、题目用阿当姆斯方法解常微分方程(用龙格—库塔公式给出出发值)1、'2010.101yyxhy2、'320.1010.101yxyxhy第一题的解第二题的解求微分方程'2(01)(0)1yyxy的精确解2dyydt21dydty公式两边分别对y和t求积分1ytc，又因为(0)1y，得到1c，所以11yt，t的取值范围为[0,1]。得到精确解：当t=0.1时，y=1.1111当t=0.2时，y=1.2500当t=0.3时，y=1.4286当t=0.4时，y=1.6667当t=0.5时，y=2当t=0.6时，y=2.5当t=0.7时，y=3.3334当t=0.8时，y=5当t=0.9时，y=10当t=1时，y取不到值二、理论设微分方程的初值问题为：'00(,)()yftyyty，将微分方程(,)dyftydt化为(,)dyftydt，然后再小区间1,kktt上对上式的两边进行积分，即11(,)kkkkytytdyftydt则微分方程就化为11(,)kktkktyyftydt。其中积分1(,)kkttftydt可以用数值积分的方法求得。通常采用不同的数值积分公式可以得到微分方程数值解的不同公式。如果对于被积函数(,)fty，取3kt，2kt，1kt，kt作为插值结点，则得到插值多项式333()()(,)()kkjiiikjkijjittPtftytt用插值多项式3()Pt的积分来近似代替(,())ftyt在区间1,kktt上的积分，即1113(,)()kkkkttkkkttyyftydtyPtdt最后得到如下公式1123(5559379)24kkkkkkhyyfff，这是一个显示阿当姆斯公式。如果取3kt，2kt，1kt，kt作为插值结点，则得到插值多项式11322()()(,)()kkjiiikjkijjittPtftytt代替(,())ftyt在区间1,kktt上的积分得到1112(9195)24kkkkkkhyyffff，这是隐式阿当姆斯公式。四阶阿当姆斯预报——校正系统。预报显示公式：1123ˆ(5559379)24nnnnnnhyyffff；校正隐式公式：1112(9195)24nnnnnnhyyffff；其中111ˆ(,)nnnffxy。阿当姆斯方法是线性多步法，当计算1ky时，需要用到前面四个点上的解函数值，所以用经典四阶龙格--库塔方法，给出初始值，公式如下：112341213243(22),6(,),(,),22(,),22(,).nnnnnnnnnnhyyKKKKKfxyhhKfxyKhhKfxyKKfxhyhK根据题目给出的第一个值0y，步长h,算出1y，2y，3y，代入阿当姆斯公式得到解。三、算法与程序设计#includestdio.h#includemath.h#includestdlib.hdoubletf(doublex,doubley){doubledy=y*y;returndy;}voidrkt4(doublet,doubleh,doubley[]){inti;for(i=0;i3;i++){doublek1,k2,k3,k4;k1=tf(t+i*h,y[i]);k2=tf(t+i*h+1.0/2*h,y[i]+1.0/2*h*k1);k3=tf(t+i*h+1.0/2*h,y[i]+1.0/2*h*k2);k4=tf(t+i*h+1.0/2*h,y[i]+h*k3);y[i+1]=y[i]+1.0/6*h*(k1+2*k2+2*k3+k4);}}voidadmas(doublet,doubleh,doubley[]){intm;doublex[11];for(m=0;m11;m++)x[m]=t+m*h;rkt4(t,h,y);for(m=0;m7;m++){y[m+4]=y[m+3]+1.0/24*h*(55*tf(x[m+3],y[m+3])-59*tf(x[m+2],y[m+2])+37*tf(x[m+1],y[m+1])-9*tf(x[m],y[m]));y[m+4]=y[m+3]+1.0/24*h*(9*tf(x[m+4],y[m+4])+19*tf(x[m+3],y[m+3])-5*tf(x[m+2],y[m+2])+tf(x[m+1],y[m+1]));}}voidmain(){intj;doublet,h,x[11],y[11];t=0.0;h=0.1;y[0]=1.0;admas(t,h,y);for(j=0;j11;j++){x[j]=t+j*h;printf(x=%e,x[j]);printf(\t\t);printf(y=%e,y[j]);printf(\n);}}其中函数tf的功能是得到f(x,y)的值，故在上述程序中只需修改函数tf中的内容，便可以得到不同的f(x,y)的值。题目第二题只需将函数tf改为：doubletf(doublex,doubley){doubledy=0.1*(x*x*x+y*y)returndy;}即可得到第二题的解四、第一题误差分析精确解与近似解的比较0x=001y0()1yx1x=0.111.11111y1()1.1111yx2x=0.221.25000y2()1.2500yx3x=0.331.42857y3()1.4286yx4x=0.441.66668y4()1.6667yx5x=0.552.00000y5()2yx6x=0.662.49977y6()2.5yx7x=0.773.33114y7()3.3334yx8x=0.884.97832y8()5yx9x=0.999.60200y9()10yx可以看到精确值与我们算得的值非常接近，四阶龙格——库塔公式的截断误差为O(5h)，即具有五阶精度。阿当姆斯预测校正系统也具有五阶精度。五、实际案例生物种群数量问题，设某生物种群在其适应的环境下生存，试讨论该生物种群的数量变化。问题假设：1、假设该生物种群的自然生长率为常数。2、设在其适应的环境下只有该生物种群生存或其他的生物种群的生存不影响该生物种群的生存。3、假设时刻t生物种群数量为N(t)，由于N(t)的数量很大，可视为时间t的连续可微函数。4、假设在t=0时刻该生物种群的数量为N0。得到模型0()()(0)dNtNtdtNN，再利用上述方法得到我们想要的值六、参考文献数值分析与算法主编徐士良机械工业出版社