exp04w_微分方程的解

2012-8-6
13
Euler折线法举例（续）
解析解： y
2 5 3x e 2x 3 3
1/3
解析解
近似解
2012-8-6
14
Runge-Kutta 方法
为了减小误差，可采用以下方法：
让步长 h 取得更小一些； 改用具有较高精度的数值方法：
Runge-Kutta (龙格-库塔) 方法
16
四阶 R-K 方法源程序
clear; f=sym('y+2*x/y^2'); a=0; b=2; h=0.4; n=(b-a)/h+1; % n=(b-a)/h; x=0; y=1; szj=[x,y]; for i=1:n-1 % i=1:n l1=subs(f,{'x','y'},{x,y}); l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2}); l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2}); l4=subs(f,{'x','y'},{x+h,y+l3*h}); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=[szj;x,y]; end plot(szj(:,1),szj(:,2), 'dg-')
ode23t 适度刚性 采用梯形算法 适度刚性情形 多步法；Gear’s 反向数值微 若 ode45 失效时，可 ode15s 刚性
分；精度中等 尝试使用
ode23s
ode23tb
2012-8-6
刚性 刚性
单步法；2 阶Rosebrock 算 当精度较低时，计算时 法；低精度 间比 ode15s 短 梯形算法；低精度 当精度较低时，计算时 间比ode15s短
2012-8-6
6
dsolve 举例
例： [x,y]=dsolve('Dx+5*x=0','Dy-3*y=0', ...
'x(0)=1', 'y(0)=1','t')
r = dsolve('Dx+5*x=0','Dy-3*y=0', ... 'x(0)=1', 'y(0)=1','t') 这里返回的 r 是一个 结构类型 的数据 r.x r.y %查看解函数 x(t) %查看解函数 y(t)
2012-8-6 22
数值求解举例
如果需求解的问题是高阶常微分方程，则需将其化为一阶常 微分方程组，此时需用函数文件来定义该常微分方程组。
d y dy 2 (1 y ) y 0 2 dt dt y ( 0 ) 1 , y '( 0 ) 0 , 7
2012-8-6 17
Runge-Kutta 方法
2012-8-6
18
Euler 法与 R-K法误差比较
2012-8-6
19
Matlab函数数值求解
[T,Y] = solver(odefun,tspan,y0)
其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (向量) 中返回的是分割点的 值(自变量)，Y (向量) 中返回的是解函数在这些分割点上的函 数值。
Matlab基础
第六讲
求微分方程的解
2012-8-6
1
问题背景和课程要求
自牛顿发明微积分以来，微分方程在描述事物运 动规律上已发挥了重要的作用。实际应用问题通过 数学建模所得到的方程，绝大多数是微分方程。 由于实际应用的需要，人们必须求解微分方程。 然而能够求得解析解的微分方程十分有限，绝大多 数微分方程需要利用数值方法来近似求解。 本节主要介绍求解微分方程（组）的 Matlab 命令. 研究如何用 Matlab 来计算微分方程（组）的数值解， 并重点介绍一个求解微分方程的基本数值解法－－ Euler折线法。 2012-8-6 2
Euler 折线法源程序
clear f=sym(‘y+2*x/y^2'); a=0; b=2; h=0.4; n=(b-a)/h+1; % n=(b-a)/h; x=0; y=1; szj=[x,y]; for i=1:n-1 % i=1:n y=y+h*subs(f,{'x','y'},{x,y}); x=x+h; szj=[szj;x,y];%每循环一次矩阵szj增加一行 end szj plot(szj(:,1),szj(:,2),'or-')
其中
L1 L2 L3 L 4
2012-8-6
f ( x k , yk ) f ( x k h / 2 , y k h L1 / 2 ) f ( x k h / 2 , yk h L2 / 2 ) f ( x k h , y k h L3 )
10
dy dx x k
2012-8-6

y ( x k 1 ) y ( x k ) h
O (h)
初值问题的Euler折线法
具体步骤：分割求解区间，差商代替微商，解代数方程
分割求解区间
x k k 0 为分割点
n
等距剖分： a x 0 x 1 x 2 x n 1 x n b 步长：h x k 1 x k ( b a ) / n , 差商代替微商
Matlab提供的ODE求解器
求解器 ODE类型 特点 说明
ode45
非刚性
单步法；4，5 阶 R-K 方法；大部分场合的首选方法 累计截断误差为 (△x)3
ode23
ode113
非刚性
非刚性
单步法；2，3 阶 R-K 方法；使用于精度较低的情形 累计截断误差为 (△x)3
多步法；Adams算法；高低精 计算时间比 ode45 短 度均可到 10-3～10-6
x 0 0 , y0 1 y k 1 y k h f ( x k , y k ) y k h ( y k 2 x k / y k ) x k 1 x k h
当 h=0.4，即 n=5 时，Matlab 源程序见 fulu1.m
2012-8-6 12
solver 为Matlab的ODE求解器（可以是 ode45、ode23、
ode113、ode15s、ode23s、ode23t、ode23tb）
没有一种算法可以有效地解决所有的 ODE 问题，因此 MATLAB 提供了多种ODE求解器，对于不同的ODE， 可以调用不同的求解器。
2012-8-6 20
y ( x ) y ( x k ) ( x x k ) y '( x k ) O ( x )
2
y ( x k 1 ) y ( x k ) h y '( x k ) O ( h )
2
h x k 1 x k
y ( x k 1 ) y ( x k ) h
求解微分方程(组)的Matlab 命令
用 Maltab自带函数 求解 求解析解：dsolve 可求通解,也可求特解 求数值解：
ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb
2012-8-6
3
dsolve 求解析解
dsolve 的使用
围为 [0,0.5]
fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1);
注：也可以在 tspan 中指定对求解区间的分割，如：
[x,y]=ode23(fun,[0:0.1:0.5],1); %此时 x=[0:0.1:0.5]
y(x 1 ), y ( x 2 ), , y ( x n ) 的相应近似值
2012-8-6
9
Euler 折线法
考虑一维经典初值问题
dy dx f ( x , y ) , y ( x 0 ) y0 , x [ a , b ]
基本思想：用差商代替微商
根据 Talyor 公式，y(x) 在点 xk 处有
2012-8-6
k = 0, 1, 2, ..., n-1 yk 是 y (xk) 的近似
11
Euler 折线法举例
例：用 Euler 法解初值问题
2x dy y 2 dx y y(0) 1 x [0, 2]
解：取步长 h = (2 - 0)/n = 2/n，得差分方程
y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 其中 y 为输出， eq1、eq2、...为微分方程，cond1、 cond2、...为初值条件，v 为自变量。
例 1：求微分方程
dy dx
2 xy xe
x
2
的通解，并验证。
>> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') >> syms x; diff(y)+2*x*y - x*exp(-x^2)
t
％ dy/dx = 2x ％ dy/dt = 2x