试验4常微分方程的数值解.doc

实验5 常微分方程的数值解 概要：将装满放射性废物的圆桶扔到水深300ft的海底，圆桶体积55gal，装满废料的桶重为527.436lbf，在海中浮力为470.327lbf。此外，下沉时受到的阻力与速度成正比，比例系数为0.08lbf/s。实验发现当圆桶速度超过40ft/s时，就会因与海底冲撞而破裂。 要求：（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求解微分方程，并回答谁赢得了官司。 模型建立 由牛顿第二定律可列出圆桶下沉速度的微分方程： 其中G为圆桶重量，F为浮力，b为下沉阻力与速度关系的比例系数。换算到国际单位制，dept=300*0.3048=91.4400 海深（m） ve=40*0.3048=12.1920 速度极限（超过就会使圆筒碰撞破裂）(m/s) G=527.436*0.4536*9.8＝2344.6 圆筒重量（N） F=470.327*0.4536*9.8＝2090.7 浮力(N) m=527.436*0.4536＝239.24 圆筒质量（kg） b=0.08*0.4536*9.8/0.3048＝1.1667 比例系数(Ns/m) 模型求解 一．求数值解 Matlab程序如下： sd.m: function dx=sd(t,x,G,F,m,b) dx=[(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阶龙格—库塔公式计算 %[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.2 y=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=300.*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)); ? 求解结果如下图： 速度—时间曲线： 可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常值。那时重力，浮力和阻力达到平衡。 对速度积分可得时间－深度曲线（估计20s内桶会落到300ft的海底，故图中时间轴上数值到20）： 再画出速度－时间曲线，同时在图中标出了海深dept和速度极限值Vmax 从图中可以观察到圆桶落到海底时（即dept（v）＝dept时），速度已经超过Vmax，故圆桶会因为速度过快而与海底碰撞破裂。为求处此速度，将上图局部放大： 可以看到当圆桶沉到海底时速度为13.6344m/s，超过了12.1920的极限速度。 ? 二．解析解 用Mathematica求解模型的微分方程： 对v(t)积分得到下沉深度de(t)= 44624.996*exp(-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),所以圆桶在海底将会