华南理工大学数学实验实验三

dy
（21）
则微分方程为
d2x

y
dy2
=
1+ ( dx )2 dy
，c = a v

x( y0 ) = 0 ，x '( y0 ) = 0
在 MATLAB 中使用 ode45 函数需要将上述方程化为
（22）

x1
=
x2

x2
=
c y
1+
x12
,
c
=
a v
（23）
在函数代码中，待解变量、待解变量其导数、初始值都用数组来表示：
两边积分
������������ √1+������2
=
������
������������ ������
ln|������ + √1 + ������2| = ������ ∙ ln⁡|������| + ������ 则
������ + √1 + ������2 = ������1������������
−y
d2x dy2
=
a
dt dy
又
4
（4） （5） （6） （7） （8） （9） （10） （11） （12） （13） （14） （15） （16）
（17） （18）
dt = dt ds dy ds dy
（19）
而
v = ds dt
（20）
则
y d 2 x = a 1+ ( dx )2
dy2 v
5
plot(x(:,1),y,color(i),'LineWidth',2); %不同速度下的缉私舰 运动轨迹
text(x(180+70*i),y(150+80*i),['\leftarrow k=',num2str(c)]);%固定位置设置箭头
end y1=0:0; %走私船的 y 轴位置 x1=10:5:80;%走私船的 x 轴位置 plot(x1,y1,'m^','linewidth',2);%走私船运动轨迹 xlabel('x');ylabel('y'); hold off; %关闭保持图形功能 end
������
设
������ = ������������
������������
则
������������ ������������
=
������2������ ������������2
分离参数
������ ������������ = ������√1 + ������2
������������
6
代码的完整性、详细清楚的注释和图形的美化都是十分值得学习的地方，今后也 要继续努力。
7
2． 实验任务 问题1
用 dsolve 函数求解下列微分方程
（2）

y(x) = y(x) + 2 y(x) y(0) = 1, y(0) = 0
问题2
我辑私雷达发现，距离 d 处有一走私船正以匀速 a 沿直线行驶，缉私舰 立即以最大速度（匀速 v）追赶。若用雷达进行跟踪，保持船的瞬时速度方 向始终指向走私船，则辑私舰的运动轨迹是怎么的？是否能够追上走私船？
[t,Y]=ode45(odefun,tspan,y0) 说明：
（1） odefun为待解一阶微分方程或方程组的句柄，对应一个M函数文
件，在函数文件中定义微分方程或方程组的结构；
（2） tspan求解区间，y0为初值（初始条件）；
（3） 返回值t为自变量的数据列；
（4） 返回值Y一般是矩阵，每列对应一个待解变量的数据列；
求倒数得
������ − √1 + ������2 = −1 ������−������
������1
以上两个式子相加
������
=
1 2
(������1������������
−
1 ������1
������ −������ )
当 k=1 时，积分得
������
=
1 4
������1������2
15 k=0.6
10
k=0.8 5
0
0
10
20
30
40
50
60
70
80
x
图 2.问题二运行结果
y
4. 实验总结和实验感悟
本次实验的难点在于微分方程的建立，在草稿纸上推演和运算花了大量的时 间，犯了不少错误浪费了不少时间，但也复习了以前学过的微分方程的知识，颇 有收获。另外，学会了 ode45 函数的使用方法，虽然微分方程的句柄定义有点绕 而一开始理解运用起来有点吃力，但随着对例子的熟悉和多次练习使用，还是掌 握了 ode45 函数的使用。最后，在与同学合作的过程中也学到了不少东西，对方
如果能追上，需要多长时间？
1
y M0
M(x, y)
源自文库
d

