数学实验4-常微分方程数值解-安振华-2012011837

1、清华大学数学实验报告实验4：常微分方程数值解化学工程系 分2 安振华 2012011837【实验目的】1. 练习数值微分的计算；2. 掌握用MATLAB 软件求微分方程初值问题数值解的方法；3. 通过实例学习用微分方程模型解决简化的实际问题；4. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式及稳定性等概念。【实验内容】【实验四：习题5】射性废物的处理：有一段时间，美国原子能委员会（现为核管理委员会）处理浓缩放射性废物时，把它们装入密封性很好的圆桶中然后扔到水深300ft的大海中。这种做法是否会造成放射性污染，自然引起生物学家及社会各界的关注。原子能委员会一再保证，圆桶非常坚固，绝不会破漏，

2、这种做法是绝对安全的。然而一些工程师们却对此表示怀疑，认为圆桶在和海底相撞时有可能发生破裂。于是双方展开了一场笔墨官司。究竟谁的意见正确呢？原子能委员会使用的是55gal的圆桶，装满放射性废物时的圆桶重量为527.436 lbf,在海水中受到的浮力为470.327 lbf。此外，下沉时圆桶还要受到海水的阻力，阻力与下沉速度成正比，工程师们做了大量实验，测得其比例系数为0.08lbf s/ft。同时，大量破坏性试验发现当圆桶速度超过40ft/s时，就会因与海底冲撞而发生破裂。（1）建立解决上述问题的微分方程模型；（2）用数值和解析两种方法求解微分方程，并回答谁赢了这场官司。【分析】为了便于分析，

3、首先，假设圆桶下沉过程无障碍且做直线运动，圆桶质量为m，装满放射性废物时的圆桶重量为M,下沉速度为v，极限速度为u，加速度为a, 重力加速度为g，海水深度为H，圆桶下沉时距海面距离为h.，比例系数为k。GFf受力分析：圆桶下落过程受重力G、浮力F、阻力f三个力的作用，受力情况如图所示圆桶下落过程分析：由于圆桶刚开始时所受重力和浮力均不变，而所受阻力与速度成正比关系，速度增加时，阻力也随之增加，因此圆桶在此过程做加速度减小的加速直线运动。当阻力f增加到G-F-f=0，即加速度a=0时，阻力和速度同时达到最大，此后速度不变，即圆桶做匀速直线运动。碰撞过程分析：通过以上分析我们可以猜想圆桶与海底碰撞

4、时的速度有两种可能。第一，圆桶在与海底发生碰撞前就已达到最大速度，也就是说，此时圆桶是以最大速度与海底发生碰撞的；第二，圆桶有可能未达到最大速度就已海底发生碰撞。【解答】由已知条件换算得：M=527.436 lbf= 239.241kg，m=470.327lbf=213.337kg，g=9.800m/s2，k=0.08 lbf.s/ft=0.119 kg.s/m，u=40ft/s=12.192m/s，H=300ft=91.440m 通过以上方程，我们可以计算出圆桶到达海底时的速度v，然后与极限速度u作比较，当v>u时，圆桶与海底发生碰撞后破裂，当v<u时，圆桶与海底发生碰撞后不会破

5、裂。构建函数：function dv=fun1(t,v) %构造“时间-速度”函数M=239.241;m=213.337;g=9.8;k=0.119; %根据题意选取模型参数dv=(M-m-k*v)*g/M; %表示微分方程（1）MATLAB程序:clf ;ts=linspace(0,1500,15000); %考察圆桶下落25分钟内速度变化情况v0=0; %假设初速度为0t,v=ode45(fun1,ts,v0); %解微分方程（1）t,v; %输出结果plot(t,v,'b'); %作出“时间-速度”图象grid on; %加网格线xlabel('t-时间'

6、); %横坐标表示时间ylabel('v-速度'); %纵坐标表示速度title('圆桶下落速度随时间变化情况');运行结果：结果分析：图像显示结果与前面分析大致相同，即圆桶在此过程做加速度减小的加速运动，直到当a=0时，v达到最大。即 此时 这就是前面对圆桶与海底碰撞过程分析猜想的第一种可能情况，即圆桶有可能以最大速度与海底发生碰撞，此时圆桶的速度远远大于极限速度，圆桶与海底发生碰撞后会发生破裂。但由于该计算过程并未考虑海底深度，也就是说圆桶有可能并未达到最大速度就与海底发生碰撞，那么在这种情况下圆通是否能够达到极限速度还无法判断，因此还必须分析速度与圆桶下沉

