实验二常微分方程的数值解
土木系结23李会平
【实验目的】
1、掌握用matlab解数值微分方程
2、了解龙格-库塔方法的基本原理
3、用这些手段解决一些实际的问题
【实验内容】
4-5核废料问题
首先列出问题的运动学方程，由牛顿第二定律，
md 2s/dt2=G-F-f
其中m=G/g, f=kv,由于我们熟悉的单位是公制的，所以在定义函数的时候进行了单位
436*0.45^6*9.8;
F=470.327*0.4536*9,8
jn-527. 136*0. 1536.
k=0. 081*0. 4536^9. 3/0. 3fl46.
I d^Kt-F-k^djJ/M xd)]; 谥先应该走又G弄宾里，苦咄益比错;貞"就義噬度，代表潇償
这其实是一个关于s的二阶常微分方程，需要定义两个变量x(1),x(2)将其化为一阶微分方程组，此处x(1),x(2)实际上分别代表速度和深度，相应的微分方程如代码中所示:
dx=[(G-F-k*x(1))/m;x(1)]
(如课堂提醒，中间应该是分号，这点容易出错。
) 接写来进行m文件的命令编写，如下所示:
Ejditci - C ;\U SET s\ I enovo\Deii ktc p ,ri atl J b te a di er^i I m rvttr k\ho TI r ■ k2\Fuiig,rrFeiwaE
- t=0:2000 :
2- ^[3,0].
S—[t, nJ =D de4 S e IW U L t, irO):
4- pl nt (t, i f;, 1)/ b J)f frid;
5- hold on：
C - v>ax-40*0. Jilts (leiiEth^t), L ■.
；- plilt (t, VMI, 1J> ；可画出最犬逋度允许望
G - hold cff
9 - xlabelCt (i)J/ fGntEi:e'>16)
it - ylab^ljc v(n/s''.L toutsiz?L,.16?
H - t"l亡广速度的吏北观韋
12- P3US 9
13- pint <*(!, 1), a |!, 2) / b?grid：雷画二募縻肖逋度的美童
14- hold on;
15- 5»ax^30Q<Q» SQ^S^cnesd^nffthCx 1)), 1 ; $—走聲特常数转'七为数绢，'才肖作用！
1C - plal <A<：J l)j inai}i ') *画出通丸洙區
1~ -plot ( [vhM^ vtnMlj [0,2001^^ t-)兔画一举垂直孑橫柚的直鉞■用翔坐标儿m列km ■级由从□至|]二加IF - iclab^l <J vfn/s |J, ontsiz^^ □，]n
IS - ylab?l C s(n)" / fcntsi:e J, 16)
2D - tit It C 耒度匚連度的殖系’J f Wit Sire^j lfl),
21 —pauss
於-^i S<[10 is 90 IOO])：夷除制坐垛范匿・细北局部
23- et ext 广V™ar J /f srrtsize' , 18 ■, pt ext <'Smax' / f ontsise ・ 13)
24- hold eft
(1)执行pause之前的代码，得到的结果如下图所示:
该图代表了假设水无限深的情况下，物体速度的变化情况，容易看出速度是有一个上限
的，这也符合直观的感受，因为速度不可能无限大，否则阻力无限大，物体将无法继续运动。
(2 )红线为题给的速度阈值，可见初步判断速度是有可能超过该阈值的，需要进一步
判断沉底的时候是否超过该速度，于是需要做出s和v的关系曲线，即x(2),x(1)的关系曲线, 如下所示:
从图像中看出，s趋于g时候，v也趋于极限值，题给的smax=300*0.3048m, 在编程的
时候，如果要绘制该水平直线，需要将其转化为数组才能进行绘图，否则会出现错误。
(3)但是从这幅图像中，由于s非常大，红线触底，无法直接看出s,v的相交情况，需要进一步限定坐标轴范围进行细化，如下所示:
从图中明显读出，当s达到smax时候，其相应的v>vmax,所以从题给的情况看，工程师
们的说法是更有道理的。
在这个例题中，加深了对微分方程数值解的理解，同时在用matlab绘图时也有了更多
的心得。
绘图时候尤其要注意为数组和数组的对应！另外在通过查询相关知识知道，如果要
绘制直接过（a,b ）和（c,d ）两点的直线，可以直接用plot（[a,c],[b,d]）命令快速实现，方便了一些作图过程。
4-6小船渡河问题
首先给出问题的理论解法：
建立小船运行的数学模型，用极坐标，以B为原点进行计算。
取定BA为极轴。
小船位置设为C,可以得到方程组:
可以解得:
L 1 + sina 鸟 cosa ' 1 一
sina^ 2
接下来用程序进行数值解和理论解答的分析：
(1 )用matlab 运算时，采用普通坐标更加方便，所以以 A 为原点建立直角坐标系,
容易得到
r dX X
=V 1-V 1 - 血
Vx 2 + Crf-y): 据此定义函数m 文件如下所示:
'Ed rtor - C :\U se r^\len pvp\Des kto p\mat I a bteach 已rfillE'myli o me work \h om e v/o-rk2\xi aoc h u an.m
File Edrt Text Go Cell Tools Debug Desktop Window Help
□ ■怡■胃专0劭二：▼此■■已 归*村电
10
function ±f=xi3ochuan(tjx)
d=iOO: vl=1: v2=2：
dx= [vl —v2*x (D/sqr!: (x (1) ' 2 + (d-z (2)「2) ; v2* * d_X * 2 T f sqrt 厲⑴"2~^ (d.-x (2) ^2) J (2)相应的主程序为:
= v^sma — v 2
=卩二CDsa/p
消去t
dp Vising — r a -- = --- ffl£ p v^eosa
积分得到:
111朗i 禽-知
1,1
其中pause前的代码表示进行v1=1 , v2=2时候的数值运算求解绘图。
pause后的代
码表示分别绘制出v仁005,1,1.5,2 时候的运动轨迹并实现题目（2）中要求的数值解与理论解的对比情况，运行的结果如下所示:
① v1=1m/s,v2=2m/s 时候的轨迹图:
'O Figure 1 . @ file Edit View Insert Tools Qesktop Window Help 勺□ 輪|4|鬓只甥®膜尿|□盂l|
小船计算运动轨迹
(2)然后是v1=0,0.5m/s,1m/s,1.5m/s,2m/s 的理论曲线，并实现v1=1m/s 时候的
理论曲线和实际曲线的对比：
该图中有两点需要注意的地方：
1）标注均在曲线的右方，而3=0时候的曲线实际上和y轴重合，所以绘出的曲线没有
显现出来；
2 ）注意到v1=1m/s 时候的曲线，实际上此处还有一条蓝色曲线，图上可以看到蓝色曲线的痕迹，这是v1=1m/s 时候的数值解曲线，可以看出二者的重合度是非常高的，几乎完全重合，这说明之前的数值解的结果是可信的。
另外我们注意到
1 ）当v1=2m/s 的时候，船已经不可能正好到达河的正对岸了，从物理规律也可以直观
的理解，因为当v1 >v2时候，v2的任何一个分量都不可能大于v1，这样一定会随着水流向下游走一定的距离，而不可能完全到达对岸;
2）当水流速度逐渐增大的时候，路径越来越凸向右侧，即要到达正对岸会更加’费力气’，这也是符合生活规律的。
【实验总结】
本次实验总的来说思路较为清晰，通过实验进一步加深了对matlab的了解，巩固了课
堂所学的常微分方程的一些知识。
在matlab的实现过程中，有很多需要注意的细节，需要
耐心细致的去理解和注意，尤其是一些基本命令的实现，编程的时候不能凭感觉，而要符合
规范。
最后就是可以通过直观的生活规律、物理定律对所得的结果进行初步验证，体现了数
值解法的科学性。