数学实验4

实验4 常微分方程数值解实验目的：1、用MATLAB软件掌握求微分方程数值解的方法，并对结果作初步分析；2、通过实例学习用微分方程模型解决简化的实际问题。

【题目1】放射性废物的处理：有一段时间，美国原子能委员会（现为核管理委员会）处理浓缩放射性废物时，把它们装入密封性能很好的圆桶中，然后扔到水深300ft的大海中。

这种做法是否会造成放射性污染，自然引起生态学家及社会各界的关注。

原子能委员会一再保证圆桶非常坚固，绝不会破漏，这种做法是绝对安全的。

然而一些工程师们却对此表示怀疑，认为圆桶在海底相撞时有可能发生破裂。

于是双方展开了一场笔墨官司。

究竟谁的意见正确呢？原子能委员会使用的是55gal的圆桶，装满放射性废物时的圆桶重量为527.436lbf，在海水中受到的浮力为470.327lbf。

此外，下沉时圆桶还要受到水的阻力，阻力与下沉速度成正比，工程师们做了大量实验，测得其比例系数为0.08lbf s/ft。

同时，大量破坏性实验发现当圆桶速度超过40ft/s时，就会因为与海底冲撞而发生破裂。

（1）建立解决上述问题的微分方程模型。

（2）用数值和解析两种方法求解微分方程，并回答谁赢了这场官司。

【问题分析】本题应该求出圆桶下沉到海底时的速度，比较其与40ft/s的大小，来决定谁最后赢了这场官司。

已知水深h=300ft，装满废物的圆桶质量为m=527.436lbf，在海水中受到的浮力为B=470.327lbf，圆桶受到的海水阻力D与下降速度v成正比，比例系数K=0.08lbfs/ft。

设t时刻圆桶下落深度为y，圆桶速度为v，加速度为a，所受重力为G=mg，所受浮力为f，所受阻力为K*v,则圆桶所受合力为：F=G-f-Kv依据牛顿第二定律有：ma=m dvdt=md2ydt2=G−f−Kv=G−f−Kdydt由此可以利用常微分方程求解，其中v(0)=0,y(0)=0。

【程序运算】一、数值法记y(1)=y,y(2)=v,y=(y(1),y(2))T.编写如下M文件：function dy=fangshe(t,y)%建立名为fangshe的函数M文件m=527.436*0.4536;G=527.436*0.4536*9.8;f=470.327*0.4536*9.8;Kv=0.08*0.4536*9.8*y(2)/0.3048;dy=[y(2);(G-f-Kv)/m]%以向量性格是表示方程编写如下运行程序：>> ts=[0:2000];>> y0=[0,0];%给定初值>> opt=odeset('reltol',1e-3,'abstol',1e-6);%给定精确度>> [t,y]=ode45(@fangshe,ts,y0,opt);>> plot(t,y(:,2)，‘r’ )%画出v(t)-t图象>> gtext('t');>> gtext('v(t)');运行结果如下：>> plot(t,y(:,1))，axis([0 25 0 350]),grid; %画出y(t)-t图像>> hold on, ymax=300*0.3048;>> depth=linspace(ymax,ymax,2001);>> plot(t,depth,'b')，axis([0 25 0 350]),grid;%画出最深高度图像>>gtext('t'),>>gtext('y(t)'),>>gtext('ymax')运行结果如下：由图像知，t 大概等于13s 时，圆桶到达海底，结合以上两图可以得知，圆桶落入海底时，会因为速度过大而破裂。

二、解析法将模型分离变量得： mvdv mg −B −Cv =dy 对两边分别不定积分得：y =−vm C −gm 2−bm c 2ln (g −B m −Cv m )+C1 代入已知数据得：y =−v 4.87676×10−3− 1.06111（−4.87676）2ln (1+−4.87676×10−3v 1.06111) 当v=40ft/s 时，可利用上式算出y 值（即v 的一阶常微分方程）：齐次解的形式为v =Ae qt ，其中q=-4.87676×10−3 A =-217.59 ;特解为v=217.59 ；综上，v (t )=−217.59e −4.87676×10−3t +217.59 故有：y (t )=∫v (t )dt =∫(−217.59e −4.87676×10−3t +217.59)dt =91.44t 0t 0V(t)=13.8m/s=44.0ft/s【结果分析】结合题目的约束条件作图分析，判断下沉距离S=300ft时速度是否大于40ft/s或判断速度达到40ft/s时下沉的距离是否大于300ft（我选择用速度是否大于40ft/s判断），从而判断圆桶与海底冲撞时是否会破裂。