7、距离的关系，求出圆桶下落到海底时的速度，然后与极限速度作比较，进而得出结论。clf;ts=linspace(0,300,3000); %考察圆桶下落5分钟内速度变化情况v0=0; %假设初速度为0t,v=ode45(fun1,ts,v0); %解微分方程（1）t,v; %输出结果h=cumtrapz(t,v); %利用梯形积分公式求圆桶距海面距离H=91.44;u=12.192; %根据题意选取模型参数plot(h,v,'r'); %作出“速度-距离”图象hold on;plot(0,100,12.192,12.192,'-g'); %设定水平线-表示极限速度p

8、lot(91.44,91.44,0,15,'-b'); %设定铅垂线-表示海底grid on; axis(0,100,0,15); %设定图像范围text(91.44,6,'海底'); %标注“海底”的位置text(35,12.192,'极限速度'); %标注“极限速度”的位置xlabel('h-圆桶到海面的距离'); ylabel('v-速度'); title('速度随距离变化情况');MATLAB程序:运行结果：结果分析：由图像可知，到达海底时，圆通速度已经超过了极限速度，因此，圆桶与海底碰撞后

9、会发生破裂，里面的放射性废物会泄露。 从而解得这与数值法求解微分方程的结果相一致。另外,由（4）式，令h(t)=91.440m，可解出t=13.2697s，代入（3）式解得，v=13.6348m/s>12.192m/s，即圆桶到达海底时的速度大于极限速度，因此圆桶与海底碰撞后会发生破裂。由上面的分析可知，工程师们的结论是正确的。【结论】通过以上两种方法得出结论：圆桶到达海底时的速度大于极限速度，因此圆桶与海底碰撞后会发生破裂，所以工程师将赢得这场官司。 【实验四：习题6】一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。已知河水流速v1 与船在静水中的速度v2之比为k。（1

10、） 建立描述小船航线的数学模型，求其解析解；（2） 设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。（3） 若流速v1=0，0.5，1.5，2（m/s），结果将如何？【分析】小船运动速度矢量分解图如下：vV2V1根据实际情况，分析可得到较为符合实际的数学模型：船并不知道水速，否则只要调整合适角度，即可沿直线通过即可。现小船不知道水速，应该通过船的控制系统对速度进行xy两个方向的分解，可列出常微分方程如下：dxdt=v1-v2xx2+y2dydt=-v2yx2+y2其初值条件为x(0)=0,y(0)=-d.1、常微分

