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

清华大学数学实验报告 实验4：常微分方程数值解 化学工程系 分2 安振华 2012011837 【实验目的】 1. 练习数值微分的计算； 2. 掌握用MATLAB 软件求微分方程初值问题数值解的方法； 3. 通过实例学习用微分方程模型解决简化的实际问题； 4. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式及稳定性等概念。 【实验内容】 【实验四：习题5】 射性废物的处理：有一段时间，美国原子能委员会（现为核管理委员会）处理浓缩放射性废物时，把它们装入密封性很好的圆桶中然后扔到水深300ft的大海中。这种做法是否会造成放射性污染，自然引起生物学家及社会各界的关注。原子能委员会一再保证，圆桶非常坚固，绝不会破漏，这种做法是绝对安全的。然而一些工程师们却对此表示怀疑，认为圆桶在和海底相撞时有可能发生破裂。于是双方展开了一场笔墨官司。 究竟谁的意见正确呢？原子能委员会使用的是55gal的圆桶，装满放射性废物时的圆桶重量为527.436 lbf,在海水中受到的浮力为470.327 lbf。此外，下沉时圆桶还要受到海水的阻力，阻力与下沉速度成正比，工程师们做了大量实验，测得其比例系数为0.08lbf s/ft。同时，大量破坏性试验发现当圆桶速度超过40ft/s时，就会因与海底冲撞而发生破裂。 （1） 建立解决上述问题的微分方程模型； （2） 用数值和解析两种方法求解微分方程，并回答谁赢了这场官司。【分析】 为了便于分析，首先，假设圆桶下沉过程无障碍且做直线运动，圆桶质量为m，装满放射性废物时的圆桶重量为M,下沉速度为v，极限速度为u，加速度为a, 重力加速度为g，海水深度为H，圆桶下沉时距海面距离为h.，比例系数为k。 G F f 受力分析：圆桶下落过程受重力G、浮力F、阻力f三个力的作用，受力情况如图所示 圆桶下落过程分析：由于圆桶刚开始时所受重力和浮力均不变，而所受阻力与速度成正比关系，速度增加时，阻力也随之增加，因此圆桶在此过程做加速度减小的加速直线运动。当阻力f增加到G-F-f=0，即加速度a=0时，阻力和速度同时达到最大，此后速度不变，即圆桶做匀速直线运动。 碰撞过程分析：通过以上分析我们可以猜想圆桶与海底碰撞时的速度有两种可能。第一，圆桶在与海底发生碰撞前就已达到最大速度，也就是说，此时圆桶是以最大速度与海底发生碰撞的；第二，圆桶有可能未达到最大速度就已海底发生碰撞。 【解答】 由已知条件换算得： M=527.436 lbf= 239.241kg， m=470.327lbf=213.337kg， g=9.800m/s^2， k=0.08 lbf.s/ft=0.119 kg.s/m， u=40ft/s=12.192m/s， H=300ft=91.440m 通过以上方程，我们可以计算出圆桶到达海底时的速度v，然后与极限速度u作比较，当vu时，圆桶与海底发生碰撞后破裂，当vu时，圆桶与海底发生碰撞后不会破裂。 构建函数： 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; %假设初速度为0 [t,v]=ode45(@fun1,ts,v0); %解微分方程（1） [t,v]; %输出结果 plot(t,v,b); %作出“时间-速度”图象 grid on; %加网格线 xlabel(t-时间); %横坐标表示时间 ylabel(v-速度); %纵坐标表示速度 title(圆桶下落速度随时间变化情况); 运行结果： 结果分析： 图像显示结果与前面分析大致相同，即圆桶在此过程做加速度减小的加速运动，直到当a=0时，v达到最大。 即 此时 这就是前面对圆桶与海底碰撞过程分析猜想的第一种可能情况，即圆桶有可能以