从依据数值解图的中可见，圆桶还没沉到海底时，其速度已达到规定的圆桶不破裂的最大速度。

由解析解算出的数据也说明，当下落高度达到300ft时，速度已达44.0ft/s超过规定速度，故圆桶会破裂。

工程师们赢了官司。

从速度-时间曲线可以看到，开始时圆桶的速度随着时间逐渐增大，加速度逐渐减小，可见圆桶做的是加速度逐渐减小的变加速运动。

时间充分大后（t>1500s）后，圆桶的加速度减至0，速度趋于一个定值，约为275m/s。

【题目2】一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。

已知河水流速为v1 与船在静水中的速度之比为k。

（1）建立描述小船航线的数学模型，求其解析解；（2）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。

（3）若流速v1=0，0.5，1.5，2（m/s），结果将如何？【问题分析】假设航行者并不知道水流速度。

设时间为t，小船坐标为（x,y）。

据题目所给公式，可以列如下方程：dx dt =v1−v x√x2+y2dy dt =v y√x2+y2初值为：（x,y）=(0,d) 【程序运算】（1）运用解析法将上两式相除得：dxdy =v1v2√(xy)2+1+xy令x/y=p<0, v1/v2=k有：dpdyy+p=k√p2+1+p整理得：√p2+1=k dyy求解得：ln(√p2+1+p)=k ln(cy)p=12(cy)−k−12(cy)kx=12(cy)−k y−12(cy)k y代入初值（x,y）=(0,d)可求得：c=−1 dfunction x=xiaochuan(y,k)x=0.5*(-0.01)^(-k)*y.^(1-k)-0.5*(-0.01)^k*y.^(1+k);%建立函数文件>> y=[0:-0.01:-100];>> x=xiaochuan(y,0.5);>> plot(x,y);%画出图像>> gtext('x/m')>> gtext('y/m')(2) 运用解析法记x(1)=x,x(2)=y,x=(x(1),x(2))T。

function dx=xiaochuan2(t,x,v1,v2)s=(x(1)^2+x(2)^2)^0.5;dx=[v1-x(1)/s*v2;-x(2)/s*v2]; %建立函数文件>> ts=[0:0.1:100];>> x0=[0,-100];%给定初始值>> opt=odeset('reltol',1e-3,'abstol',1e-9);%给定精确度>> [t,x]=ode15s(@xiaochuan2,ts,x0,opt,1,2);>> plot(t,x)%画出函数图像>>gtext('t/s')>>gtext('x,/m')>>gtext(‘y/m’)t/s x(t)/m y(t)/m0 0 -1000.1 0.099885 -99.80.2 0.199576 -99.60.3 0.29907 -99.40.4 0.398364 -99.20.5 0.497459 -990.6 0.596352 -98.80.7 0.695043 -98.60.8 0.793532 -98.40.9 0.891817 -98.21 0.989899 -981.1 1.087753 -97.8001·········66.6 0.23916 -0.0022666.7 0.139166 -0.0007666.8 0.039167 -6.00E-05从上图可以看出，当t=66.8s时，y接近0，故可认为小船到岸。

再运行如下程序：>> plot(x(:,1),x(:,2));>> xlabel('x/m')>> ylabel('y/m')（3）改变流速v1，得到以下几种情况：① v1=0，此时为静止的湖面。

显然船将沿直线直接从A 到B，渡河所需时间为100/2=50s② v1=0.5m/s 时，结果如下：跟v1=1 时的图形基本相似，仅在趋势上有微小不同。

在t=53.3 s时，y=-0.0099m可以认为船到达对岸。

③v1=1.5m/s时，运算结果如下：当t=114.0s时，y=0.000m，认为船到达对岸。

与v1=0.5，1m/s 时比较，小船几乎都是在[30s，40s]时间范围内到达x 方向的最大值，然后逆流驶向B 点。

v1 越小，小船在逆水回到B 的时间就越短。

④v1=2m/s 时，运算结果如下：可以发现曲线的变化趋势明显不同于前几个，且最终小船无法到达B 点，而到达了距离B点50 米的位置。

因为v1=v2，故k=1,将k代入x=12(cy)−k y−12(cy)k y得一抛物线，将 d=100,y=0 代入解得x=50，这与实验结果相符。

【结果分析】通过改变v1的值我们可以得到以下结论：1、在0<v1<v2 的情况下，随着v1 的增加，小船到达B 点的时间增加，与此同时，在整个过河的过程中，小船的最大水平距离也将增加。