11、方程解析解两式相除得dxdy=-v1x2+y2v2y+xy=-kxy2+12+xy令u=xy 则有dxdy=ydudy+u代入ydudy=-ku2+12分离变量积分可得y-k=c(u2+12+u)代入初值，可得c=-d-k整理得x=-d2(y-d)1-k-(y-d)1+k2、 MATLAB数值解【解答】用matlab求上述方程组的数值解。记x(1)=x, x(2)=y, x=(x(1), x(2)T。编写M文件如下：function dx=guohe(t,x)v1=1;v2=2;s=(x(1)2+x(2)2)0.5;dx=v1-x(1)*v2/s;-x(2)*v2/s;endMATLAB程序：

12、x0=0,-100;ts=0:0.1:70; %经过试验设定终值为70t,x=ode15s(guohe,ts,x0); %求解方程组plot(t,x),grid; %画图同时做出x、y方向的位移title('图7.分位移-时间')xlabel('t/s');ylabel('x,y/m');gtext('x(t)');gtext('y(t)');pause;plot(x(:,1),x(:,2),grid; %画出航行曲线title('航行曲线');xlabel('x');ylabel(

13、'y');t,x(:,1),x(:,2) %输出数据表运行结果： 由图7可以看出，在前30s中，y(t)有较好的的线性，因此，在这个时间段内，小船沿y轴方向的为匀速。从图8可以看出，船合速度的方向为轨迹的切线方向。在刚开始阶段，可以认为传顺水而行，小船在刚开始的速度方向的确变化不大，与上图相吻合。小船在达到x最大值时开始转向，着这一时刻附近的合速度最小。部分数据：t/sx/my/m00-100.00000.10000.0999-99.80000.20000.1996-99.60000.30000.2991-99.40000.40000.3984-99.20000.50000.4

14、975-99.00000.60000.5964-98.80000.70000.6951-98.60000.80000.7936-98.40000.90000.8919-98.20001.00000.9900-98.00001.10001.0879-97.80001.20001.1856-97.60001.30001.2831-97.40001.40001.3804-97.20001.50001.4775-97.00011.60001.5744-96.80011.70001.6710-96.60011.80001.7675-96.400165.10001.5882-0.098765.20001.

15、4886-0.086665.30001.3889-0.075465.40001.2892-0.064965.50001.1894-0.055265.60001.0896-0.046365.70000.9898-0.038265.80000.8899-0.030865.90000.7900-0.024366.00000.6901-0.018566.10000.5902-0.013566.20000.4902-0.009366.30000.3902-0.005966.40000.2903-0.003366.50000.1903-0.001466.60000.0903-0.000366.70000.

16、0000066.80000.00000解析解与数值解的比较MATLAB程序：y0=0:-1:-100;x0=50.*(y0./(-100).0.5-(y0./(-100).1.5);plot(x(:,1),x(:,2),'r',x0,y0,'b'),grid;legend('数值解','解析解')title('航行曲线比较');xlabel('x');ylabel('y');运行结果：由图9可以看出，数值解得到的曲线与解析解的到的曲线几乎完全重合，说明数值解得到的结果还是很准确的。改

17、变流速（1）v1=0m/s时,更改程序中的数值，可得如下分位移-时间图像及航行曲线：由图像可知，v1=0m/s时，水的流速为0，即小船的合速度等于v2，因此小船沿y轴方向匀速前进。此时小船的位置的解析解为x=0.5*y-0.5*y=0，所以小船过河轨迹为沿y轴方向的一条直线。小船沿y方向的速度dydt=-v2yx2+y2=-v2y|y|=v2，故小船匀速过河。过河所用时间为50s。（2）v1=0.5s时，由数据表可知，t=53.3s时，y=-0.0131,可近似认为已到达岸边。将V1=0.5m/s与v1=1m/s时相比，河水流速减小后，船在x方向上的分速度减小了，因此小船在x方向上行驶的最大值

18、减少。这与途中x最大值减小了一半左右相符。与v1=0时比较，小船过河花的时间并没有增加很多，但路程长了很多，这说明小船的平均合速度大于v2。这是因为实际速度为河水与船水速度的矢量叠加结果。（3）v1=1.5m/s时，由数据表可知，t=114.3s时，x=0.0007m,y=0.000m，可认为小船到达对岸。从图中可以看出，小船在30s，40s时间范围内达到x方向的最大值，然后逆流驶向B点。所以v1越小，船回到B的时间就越短，这一点与常识相符。将解析式x=d2(y-d)1-k-(y-d)1+k求导并令其为零，可以求出x在y=-100*2k1-k1+k处，取得最大值，所以k在0,1区间内越大，x取得最大值的点就越靠近x轴。（4）v1=2m/s时利用程序可以求出，小船最终在下游50m处靠岸，不能到达正对岸。也就是说，无论多长时间小船都无法到达正对岸。由解析解也可以证明这一点。dxdt=v1-