《实验5常微分方程的数值解》由会员分享，可在线阅读，更多相关《实验5常微分方程的数值解（27页珍藏版）》请在金锄头文库上搜索。
1、实验实验 5 常微分方程的数值解常微分方程的数值解概要：将装满放射性废物的圆桶扔到水深 300ft 的海底，圆桶体积 55gal，装满废料的桶重 为 527.436lbf，在海中浮力为 470.327lbf。此外，下沉时受到的阻力与速度成正比，比例系 数为 0.08lbf/s。实验发现当圆桶速度超过 40ft/s 时，就会因与海底冲撞而破裂。 要求：（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求解微分方程， 并回答谁赢得了官司。 模型建立模型建立 由牛顿第二定律可列出圆桶下沉速度的微分方程： dvGFfGFbv dtmm其中 G 为圆桶重量，F 为浮力，b 为下沉阻力与速度关
2、系的比例系数。换算到国际单 位制，dept=300*0.3048=91.4400 海深（m） ve=40*0.3048=12.1920 速度极限（超过就会使圆筒碰撞破裂） (m/s) G=527.436*0.4536*9.82344.6 圆筒重量（N） F=470.327*0.4536*9.82090.7 浮力(N) m=527.436*0.4536239.24 圆筒质量（kg） b=0.08*0.4536*9.8/0.30481.1667 比例系数(Ns/m) 模型求解模型求解 一求数值解 Matlab 程序如下： sd.m: function dx=sd(t,x,G,F,m,b) dx=(
3、G-F-b*x)/m;%微分方程sddraw.m: clear; G=527.436*0.4536*9.8;%圆筒重量（N） F=470.327*0.4536*9.8;%浮力(N) m=527.436*0.4536;%圆筒质量（kg） b=0.08*0.4536*9.8/0.3048%比例系数(Ns/m) h=0.1;%所取时间点间隔 ts=0:h:2000;%粗略估计到时间 2000 x0=0;%初始条件 opt=odeset(reltol,1e-3,abstol,1e-6);%相对误差 1e-6，绝对误差 1e-9 t,x=ode45(sd,ts,x0,opt,G,F,m,b);%使用 5
4、 级 4 阶龙格库塔公式计算 %t,x%输出 t,x(t),y(t) plot(t,x,-),grid%输出 v(t)的图形 xlabel(t); ylabel(v(t);%用辛普森公式对速度积分求出下沉深度 T=20;%估计 20s 以内降到海底 for i=0:2:10*T%作图时间间隔为 0.2y=x(1:(i+1);k=length(y);a1=y(2:2:k-1);s1=sum(a1);a2=y(3:2:k-1);s2=sum(a2);z4(i+2)/2)=(y(1)+y(k)+4*s1+2*s2)*h/3;%辛普森公式求深度 end i=0:2:10*T;figure; de=30
5、0.*0.3048.*ones(5*T+1,1);%海深 ve=40.*0.3048*1 1;%速度极限值（超过就会使圆筒碰撞破裂） plot(x(i+1),z4,x(i+1),de,ve,0 z4(5*T+1);%作出速度深度图线，同时画出海深和速度要求 grid; gtext(dept),gtext(Vmax); xlabel(v); ylabel(dept(v);figure; plot(i/10,z4);%作出时间下降深度曲线 grid; xlabel(t); ylabel(dept(t);求解结果如下图： 速度时间曲线：可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常
6、值。那时重 力，浮力和阻力达到平衡。 对速度积分可得时间深度曲线（估计 20s 内桶会落到 300ft 的海底，故图中时间轴上 数值到 20）：再画出速度时间曲线，同时在图中标出了海深 dept 和速度极限值 Vmax从图中可以观察到圆桶落到海底时（即 dept（v）dept 时） ，速度已经超过 Vmax， 故圆桶会因为速度过快而与海底碰撞破裂。为求处此速度，将上图局部放大：可以看到当圆桶沉到海底时速度为 13.6344m/s，超过了 12.1920 的极限速度。二解析解 用 Mathematica 求解模型的微分方程：对 v(t)积分得到下沉深度 de(t)= 44624.996*exp(
7、-0.004877*t)+217.6224*t-44624.996,令 de(t)91.4400 得 t=13.2812,代入 v(t)得 v(13.2812)= 13.6507（m/s）Vmax=12.1920(m/s),所 以圆桶在海底将会碰撞破裂。 解析解 13.6507 与数值解 13.6344 相差甚小，很好地验证了数值解得正确性。结论结论：经过计算表明圆桶到达海底时的速度将超过 40ft/s，即圆桶会破裂。工程师们将赢得 官司。第五题第五题已知已知：一只小船度过宽为 d 的河流，目标是起点 A 正对着的 B 点，已知河水流速1v与船在静止的水中的速度2v之比为 k。求求：（1）建立
8、描述小船航线的微分方程模型。 （2）设 d=100m, 1v=1m/s, 2v=2m/s，用数值 方法求渡河所需时间，任意时刻小船位置及航行曲线，作图并与解析解比较；（3）若流速1v=0, 0.5, 1.5, 2（m/s） ，结果将如何？ 模型建立：模型建立：如图以 B 点为原点建立直角坐标系。假设驾驶船的人不知道水流速度（如果知道那可 以走直线从 A 到 B 点） ，他驾船使船头的方向始终对着目标 B 点，如图，可得速度 V 的 x，y 方向的分量为：2 122222v xdxvdtxyv ydy dtxy 这是一个微分方程模型，初始条件为( , )(0, 100)x y 。 模型求解：模型
9、求解： 用龙格库塔方法求解此微分方程，程序如下： gh.m: function dx=gh(t,x,v1,v2) s=(x(1)2+x(2)2)0.5; dx=v1-x(1)/s*v2;-x(2)/s*v2;%以向量形式表示微分方程ghdraw.m:h=0.01;%所取时间点间隔 ts=0:h:100;%粗略估计到达目标点时间在 100 以内 x0=0,-100;%初始条件 opt=odeset(reltol,1e-6,abstol,1e-9);%相对误差 1e-6，绝对误差 1e-9 t,x=ode15s(gh,ts,x0,opt,1,2);%使用解刚性方程得龙格库塔公式计算,1,2 是给
10、gh 函数的 参数 t,x%输出 t,x(t),y(t) plot(t,x,-),grid%输出 x(t),y(t)的图形 gtext(x(t),gtext(y(t),pause plot(x(:,1),x(:,2),-),grid,%作 y(x)的图形 gtext(x),gtext(y); 说明：说明：1由于 v2-v1=1,d=100 所以过河时间一定在 100 以内2一开始我使用了 ode45()公式求解，但计算时间太长，后来使用了求解刚性方程 的函数 ode15s()，计算速度较快，可以推断本模型的微分方程是刚性刚性的。 t,x命令输出的 t, x(t), y(t) 有 10001 组
11、，不可能在此一一列举，并且从图上可以观察 数据变化趋势，在此列出开始和最后的一些数据：有此可以看出到 t=66.66(s)时，x(t)=0.0070(m),y(t)=0(m)，可以认为船已经到达目标。 下面给出 x(t),y(t)，和 y(x)的图形：从上图可以看出船的轨迹形状。 为了看出船在轨迹上各点的情况，将时间间隔 h 改成 0.5 重新画图：图中每一个圆圈标记的时间间隔为 0.5，圆圈密集的地方船的速度较慢，圆圈稀疏的地圆圈密集的地方船的速度较慢，圆圈稀疏的地 方船速较快方船速较快。可以预测船头方向与水流方向夹角越小时，船速越快，这在图中有所反映。结合解析解的讨论结合解析解的讨论：由本
12、题的微分方程表达式：2 122222v xdxvdtxyv ydy dtxy 可以推出 y-x 的解析表达式，如下：1111022kkkkxc ycy，其中 c 为待定常数，12vkv 为速度比。代入初值条件000,100xy 得0.01c ，所以有1111( 0.01)( 0.01)022kkkkxyy ，可以表示为 x(y)的显示表达式：1111( 0.01)( 0.01)22kkkkxyy 。 用 matlab 画出此函数的曲线，程序如下： xy.m: function x=f(y) k=0.5; x=-0.5.*(-0.01).k.*y.(k+1)+0.5.*(-0.01).(-k).
13、*y.(-k+1);xyplot.m： clear; y=0:-0.1:-100; for i=0:1:1000; x(:,i+1)=xy(-i/10); end plot(x,y); grid; gtext(x);gtext(y);结果如下：与上面数值解画出的曲线相同。其他情况讨论其他情况讨论（为方便起见，下面均为时间间隔为 0.5 的结果）：（1）流速1v=0(m/s)时 这时水是静止的，船延直线开到 B 点，所用时间为21/()100/(2 1)50dvv（s） 程序运行结果如下(由 ode23s 算得)： t,x:可以看到 t=50 时，到达 B 点，与解析解相符。 下面给出 x(t)
14、,y(t)，和 y(x)的图形：（2）流速1v=0.5(m/s)时 程序运行结果如下(由 ode15s 算得)： t,x:可以看到 t53 时，到达 B 点。下面给出 x(t),y(t)，和 y(x)的图形：图形和 V1=1 时十分相似。（3）流速1v=1.5(m/s)时 程序运行结果如下(由 ode15s 算得)：t,x:可以看到 t114 时，到达 B 点。 下面给出 x(t),y(t)，和 y(x)的图形：右图可见由于河水速度的上升，船的轨迹绕的弯增大了。（3）流速1v=2(m/s)时 预测：由于12vv，所以2 112222(1)0,(0)v xdxxvvxdtxyxy ，122 0(
15、1)0,0TxvdtT xy ，船不可能到达正对岸的 B 点。 程序运行结果如下(由 ode15s 算得)：t,x:可以看到船无法到达目标点，经过足够长时间后会停在（50，0）处 下面给出 x(t),y(t)，和 y(x)的图形：这个结果与预测相统一，船最后停在（50，0）处。第第 8 题题 问题：问题：两种群相互竞争模型如下：11 1222 21(1)(1)dxxyrxsdtnndyyxr ysdtnn其中 x(t),y(t)分别为甲乙两种群的数量，1r，2r为它们的固有增长率，1n，2n为它们的最大容量。1s的含义是，对于供养甲的资源来说，单位数量的乙（相对2n）的消耗为单位数量甲（相对1n）消耗的1s倍，对2s可以作相应解释。 该模型无解析解，试用数值方法研究以下问题： （1）（1） 设 r1=r2=1,n1=n2=100,s1=0.5,s2=2,初值 x0=y0=10,计算 x(t),y(t),画 出它们的图形及图（x,y） ，说明时间 t 充分大了以后 x(t),y(t)的变化趋势（人 们今天看到的已经是自然界长期演变的结局） 。 （2）（2） 改变 r
