清华大学谢金星数学实验-作业6

实验二常微分方程的数值解

土木系结23 李会平

【实验目的】

1、掌握用matlab解数值微分方程

2、了解龙格-库塔方法的基本原理

3、用这些手段解决一些实际的问题

【实验内容】

4-5 核废料问题

●首先列出问题的运动学方程，由牛顿第二定律，

md2s/dt2=G-F-f

其中m=G/g, f=kv,由于我们熟悉的单位是公制的，所以在定义函数的时候进行了单位转换。以下是函数的定义：

这其实是一个关于s的二阶常微分方程，需要定义两个变量x(1),x(2)将其化为一阶微分方程组，此处x(1),x(2)实际上分别代表速度和深度，相应的微分方程如代码中所示：

dx=[(G-F-k*x(1))/m;x(1)]

（如课堂提醒，中间应该是分号，这点容易出错。）

●接写来进行m文件的命令编写，如下所示：

（1）执行pause之前的代码，得到的结果如下图所示：

该图代表了假设水无限深的情况下，物体速度的变化情况，容易看出速度是有一个上限的，这也符合直观的感受，因为速度不可能无限大，否则阻力无限大，物体将无法继续运动。

（2）红线为题给的速度阈值，可见初步判断速度是有可能超过该阈值的，需要进一步判断沉底的时候是否超过该速度，于是需要做出s和v的关系曲线，即x(2),x(1)的关系曲线，如下所示：

从图像中看出，s趋于∞时候，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，可以得到方程组：

消去t

积分得到：

可以解得：

接下来用程序进行数值解和理论解答的分析：

（1）用matlab运算时，采用普通坐标更加方便，所以以A为原点建立直角坐标系，容易得到

据此定义函数m文件如下所示：

（2）相应的主程序为：

其中pause前的代码表示进行v1=1，v2=2时候的数值运算求解绘图。pause后的代码表示分别绘制出v1=0,0.5,1,1.5,2时候的运动轨迹并实现题目（2）中要求的数值解与理论解的对比情况，运行的结果如下所示：

①v1=1m/s,v2=2m/s时候的轨迹图：

（2）然后是v1=0,0.5m/s,1m/s,1.5m/s,2m/s的理论曲线，并实现v1=1m/s时候的理论曲线和实际曲线的对比：

该图中有两点需要注意的地方：

1）标注均在曲线的右方，而v1=0时候的曲线实际上和y轴重合，所以绘出的曲线没有显现出来；

2）注意到v1=1m/s时候的曲线，实际上此处还有一条蓝色曲线，图上可以看到蓝色曲线的痕迹，这是v1=1m/s时候的数值解曲线，可以看出二者的重合度是非常高的，几乎完全重合，这说明之前的数值解的结果是可信的。

另外我们注意到

1）当v1=2m/s的时候，船已经不可能正好到达河的正对岸了，从物理规律也可以直观的理解，因为当v1≥v2时候，v2的任何一个分量都不可能大于v1，这样一定会随着水流向下游走一定的距离，而不可能完全到达对岸；

2）当水流速度逐渐增大的时候，路径越来越凸向右侧，即要到达正对岸会更加‘费力气’，这也是符合生活规律的。

【实验总结】

本次实验总的来说思路较为清晰，通过实验进一步加深了对matlab的了解，巩固了课堂所学的常微分方程的一些知识。在matlab的实现过程中，有很多需要注意的细节，需要耐心细致的去理解和注意，尤其是一些基本命令的实现，编程的时候不能凭感觉，而要符合规范。最后就是可以通过直观的生活规律、物理定律对所得的结果进行初步验证，体现了数值解法的科学性。