《常微分方程》综合性实验实验报告实验班级05应数（3）学生姓名江晓荣学生学号200530770314指导老师方平华南农业大学理学院应用数学系实验 微分方程在数学建模中的应用及数值解的求法一、实验目的1．了解常微分方程的基本概念。

2．常微分方程的解了解析解和数值解。

3．学习、掌握MA TLAB 软件有关求解常微分方程的解析解和数值解的有关命令。

4. 掌握微分方程在数学建模中的应用。

二、实验内容1．用MA TLAB 函数dsolve 符号求解常微分方程的通解和特解。

2．用MA TLAB 软件数值求解常微分方程。

三、实验准备1．用MA TLAB 求常微分方程的解析解的命令用MA TLAB 函数dsolve 求常微分方程()(,,,,,,)0n F x y y y y y ''''''= （7.1）的通解的主要调用格式如下：S=dsolve('eqn', 'var')其中输入的量eqn 是改用符号方程表示的常微分方程(,,,2,)0F x y Dy D y Dny = ，导数用D 表示，2阶导数用D2表示，以此类推。

var 表示自变量，默认的自变量为t 。

输出量S 是常微分方程的解析通解。

如果给定常微分方程（7.1）的初始条件()00010(),(),,()n n y x a y x a y x a '=== ，则求方程（7.1）的特解的主要调用格式如下：S=dsolve('eqn', 'condition1 ',…'conditonn ','var')其中输入量eqn ，var 的含义如上， condition1,…conditonn 是初始条件。

输出量S 是常微分方程的特解。

2．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。

考虑一阶常微分方程初值问题000()(,()),()f y t f t y t t t t y t y '=<<⎧⎪⎨=⎪⎩ 其中12(,,,)T m y y y y = ，12(,,,)T m f f f f = ，010200(,,,)T m y y y y = 。

所谓数值解法，就是寻求()y t 在一系列离散节点01n f t t t t <<<≤ 上的近似值,0,1,,k y k n = ，称1k k k h t t +=-为步长，通常取为常量h 。

最简单的数值解法是Euler 法。

Euler 法的思路很简单：在节点处用差商近似代替导数1()()()k k k y t y t y t h+-'≈这样导出计算公式 1(,),0,1,2,k k k k y y hf t y k +=+=它能求解各种形式的微分方程。

Euler 法也称折线法。

Euler 法只有一阶精度，改进方法有二阶Runge-Kutta 法、四阶Runge-kutta 法、五阶Runge-kutta-Felhberg 法和线性多步法等，这些方法可用于解高阶常微分方程（组）初值问题。

边值问题采用不同方法，如差分法、有限元法等。

数值解法的主要据点是它缺乏物理意义。

3．用MA TLAB 数值求解微分方程的命令MATLAB 中主要用函数ode45，ode23，ode15s 求常微分方程的数值解。

其中ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。

ode23与ode45类似，只是精度低一些。

ode15s 用来求解刚性方程组。

这三个函数的调用格式基本相同，下面介绍ode45的调用格式。

[tout,yout]=ode45(‘yprime ’,[t0,tf],y0)此命令是采用变步长四阶Runge-kutta 法和五阶Runge-kutta-Felhberg 法求数值解，yprime 是表示f(t,y)的M 文件名，t0表示自变量的初始值，tf 表示自变量的终值，y0表示初始向量值。

输出向量tout 表示节点0,1(,,)T n t t t ，输出矩阵yout 表示数值解，每一列对应y 的一个分量。

若无输出参数，则自动作出图形。

四、实验方法与步骤1求下列常微分方程的通解(1)dx at dt =- 22(2)(s i n c o s )d y d y a b x x d x d x +=+ 2(3)(1)dy y y dx =- （4）22sin d y y dx= 求通解相应的MATLAB 代码为>> x1=dsolve('Dx=-a*t')>> y2=dsolve('D2y+a*Dy=b*(sin(x)+cos(x))','x')>> y3=dsolve('Dy=y^2*(1-y)','x')>> y4=dsolve('D2y=sin(y)','x')运行后屏幕显示常微分方程（1）、（2）、（3）、（4）的通解x1、y1、y2、y3、y4依次为： x1 =-1/2*a*t^2+C1y2 =b*(-a*cos(x)+a*sin(x)-sin(x)-cos(x))/(a^2+1)+C1+C2*exp(-a*x)y3= x+1/y-log(y)+log(-1+y)+C1=0y4=[ Int(1/(-2*cos(a)+C1)^(1/2),a=``..y)-x-C2=0,-Int(1/(-2*cos(a)+C1)^(1/2),a=``..y)-x-C2=0]2求下列常微分方程在给定初始条件下的特解。

