实验八常微分方程初值问题数值解法报告.docVIP

实验八 常微分方程初值问题数值解法 一、基本题 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题： （1） 分别取h=0.1,0.2,0.4时数值解。 初值问题的精确解。 （2） 用r=3的Adams显式和预 - 校式求解 取步长h=0.1，用四阶标准R-K方法求值。 （3） 用改进Euler法或四阶标准R-K方法求解 取步长0.01,计算数值解，参考结果 。 （4）利用四阶标准R- K方法求二阶方程初值问题的数值解 （I） (II) (III) (IV) ? 二、应用题 1. 小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2. 建立火箭升空过程的数学模型（微分方程）； 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度． 2. 小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k，火箭升空过程的数学模型为 其中为火箭在时刻t的高度，m=1200-15t为火箭在时刻t的质量，T（=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t1 (=900/15=60秒)为引擎关闭时刻． 今测得一组数据如下（t~时间（秒），x ~高度（米），v~速度（米/秒））： t 10 11 12 13 14 15 16 17 18 19 20 x 1070 1270 1480 1700 1910 2140 2360 2600 2830 3070 3310 v 190 200 210 216 225 228 231 234 239 240 246 现有两种估计比例系数k的方法： 1．用每一个数据(t,x,v)计算一个k的估计值（共11个），再用它们来估计k。 2．用这组数据拟合一个k　． 请你分别用这两种方法给出k的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）． 三、要求 1、 根据初值问题数值算法，分别选择二个初值问题编程计算； 2、 试分别取不同步长，考察某节点处数值解的误差变化情况； 3、 试用不同算法求解某初值问题，结果有何异常； 4、 分析各个算法的优缺点。 四、目的和意义 1、 熟悉各种初值问题的算法，编出算法程序； 2、 明确各种算法的精度与所选步长有密切关系； 3、 通过计算更加了解各种算法的优越性。 五、实验学时：2学时 六、实验步骤： 1．进入C或matlab开发环境； 2．根据实验内容和要求编写程序； 3．调试程序； 4．运行程序； 5．撰写报告,讨论分析实验结果. 解： 1、欧拉法程序 function y=Euler(a,b,M,y0) %a=1,b=2,M=10,f=t*y^(1/3),y0=1; h=(b-a)/M; t=zeros(1,M+1); t=a:h:b; y=zeros(1,M+1); yy=zeros(1,M+1); y(1)=y0; for k=1:M y(k+1)=y(k)+h*t(k)*y(k)^(1/3); end yb=y(M+1); yy=((t.^2+2)./3).^1.5; det=yy-y; plot(t,y,r-,t,yy,b:,t,det); 改进欧拉法程序 function H=heeuler(a,b,M,ya,f) %a=0,b=1,M=10,f=t*t+t-y,y0=0; h=(b-a)/M; t=zeros(1,M+1); y=zeros(1,M+1); p=0;q=0; t=a:h:b; y(1)=ya; for k=1:M p=feval(f,t(k),y(k)); q=feval(f,t(k+1),y(k)+h*p); y(k+1)=y(k)+0.5*h*(p+q); end yy=t.*t-t+1-exp(-t); det=yy-y; plot(t,y,r-,t,yy,b:,t,det); H=[t,y,yy,det] function f=ff(t,y); f=t.^2+t-y; 3、四阶龙格-库塔法