S0 图 2-14 追缉模型 S
x
3. 实验过程 3.1问题1 3.1.1实验原理 在MATLAB中，dsolve()函数可用来求解微分方程的解析解，格式为： Dsolve(‘equation’,’condition’,’v’) 其中，equation为方程式，condition为条件，v是自变量（缺省为t），一 般的，一阶导数 y 可表示为Dy，二阶导数 y 可表示为D2y，以此类推； 当有多个方程或多个条件时，写多个相应参数即可。 Dsolve()函数求解规则如下 （1）若不带条件，则所得的解中会有积分常数； （2）如果没有显式解，则系统会尝试求隐式解； （3）如果连隐式解都没有，则返回空符号。 3.1.2算法与编程 3.1.2.1算法 直接使用dsolve函数即可。 3.1.2.2编程
实验三 微分方程
地
点：
计算中心 202 房； 实验台号：
03，04
实验日期与时间：
2018 年 4 月 11 日
评 分：
预习检查纪录：
实验教师：
刘小兰
电子文档存放位置：
电子文档文件名： 信息工程 4 班-03，04-陈邦栋，陈冠宏实验三.doc(x)
批改意见：
1． 实验目的
- 了解求微分方程解析解的方法 - 了解求微分方程数值解的方法 - 学习单自由度阻尼系统，观察阻尼系数对系统的影响 - 模拟弹簧振动，学习实时动画编程的原理 - 了解 dsolve,ode45 指令的使用方法 - 了解 Simulink 仿真的设计思想
待解变量 x = x(1); x(2) = x1; x2 ，
其导数 dx = dx(1); dx(2) = x1; x2 ，
初始值 x0 = x( y0 ); x( y0 ) = 0, 0 。
3.2.2.2编程
代码 1：
function dx=odefun1(y,x) global c; %设置全局变量 dx=zeros(2,1);% a column vector dx(1)=x(2); %第一个方程右边 dx(2)=(c*sqrt(1+(x(2))^2))/y;%第二个方程右边
（5） 对方程组，待解变量、其导数、初始值等，全部用数组表示。
3.2.2算法与编程
3.2.2.1算法
手工解： 某一点的切线导数为：
化简的：
������������ = ������0+������������−������
������������
0−������
������
������������ ������������
=
������
−
������������
−
������0
又因为：
������
������2������ ������������2
=
−������
������������ ������������
（1）
（2） （3）
3
������������ = ������������ ∙ ������������ = − 1 ∙ √(������������)2+(������������)2 = − 1 ∙ √1 + (������������)2
−
ln ������ ������1
+
������2
当 k≠1 时，积分得
数值解： 列出方程式
������
=
1 2
( ������1
������+1
������1+������
−
1 ������1(1−������)
������1−������ )
+
������2
dy dy y = − dx x + dx (x0 + at)
代码 3（主代码）：
clear,clc; %清除工作区和命令行窗口 close all; %关闭所有已经打开的图片 vc=[0.2,0.4 0.6 0.8];%全局变量，代表走私船与缉私船的相对速度 odefun2(vc);%传入不同的相对速度
3.2.3实验结果
50
45
40
35
30
25 k=0.2
20 k=0.4
������������ ������������ ������������
������
������������
������
������������
则
������
������2������ ������������2
=
������
∙
√1
+
(������������)2
������������
������ = ������
2
即，结果为 3.2问题2
图 1.问题一运行结果
������
=
2 3
������ −������
+
1 3
������ 2������
3.2.1实验原理
MATLAB对一阶常微分方程的数值解法，一般是基于龙格-库塔法，对
应的函数为ODE，ODE函数有多种类型分别适用于不同类型的常微分方程
的求解，这里使用最常用的ode45函数，调用格式如下：
clear,clc; %清除工作区和命令行窗口 close all; %关闭已经打开的所有图片 y=dsolve('D2y=Dy+2*y','y(0)=1','Dy(0)=0','x'); disp('特解 y=');%在命令行窗口显示固定字符 disp(y);%命令行窗口显示结果
3.1.3实验结果
end
代码 2：
function odefun2(vc) global c;%设置全局变量 color='mrbg';%各条曲线的颜色 hold on; %保持原有的图形 tspan=50:-0.1:0.1;%求解区间 for i=1:length(vc) c=vc(i);%缉私舰和走私船的相对速度 [y,x]=ode45('odefun1',tspan,[0,0]);%求解数值解