MATLAB求解微分方程微分方程求解

然后继 y i 2的 续计 下算 一。 步
此即改进的欧拉法。
12
3、使用泰勒公式
以此方法为基础，有龙格-库塔法、线性多步法等方 法。
4、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时
（k为正整数，h为步长），称它是一个k阶公式。
k越大，则数值公式的精度越高。
•欧拉法是一阶公式，改进的欧拉法是二阶公式。 •龙格-库塔法有二阶公式和四阶公式。 •线性多步法有四阶阿达姆斯外插公式和内插公式。
dy=zeros(2,1);
2 1.5
1 0.5
0 -0.5
-1 -1.5
-2 -2.5
0
500
1000
1500
2000
2500
3Hale Waihona Puke Baidu00
dy(1)=y(2);
dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
2、取t0=0，tf=3000，输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]);
x=simple(x) % 将x简化 y=simple(y) z=simple(z)
结 果 为：x =C3*exp(2*t)+exp(-t)*C1 y =C2*exp(-2*t)+C3*exp(2*t)+exp(-t)*C1 z =C2*exp(-2*t)+C3*exp(2*t)
9
微分方程的数值解
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
0
dy=zeros(3,1);
-0.2
dy(1)=y(2)*y(3);
-0.4
dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);
-0.6
-0.8
-1
0
2
4
6
8
10
12
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]);
实验
Experiments in Mathematics 微分方程求解
1
实验目的 实验内容
1、学会用Matlab求简单微分方程的解析解. 2、学会用Matlab求微分方程的数值解.
1、求简单微分方程的解析解. 2、求微分方程的数值解.
实验软件 MATLAB
2
微分方程的解析解
求微分方程（组）的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
此即欧拉法。
11
2、使用数值积分
对方程y’=f(x,y), 两边由xi到xi+1积分，并利用梯形公式，有：
y ( x i 1 ) y ( x i) x x i i 1 f( t ,y ( t ) d )x i t 1 2 x i[ f( x i,y ( x i) ) f( x i 1 ,y ( x i 1 ))]
（2）编写程序如下:（vdj.m）
[t,y]=ode23('vdpol',[0,20],[3,0]); y1=y(:,1); % 原方程的解 y2=y(:,2); plot(t,y1,t,y2, '--' ) % y1(t),y2(t) 曲线图 pause, plot(y1,y2),grid % 相轨迹图，即y2(y1)曲线
令 y1=x (t), y2 = x’(t)
y1' y2
y2 '
(1
y12 )
y2
y1
y1(0) 3 y2 (0) 0
19
（1）编写M文件 ( 文件名为 vdpol.m):
function dy = vdpol(t,y); dy=zeros(2,1); dy(1)=y(2); dy(2)=(1-y(1)^2)*y(2)-y(1); % 或 dy=[y(2);(1-y(1)^2)*y(2)-y(1)];
无解析表达式！
6
MATLAB软件求解
例4 非线性微分方程
x'(t)2x(t)21 ,x(0)0
x=dsolve('(Dx)^2+x^2=1','x(0)=0')
x = sin(t) -sin(t)
若欲求解的某个数值解，如何求解？ t=pi/2; eval(x)
7
MATLAB软件求解
例5
输入：
dx dt
x
0
x(0) 2; x'(0) 0
解: 令 y1= x，y2= y1’= x’
则微分方程变为一阶微分方程组：
y1' y2 '
y2 1000(1
y12 )
y2
y1
y1(0) 2, y2 (0) 0
1、建立m-文件vdp1000.m如下：
function dy=vdp1000(t,y)
plot(T,Y(:,1),'-')
3、结果如图
17
例 8 解微分方程组.
y1' y2 y3
y2 y3
' '
y1 y3 0.51
y1
y2
y1(0) 0, y2 (0) 1, y3(0) 1
1 0.8
0.6
解 1、建立m-文件rigid.m如下：
0.4
function dy=rigid(t,y) 0.2
返回
13
（三）用Matlab软件求常微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 函数 量值 值
ode45 ode23 ode113 ode15s ode23s
由待解 方程写 成的m-
ts=[t0,tf]， t0、tf为自 变量的初
函数的 初值
文件名 值和终值
plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图 图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线. 18
例9 Van der pol 方程：
x''(t)(1x(t)2)x'(t)x(t)0 x(0)3,x'(0)0
该方程无解析解！
输出：y= tan(t-C1)
（通解）
y1= tan(x+1/4*pi) （特解）
例2 常系数的二阶微分方程
y '' 2 y ' 3 y 0 , y ( 0 ) 1 ,y '( 0 ) 0
输入： y=dsolve('D2y-2*Dy-3*y=0','x') y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')
结果：
y = C1*exp(-x)+C2*exp(3*x) y = 3/4*exp(-x)+1/4*exp(3*x)
5
例3 非常系数的二阶微分方程
x ''( t ) ( 1 x 2 ( t ) x '( t ) ) x ( t ) 0 , x ( 0 ) 3 ,x '( 0 ) 0
x=dsolve('D2x-(1-x^2)*Dx+x=0', ' x(0)=3,Dx(0)=0')
其中参数β=8/3，σ=10，ρ=28 解：1）编写M函数文件(lorenz.m);
2) 数值求解并画三维空间的相平面轨线； (ltest.m)
23
1、 lorenz.m function xdot=lorenz(t,x) xdot=[-8/3,0,x(2);0,-10,10;-x(2),28,-1]*x; 2、ltest.m x0=[0 0 0.1]'; [t,x]=ode45('lorenz',[0,10],x0); plot(t,x(:,1),'-',t,x(:,2),'*',t,x(:,3),'+') pause plot3(x(:,1),x(:,2),x(:,3)),grid on
记号: 在表达微分方程时，用字母 D 表示求微分，D2、D3 等 表示求高阶微分。任何 D 后所跟的字母为因变量，自变量可以
指定或由系统规则选定为缺省。
注意：① y ' Dy，y ' '
D2y
② 自变量名可以省略，默认变量名 ‘t’。
3
MATLAB软件求解
例1
dy1y2, y(0)1
dx
输入：y=dsolve ('Dy=1+y^2') y1=dsolve('Dy=1+y^2','y(0)=1','x')
3x
4y
dy dt
4x
3y
x(0) 0 y(0) 1
[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')
[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1')
输出： x =-exp(3*t)*(C1*cos(4*t)-C2*sin(4*t)) y =exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)) x =exp(3*t)*sin(4*t) y =exp(3*t)*cos(4*t)
xdot = [x2(t)；x3(t)；…；f(t, x1(t), x2(t),…xn(t))]； 2、数值计算（执行以下命令）
[t,x1,x2,…,xn]=ode45(‘fun',[t0,tf], [x1(0),x2(0),…,xn(0)])
16
例7
d 2 x dt2
1000(1
x2
)
dx dt
选择一组状态变量
x 1y ,x 2y & ,L,x ny(n 1 )
x&1 x 2 , x&2 x 3 , L
xn f (t, x1, x2 ,L , xn )
15
注意
x&1 x 2 , x&2 x 3 , L
xn f (t, x1, x2 ,L , xn ) 1、建立M文件函数
function xdot = fun(t,x,y)
14
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，
m-文件中的待解方程组应以x的分量形式写成。
2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组。
y(n) f (t, y, y&,L , y(n1)) y(0), y&(0),L , y(n1)(0)
返回
10
（二）建立数值解法的一些途径
设xi1xi h, i0,1,2,n1, 可用以下离解 散微 化分 方方 y'f(x,)y y(x0)y0
1、用差商代替导数 若步长h较小，则有
y'(x)y(xh)y(x) h
故有公式：
y yi01yy (x i0 )h(fxi,yi) i0,1,n 2-1,
20
计算结果
蓝色曲线
3
——y（1）； 2
（原方程解）
1
y(1)and y(2)
0
红色曲线 -1 ——y（2）； -2
-3 0
Van Der Pol Solution
5
10
15
Time,Second
20
21
y2
相相相相 3
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
y1
22
例10 考虑Lorenz模型： x 1(t)x1(t)x2(t)x3(t) x2(t)x2(t)x3(t) x3(t)x2(t)x1(t)x2(t)x3(t)
（一）常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多
得不出一般解。而在实际上对初值问题，一般是要求得 到解在若干个点上满足规定精确度的近似值，或者得到 一个满足精确度要求的便于计算的表达式。
因此，研究常微分方程的数值解法是十分必要的。
对常微分 ：方 yy('x0程 f)(x,y0y)， 其数值解是指x0由 开初 始始 的若干离 x值 散处 的，x即 0 x对 1x2 xn， 求出准y确 (x1)值 , y(x2),,y(xn)的相应近y1,似 y2, 值 ,yn。
24
计算结果如下图：
50
40 30
30
20
20
10
0 10
-10
0
-20
20
-10
60
0
40
-20
0
2
4
6
8
10
20 -20 0
图中，x1的图形为实线(蓝），x2的图形为“*” 线(绿)， x3的图形为“+”线（红）。取[t0，tf]=[0，10]。
故有公式：
yi1yi h2[f(xi,yi)f(xi1,yi1)] y0 y(x0)
实际应用时，与欧拉公式结合使用：
yi( 01 ) yi h(fxi,yi) yi( k1 1) yi h 2[f(xi,yi)f(xi1,yi( k1 ))]k0,1,2,
对于已， 当 给满 y 的 i( k 1 1 ) 足 y 精 i( k 1 ) 确 时 取 y 度 ， i 1y （ i k 1 1 ） ，
8
例 6 求微分方程组的通解.
dx
dt dy
dt dz
dt
2x 4x 4x
3y 5y 4y
3z 3z 2z
解 输入命令 ：
[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');