22(1)1,(0)0dy y y dx ⎛⎫+== ⎪⎝⎭2220()()()(2),(0)1,x ad f x d f x d f x a f dx dx dx π===-=求特解相应的MATLAB 代码为>> y=dsolve('(Dy)^2+y^2=1','y(0)=0','x'),>> f=dsolve('D2f=-a^2*(Df)','f(0)=1,Df(pi/a)=0','x'),运行后屏幕显示常微分方程（1）、（2）在给定初始条件下的特解依次如下：y =[ sin(x)][ -sin(x)]f =13求解微分方程1,(0)1dy y t y dt=-++=的解析解和数值解，并进行比较。

相应的MA TLAB 代码为>> y=dsolve('Dy=-y+t+1','y(0)=1'),运行结果即方程的解析解为。

y =t+exp(-t)下面再求其数值解，先编写M 文件fun1.m 。

function f=fun1(t,y)f=-y+t+1；再运行相应的MATLAB 代码>> clear;close;t=0:0.1:1;>> y=t+exp(-t);plot(t,y); %画解析解的图形>> hold on; %保留已经画好的图形，如果下面再画图，两个图形合并在一起 >> [t,y]=ode45('fun1',[0,1],1);>> plot(t,y,'ro'); %画数值解图形，用红色小圈画>> xlabel('t'),ylabel('y');运行结果画图分析解析解和数值解吻合情况4.综合实验（数学建模）：求解如下问题之一：1）．小船从河边点o 处出发驶向以对岸（两岸为平行直线），设船速为a ，船行方向始终与河岸垂直，又设河宽为h ，河中任一点处的水流速度与该点到两岸距离的乘积成正比（比例系数为k ），求小船的航行路线及船到达对岸的位置。

2）．（目标的跟踪问题）设位于坐标原点的甲舰向位于x 轴上点)0,1(A 处的乙舰发射制导导弹，导弹头始终对准乙舰．如果乙舰以最大的速度0v （0v 是常数）沿平行于y 轴的直线行驶，导弹的速度是50v ，求导弹运行的曲线方程及乙舰行驶多远时，它将被导弹击中？ （说明：以上两题选做其中之一，或者自选其它题目）一． 问题：小船从河边点o 处出发驶向以对岸（两岸为平行直线），设船速为a ，船行方向始终与河岸垂直，又设河宽为h ，河中任一点处的水流速度与该点到两岸距离的乘积成正比（比例系数为k ），求小船的航行路线及船到达对岸的位置。

二． 实验目的：1．常微分方程的解了解析解和数值解。

2．学习、掌握MA TLAB 软件有关求解常微分方程的解析解和数值解的有关命令。

三． 实验内容及要求：以船的起点为原点，水流方向为y 轴正向，船的航行方向为x 方向建立坐标系。

设在t 时刻船在x 轴方向上运行了x ，在y 轴方向上运行了y ，在很小的时间dt 内，可以认为船在y 轴方向作均速v 运动，移动了dy=vdt;在x 轴方向运动了dx=adt;其中v=k*x*(h-x),由此可建立微分方程组：**()(0)dy k x h x h dt t dx a a dt⎧=-⎪⎪≤≤⎨⎪=⎪⎩ 当t=0时，x 和y 都为0，将方程组中的t 消去得下微分方程：**()dy k x h x dx a=- ()0t h ≤≤ ()00y =输入相应的MA TLAB 代码：>> y=dsolve('Dy=(k/a)*x*(h-x)','y(0)=0','x')运行结果即方程的解析解为y =-1/3*k/a*x^3+1/2*k/a*h*x^2这也是小船的运动轨迹方程。

当x=h 时，船已运行到对岸，此时y 的值为36kh a ,则船在坐标中的位置是（h, 36kh a）. 下面就110.002**k s m --=, 2/a m s =,h=50m.画图看一下小船的运行曲线。

此时方程化为：3211**300040y x x =-+ (050x ≤≤) 下面画解析式的图其数值解图，先编写M 文件function f=fun1(x,y)f=(-1/3000)*x^3+1/40*x^2;function f=fun2(x,y)f=1/1000*x*(50-x);再运动相应的MATLAB 代码>> clear;close;x=0:1:50;fplot('fun1',[0,50]); %画解析解的图形hold on; %保留已经画好的图形，如果下面再画图，两个图形合并在一起[x,y]=ode45('fun2',[0,50],0);plot(x,y,'g*'); %画数值解图形，用绿色星号xlabel('x'),ylabel('y');运行结果画图四．结果分析：两曲线比较好地吻合在一起，此也为小船的运行曲线，由图可知小船的坐标为（50,125/6），也即是小船运动到对岸时已顺不流方向运动了20.8米。