《实验四 常微分方程的数值解综述》由会员分享，可在线阅读，更多相关《实验四 常微分方程的数值解综述（15页珍藏版）》请在金锄头文库上搜索。
1、实验四 常微分方程的数值解 化21华晓【实验目的】1. 掌握用MATLAB 软件求微分方程初值问题数值解的方法；2. 通过实例学习用微分方程模型解决简化的实际问题；3. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。【实验内容】五、概要：官司争论。将装满放射性废物的圆桶扔到水深300ft的海底，圆桶体积55gal，装满废料的桶重为527.436lbf，在海中浮力为470.327lbf。此外，下沉时受到的阻力与速度成正比，比例系数为0.08lbf/s。实验发现当圆桶速度超过40ft/s时，就会因与海底冲撞而破裂。（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求
2、解微分方程，并回答谁赢得了官司。建立模型：下落的圆桶在水中受到三个力，分别是重力，浮力和摩擦力。取向下为正方向，由牛顿第二定律可知：mg-F-f=ma 对于摩擦力，有：f=kv 同时，有：v=dxdta=d2xdt2得到微分方程：md2xdt2+kdxdt+F-mg=0对于该微分方程，有初始条件：x0=0x0=0其中F：圆桶所受浮力F=470.327*0.4536*9.82090.7Nv：圆桶速度，单位m/sK：阻力与速度的比例系数k=0.08*0.4536*9.8/0.30481.1667Ns/mm：圆桶质量 m=527.436*0.4536239.24kg可得，当圆桶速度超过12.192m
3、/s 时就会破裂；水深为91.44米。模型求解1. 数值求解首先编写微分方程function dx=yt(t,x)F=2092.11;m=239.41;k=1.1674;g=9.8;dx=x(2); -(k*x(2)+F-m*g)/m; 我们推测在20S内桶就已经沉到海底，现在我们来看看圆桶的速度与深度和时间的关系Matlab程序如下：ts=0:0.01:20;x0=0,0;t,x=ode45(yt,ts,x0);t,xsubplot(1,3,1),plot(ts,x(:,1)%画出深度关于时间的坐标图hold on ,z=91.44+0*ts;plot(ts,z ,r),hold off%在
4、此图的基础上增加一条显示深度为91.44m的线subplot(1,3,2),plot(ts,x(:,2)%画出速度关于时间的坐标图hold on ,z=12.192+0*ts;plot(ts,z,r),hold off%在此图的基础上增加一条显示极限速度的的线subplot(1,3,3),plot(x(:,1),x(:,2),grid; %画出速度与深度的关系图得到如下三个图形：取其中重要的数据如下：0000.10000.00530.10610.20000.02120.21220.30000.04770.31820.40000.08490.42410.50000.13260.5300 . .
5、.13.000087.822213.369913.100089.164213.469513.200090.516113.569013.300091.878013.668513.400093.249813.768013.500094.631613.867413.600096.023313.966713.700097.424914.066013.800098.836514.1653由数据可知，当圆桶到达海底时，速度已经接近13.6m/s,故会破裂。2.求解析解首先，我们直接求解这个微分方程初值问题对x作Laplace变换，设 Lxt=Fs 有ms2Fs+ksFs+F-mgs=0于是Fs=mg-Fs
6、(ms2+ks)=g-Fms+kms2作Laplace逆变换，得-F - mgmexpktm- m + ktk2代入数据解方程得v=13.6366m/s这个速度大于12.19m/s的临界速度，故圆桶会破裂。Matlab程序设计如下：syms v x %设定字母变量dsolve(D2x=(2344.6-2090.7-1.1667*(Dx)/239.24,x(0)=0,Dx(0)=0)得到结果为：x=217.6244*t+ 44625/exp(0.0048767*t)-44625因此可解得，当x=91.44时，t为13.2835s.而再由x的导数可以求得这个时刻速度为13.636m/s.已经超过极
7、限速度，故圆桶会撞击破裂。综上，数值解与解析解很好的符合，可得出圆桶会撞击海底而破裂，工程师们会赢得这场官司。六、一只小船度过宽为d的河流，目标是起点A正对着的B点，已知河水流速与船在静止的水中的速度之比为k。（1）建立描述小船航线的微分方程模型。（2）设d=100m,=1m/s,=2m/s，用数值方法求渡河所需时间，任意时刻小船位置及航行曲线，作图并与解析解比较；（3）若流速=0, 0.5, 1.5, 2（m/s），结果将如何？1模型建立：此种模型的前提是船并不事先知道水速，否则只要调整一个合适的角度，直接沿直线通过即可。现小船不知道水速，则它的策略应为始终使船头瞄准B点。对速度进行xy两个
8、方向的分解，可列出常微分方程如下：dxdt=v1-v2xx2+y2dydt=-v2yx2+y2其初值条件为x(0)=0,y(0)=-d.2问题分析：1解析解分析两式相除，得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数值求解首先建立微分方程：function dx=gh(t,x,v1,v2)d=100;v1=1;v2=2;s=(x(1)2+x(2)2)0
9、.5;dx=v1-x(1)/s*v2;-x(2)/s*v2;Matlab程序如下：ts=0:0.1:100; x0=0,-100; t,x=ode23s(gh,ts,x0);t xsubplot(1,3,1),plot(t,x(:,1),grid%画出x的坐标关于时间的关系图hold on ,z=0*t;plot(t,z,r),hold off%在此关系图上增加一条显示最终x的坐标的线subplot(1,3,2),plot(t,x(:,2),grid%画出小船y坐标关于时间的关系图hold on ,z=0*t;plot(t,z,r),hold off%在此图的基础上增加一条显示终点y坐标的线s
10、ubplot(1,3,3),plot(x(:,1),x(:,2),g),grid;%画出x y坐标的对应关系线得到如下图形：部分重要数据为：0000.01000.00010.01060.02000.00020.02120.03000.00050.03180.04000.00080.04250.05000.00130.0531.66.30000.3902-0.005966.40000.2903-0.003366.50000.1903-0.001466.60000.0903-0.000366.70000.0000066.80000.0000066.90000.0000067.0000-0.0000
11、0由数据可得，小船在66.7s时到达彼岸与解析解相比较：x与y的如下：x+12ckyk+1-12c-ky1-k=0其中c为未知常数，而k=v1v2另外由题目中已知可以解得k=0.5，再由初始条件x=0，y=-100，解得c=-0.01利用matlab做出相应图形，程序如下：y=-100:0.1:0;%给定初值k=0.5 ;x=-y.*(-0.5*(-0.01*y).k)+0.5*(-0.01*y).(-k);%给出相应的函数关系式plot(x,y),gridxlabel(x);ylabel(y);得到图形为：发现数值解与解析解图形重合度极好，故可得66.7s到达彼岸。(一)当V1=0m/s时,
12、function dx=gh(t,x,v1,v2)d=100;v1=0;v2=2;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;t,x=ode15s(gh,ts,x0);plot(t,x(:,1),r),gridhold on,plot(t,x(:,2),g),hold offgtext(x(t);gtext(y(t);pauseplot(x(:,1),x(:,2),g),grid;xlabel(x);ylabel(y);得到图形为：重要数据为： 0 0100.0000 0.1000 099.8000
13、0.2000 099.6000 0.3000 099.4000. 49.6000 00.8000 49.7000 00.6000 49.8000 00.4000 49.9000 00.2000 50.0000 00.0000故知，小船会在50s时到达彼岸。与解析解相比较：由数学计算可知：此时t=d/v2=100/2=50s结果一致。（二）当V1=0.5m/s时,主程序不变，只把gh.m脚本文件中v1数值改变，即得如下数据：00-1000.10.049947-99.80.20.099796-99.60.30.149545-99.40.40.199194-99.20.50.248742-99 . . .52.90.604925-0.30225530.474462-0.2161853.10.340737-0.1372453.20.203201-0.0678453.30.060869-0.01339得到图形为：水速为0.5m/s时，时间为53.3s时船的距河岸0.060869m，可以认为船到了河对岸
