Euler_法求微分方程的解

L1 L2 L3 L 4
f ( xk , yk ) f ( xk h / 2, yk hL1 / 2) f ( xk h / 2, yk hL2 / 2) f ( xk h, yk hL3 )
9
四阶 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-')
k k
k 0, 1, 2,, n 1
4
Euler 折线法举例
例：用 Euler 法解初值问题
2x dy y 2 y dx y( 0 ) 1 x [0, 2]
解：取步长 h = (2 - 0)/n = 2/n，得差分方程
x 0 0, y 0 1 2 yk 1 yk h f ( xk , yk ) yk h( yk 2 xk / yk ) x x h k k 1
20
参数说明
[T,Y] = solver(odefun,tspan,y0)
odefun 为显式常微分方程，可以用命令 inline 定义，或 在函数文件中定义，然后通过函数句柄调用。
dy 2 y 2 x 2 2 x 的数值解，求解区间为 [0,0.5] 例 ： 求 dx y( 0 ) 1
2
Euler 折线法
考虑一维经典初值问题
dy f ( x , y) , y( x0 ) y0 , x [a, b] dx
基本思想：用差商代替微商
根据 Talyor 公式，y(x) 在点 xk 处有
y( x ) y( xk ) ( x xk ) y '( xk ) O( x 2 ) y( xk 1 ) y( xk ) hy '( xk ) O ( h2 ) h xk 1 xk
23
求解微分方程小结
Matlab 函数
求解析解（通解或特解），用 dsolve 求数值解（特解），用 ode45、ode23 ...
Matlab 编程
Euler 折线法 Runga-Kutta 方法
24
上机作业
教材 P97：练习 1 ~ 7
25
y( xk 1 ) y( xk ) y( xk 1 ) y( xk ) dy O ( h) dx x h h k
3
初值问题的Euler折线法
具体步骤：分割求解区间，差商代替微商，解代数方程
分割求解区间
xk k 0 为分割点
n
等距剖分： a x0 x1 x2 xn1 xn b 步长：h xk 1 xk ( b a) / n, 差商代替微商
下的特解，并画出解函数的图形。
[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0', ... 'x(0)=1', 'y(0)=0', 't') ezplot(x,y,[0,1.3]);
注：解微分方程组时，如果所给的输出个数与方程个数相同， 则方程组的解按词典顺序输出；如果只给一个输出，则输出 的是一个包含解的结构（structure）类型的数据。
如果需求解的问题是高阶常微分方程，则需将其化为一阶常 微分方程组，此时必须用函数文件来定义该常微分方程组。
d2y 2 dy (1 y ) y 0 例 ： 求解 Ver der Pol 初值问题 dt 2 dt 7, y(0) 1, y '(0) 0 dx1 dt x2 x1 y dx2 (1 x 2 ) x x 令 1 2 1 dt dy x2 dt 7 x (0) 1, x (0) 0, 1 2
y( xk 1 ) y( xk ) dy y( xk 1 ) y( xk ) h y '( xk ) dx x h k k = 0, 1, 2, ..., n-1 y0 y( x0 ) 得方程组： yk 1 yk h f ( xk , yk ) x x h y 是 y (x ) 的近似 k k 1
>> y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') >> ezplot(y);
16
dsolve 举例
dx 5x y e t x |t 0 1 dt 例3：求微分方程组 在初值条件 y |t 0 0 dy x 3y 0 dt
dsolve 的使用
y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 其中 y 为输出， eq1、eq2、...为微分方程，cond1、 cond2、...为初值条件，v 为自变量。
dy x2 例 1：求微分方程 的通解，并验证。 2 xy xe dx
求解器 ODE类型 特点 说明
ode45
非刚性
单步法；4，5 阶 R-K 方法； 大部分场合的首选方法 累计截断误差为 (△x)3
ode23
ode113
非刚性
非刚性
单步法；2，3 阶 R-K 方法； 使用于精度较低的情形 累计截断误差为 (△x)3
多步法；Adams算法；高低精 计算时间比 ode45 短 度均可到 10-3～10-6
solver 为Matlab的ODE求解器（可以是 ode45、ode23、
ode113、ode15s、ode23s、ode23t、ode23tb）
没有一种算法可以有效地解决所有的 ODE 问题，因此 MATLAB 提供了多种ODE求解器，对于不同的ODE， 可以调用不同的求解器。
19
Matlab的ODE求解器
数学实验
实验四
求微分方程的解
1
问题背景和实验目的
自牛顿发明微积分以来，微分方程在描述事物运 动规律上已发挥了重要的作用。实际应用问题通过 数学建模所得到的方程，绝大多数是微分方程。 由于实际应用的需要，人们必须求解微分方程。 然而能够求得解析解的微分方程十分有限，绝大多 数微分方程需要利用数值方法来近似求解。 本实验主要研究如何用 Matlab 来计算微分方程 （组）的数值解，并重点介绍一个求解微分方程的 基本数值解法－－Euler折线法。
是一类求解常微分方程的数值方法 有多种不同的迭代格式
8
Runge-Kutta 方法
用得较多的是 四阶R-K方法（教材第 98 页）
y0 y( x0 ), xk 1 xk h yk 1 yk h (L1 2 L2 2 L3 L4 )/6
其中
dsolve的输出个数只能为一个 或 与方程个数相等
只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。
18
数值求解
[T,Y] = solver(odefun,tspan,y0)
其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (列向量) 中返回的是分割点 的值(自变量)，Y (数组) 中返回的是这些分割点上的近似解， 其列数等于因变量的个数。
22
数值求解举例
先编写函数文件 verderpol.m
function xprime=verderpol(t,x) global mu; xprime=[x(2); mu*(1-x(1)^2)*x(2) - x(1)];
再编写脚本文件 vdpl.m，在命令窗口直接运行该文件
clear; global mu; mu=7; y0=[1; 0]; [t,x]=ode45('verderpol', [0,40], y0); plot(t,x(:,1));
dsolve('Dy=2*x','x'); dsolve('Dy=2*x');
t
％ dy/dx = 2x ％ dy/dt = 2x
若找不到解析解，则返回其积分形式。
15
dsolve 举例
例 2：求微分方程 xy ' y e x 0 在初值条件 y(1) 2e
下的特解，并画出解函数的图形。
6
Euler折线法举例（续）
解析解： y=1/3*(-18-54*x+45*exp(3*x))^(1/3)
解析解
近似解
7
Runge-Kutta 方法
为了减小误差，可采用以下方法：
让步长 h 取得更小一些； 改用具有较高精度的数值方法：
Runge-Kutta (龙格-库塔) 方法
龙格-库塔方法
10
Runge-Kutta 方法
11
Euler 法与 R-K法误差比较
12
Matlab 解初值问题
用 Maltab自带函数 解初值问题 求解析解：dsolve 求数值解：
ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb
13
dsolve 求解析解
17
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)
ode23t 适度刚性 采用梯形算法 适度刚性情形 若 ode15s 多步法；Gear’s 反向数值微分； ode45 失效时，可 刚性
精度中等 尝试使用
ode23s
ode23tb
刚性 刚性
单步法；2 阶Rosebrock 算法； 当精度较低时，计算时 间比 ode15s 短 低精度 梯形算法；低精度 当精度较低时，计算时 间比ode15s短
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]
21
数值求解举例
>> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') >> syms x; diff(y)+2*x*y - x*eபைடு நூலகம்p(-x^2)
14
dsolve 的使用
几点说明
微分方程中用 D 表示对 自变量 的导数，如： Dy y'； D2y y''； D3y y'''
如果省略初值条件，则表示求通解； 如果省略自变量，则默认自变量为
当 h=0.4，即 n=5 时，Matlab 源程序见 fuluA.m
5
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]; end szj plot(szj(:,1),szj(:,2),'or-')