matlab解微分方程

f (t , y ( t )) dt ≈
x i +1 − x i [ f ( x i , y ( x i )) + f ( x i +1 , y ( x i +1 ))] 2
3、使用泰勒公式
4、数值公式的精度
k越大，则数值公式的精度越高.
• 龙格-库塔法有二阶公式和四阶公式. • 线性多步法有四阶阿达姆斯外插公式和内插公式.
N
U
ST -M
d 2 y dy 2 + 4 + 29 y = 0 dx dx y (0) = 0, y ' (0) = 15
at la
b
例3
[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z', 't') ; x=simple(x) % 将x化简 y=simple(y) z=simple(z)
此即改进的欧拉法.
N
然后继续下一步 y i + 2 的计算。
U
(k ) （ k +1） k +1) 对于已给的精确度 ε ， 当满足 y i(+ 时， 取 y − < = y ε y 1 i +1 i +1 i +1 ，
ST -M
at la
b
y ( x i +1 ) − y ( x i ) =
∫
xi + 1
又根据题意，弧 OP 的长度为 AQ 的 5 倍， 即
∫
x
0
1 + y ' 2 dx = 5v0 t
U
(1)
ST -M
(2)
at la
设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰 发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度 v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求 导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？ 解法一：（解析法）
1、微分方程的解析解
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
解
dt 输入命令： u=dsolve('Du=1+u^2','t')
结果：u = tan (t+c1)
N
记号: 在表达微分方程时，用字母 D 表示求微分，D2、 D3 等表示求高阶微分.任何 D 后所跟的字母为因变量，自 变量可以指定或由系统规则选定为确省. d2y 例如，微分方程 = 0 应表达为：D2y=0. 2 dx du
b
2 4 6 8 10 12
3、结果如图： 图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.
N
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
b
解即为导弹的运行轨迹 :
4 5
N
5 5 ) 处时被导弹击中. 当 x = 1时 y = ，即当乙舰航行到点 (1, 24 24 y 5 被击中时间 为 : t = = . 若 v0=1, 则在 t=0.21 处被击中. v0 24v0
U
ST -M
5 5 5 y = − (1 − x) + (1 − x) + 8 12 24
例1 求
= 1 + u 2 的通解 .
U
ST -M
at la
b
求微分方程（组）的解析解命令 :
例如求下例微分方程的特解。
N
如想画出函数在自变量x在区间 [-10,10]的函数图像，可在 命令窗口中输入： ezplot(y,[-10,10])
U
ST -M
在Matlab命令窗口中输入： y=dsolve('Dy=exp(x)','y(0)=exp(1)','x') 输出结果为： y =exp(x)-1+exp(1)
at la
dy = e x , y (0) = e dx
bLeabharlann 例2求微分方程的特解.解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 结 果 为 : y =3*exp(-2*x)*sin(5*x) 作图命令：ezplot(y,[1.0,4])
U
ST -M
3、数学建模实例
（一）导弹追踪问题
假设导弹在 t 时刻的位置为 P(x(t), y(t))，乙舰位于 Q (1, v 0 t ) .
由于导弹头始终对准乙舰，故此时直线 PQ 就是导弹的轨迹曲线弧 OP 在点 P 处的切线，
N
v0t − y 即有 y' = 1− x 即 v0 t = (1 − x) y '+ y
N
U
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差 .
ST -M
ode45 ode23 ode113 ode15s ode23s
由待解 方程写 成的m文件名
运行结果为： x = (c1-c2+c3+c2e -3t-c3e-3t)e2t , y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t , z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t .
N
U
ST -M
解 输入命令 ：
at la
求微分方程组的通解 . dx dt = 2 x − 3 y + 3 z dy = 4 x − 5 y + 3z dt dz = 4 x − 4 y + 2 z dt
N
U
1、建立m-文件vdp1000.m如下： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
y1 ' = y2 , 2 y ' = 1000(1 − y 2 1 ) y 2 − y1 , y (0) = 2, y (0) = 0. 1 2
at la
dx 2 dy 2 2 1．设导弹速 度恒为 w ，则 ( ) + ( ) = w dt dt
b
（1）
ST -M
（2 ）
3．因乙舰以速度v0沿直线x=1运动，设v0=1，则w=5，X=1，Y=t
at la
-0.4 -0.6 -0.8 -1 0
y1 ' = y2 y3 , y 2 ' = − y1 y3 , y3 ' = −0.51 y1 y 2 , y1 (0) = 0, y2 (0) = 1, y3 (0) = 1.
1 0.8 0.6 0.4 0.2 0 -0.2
at la
6 5
b
1 (1 − x) y " = 1 + y '2 , (3) 5 初值条件为: y (0) = 0 ， y ' (0) = 0 .
由(1),(2)消去t整理得模型:
解法二：(数值解) 令y1=y,y2=y1’，将方程（3）化为一阶微分方程组. 1 (1 − x) y '' = 1 + y '2 5
dx X − x dt = λ 即： ， λ >0 dy Y − y dt w dx = ( X − x) dt 2 2 ( X − x) + (Y − y ) 消去λ得： dy w = (Y − y ) 2 2 dt ( X − x) + (Y − y )
N
U
ST -M
at la
b
例4
解: 令 y1=x，y2=y1’
d 2x 2 dx 2 − 1000(1 − x ) − x = 0, dt dt x (0) = 2; x '(0) = 0.
则微分方程变为一阶微分方程组：
ST -M
2、取t0=0，tf=3000，在matlab主命令窗口中输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]); plot(T,Y(:,1),'-') 3、结果如图：
x0=0; xf=0.9999; [x,y]=ode15s('eq1',[x0 xf],[0 0]); plot(x,y(:,1),'b.') hold on; y=0:0.01:2; plot(1,y, 'r*')
结论: 导弹大致在（1，0.2）处击中乙舰.
N
U
ST -M
1. 建立m-文件eq1.m function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(2)^2)/(1-x); 2. 取x0=0; xf=0.9999，建立主程序ff6.m如下：
故有公式：
yi +1 = yi + hf ( xi , yi ), i = 0,1,2, L ,n-1. y 0 = y ( x0 ),