几种数值积分方法

如果省略初值条件，则表示求通解； 如果省略初值条件，则表示求通解； 如果省略自变量， 如果省略自变量，则默认自变量为
dsolve('Dy=2*x','x'); dsolve('Dy=2*x');
t
％ dy/dx = 2x ％ dy/dt = 2x
若找不到解析解，则返回其积分形式。 若找不到解析解，则返回其积分形式。
Runge-Kutta (龙格 库塔 方法 龙格-库塔 龙格 库塔)
龙格-库塔方法 龙格 库塔方法
是一类求解常微分方程的数值方法 一类求解常微分方程的数值方法 有多种不同的迭代格式
Runge-Kutta 方法
四阶R-K方法 用得较多的是 四阶 方法
y0 = y ( x0 ), x k +1 = x k + h yk +1 = yk + h (L1 + 2 L2 + 2 L3 + L4 )/6
Euler折线法举例（续） 折线法举例（ 折线法举例
解析解： 解析解： y=1/3*(-18-54*x+45*exp(3*x))^(1/3)
解析解
近似解
Runge-Kutta 方法
为了减小误差，可采用以下方法： 为了减小误差，可采用以下方法：
取得更小一些； 让步长 h 取得更小一些； 改用具有较高精度的数值方法： 改用具有较高精度的数值方法：
Matlab提供的 提供的ODE求解器 提供的 求解器
求解器 ODE类型 类型 特点 说明
ode45 ode23
ode113
非刚性 非刚性 非刚性
单步法； ， 方法；大部分场合的首选方法 单步法；4，5 阶 R-K 方法；大部分场合的首选方法 累计截断误差为 (△x)3 △ 单步法； ， 方法； 单步法；2，3 阶 R-K 方法；使用于精度较低的情形 累计截断误差为 (△x)3 △ 多步法；Adams算法；高低精 计算时间比 ode45 短 算法； 多步法； 算法 度均可到 10-3～10-6
dsolve 举例
例 2：求微分方程 xy ' + y − e x = 0 在初值条件 y (1 ) = 2 e ：
下的特解，并画出解函数的图形。 下的特解，并画出解函数的图形。
>> y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') >> ezplot(y);
dsolve的输出个数只能为一个 或 与方程个数相等。 的输出个数只能为一个 与方程个数相等。 只有很少一部分微分方程（ 只有很少一部分微分方程（组）能求出解析解。 能求出解析解。 大部分微分方程（ 只能利用数值方法求数值解。 数值方法求数值解 大部分微分方程（组）只能利用数值方法求数值解。
Dsolve举例 马尔萨斯人口模型 举例--马尔萨斯人口模型 举例
求微分方程的解
研究具有实际应用背景的数学函数 物理量在空间的分布情况其及随时间变化的规律。 物理量在空间的分布情况其及随时间变化的规律。
物理现象
物理定律
微分方程
数值解
建模型: 物理学家、化学家、生物学家、 建模型 物理学家、化学家、生物学家、社会学家
问题背景和实验目的
自牛顿发明微积分以来， 自牛顿发明微积分以来，微分方程在描述事物运 动规律上已发挥了重要的作用。实际应用问题通过 动规律上已发挥了重要的作用。 数学建模所得到的方程，绝大多数是微分方程。 数学建模所得到的方程，绝大多数是微分方程。 由于实际应用的需要，人们必须求解微分方程。 由于实际应用的需要，人们必须求解微分方程。 然而能够求得解析解的微分方程十分有限， 然而能够求得解析解的微分方程十分有限，绝大多 数微分方程需要利用数值方法来近似求解。 数微分方程需要利用数值方法来近似求解。 首先重点介绍一个求解微分方程的基本数值解法 －－Euler折线法。 折线法。 －－ 折线法
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) 查看解函数
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-')
单摆问题 天体力学二体问题 电磁场计算············ 电磁场计算
y(t) = exp( at + b)
人口数据拟合实验
Matlab函数数值求解 函数数值求解
[T,Y] = solver(odefun,tspan,y0)
为初值条件， 为求解区间； 其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 为求解区间 在数值求解 自动对求解区间进行分割， 向量 向量) 时自动对求解区间进行分割，T (向量 中返回的是分割点的 自变量)， 向量 向量) 值(自变量 ，Y (向量 中返回的是解函数在这些分割点上的函 自变量 数值。 数值。
某时刻人口数为y0, 在以后某一时间 t 某时刻人口数为 人口数y 的变化率正比于y 人口数 的变化率正比于 y′ = ry t > 0 y(t)=y0 e r t Thomas Malthus (1766-1834) y ( 0) = y0 弹性力学问题 MATLAB符号求解方法 符号求解方法 y=dsolve('Dy=r*y','t') y =C1*exp(r*t) 令 r = a, y0 = exp( b ) a = ?, b = ?
dsolve 求解析解
dsolve 的使用
y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 为输出， 为微分方程， 其中 y 为输出， eq1、eq2、...为微分方程，cond1、 为微分方程 cond2、...为初值条件，v 为自变量。 为初值条件， 为自变量。
dsolve 举例
dx + 5x + y = et dt x |t = 0 = 1 例3：求微分方程组 ： 在初值条件 y |t = 0 = 0 dy − x − 3y = 0 dt
下的特解，并画出解函数的图形。 下的特解，并画出解函数的图形。
[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]);
2 例 1：求微分方程 dy + 2 xy = xe − x 的通解，并验证。 ： 的通解，并验证。 dx
>> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') >> syms x; diff(y)+2*x*y - x*exp(-x^2)
dsolve 的使用
几点说明
的导数， 微分方程中用 D 表示对 自变量 的导数，如： Dy y'； D2y ； y''； D3y ； y'''
Euler 折线法
考虑一维经典初值问题 考虑一维经典初值问题
dy = f ( x , y ) , y ( x 0 ) = y0 , x ∈ [ a , b ] dx
基本思想：用差商代替微商 基本思想：
公式， 根据 Talyor 公式，y(x) 在点 xk 处有
y ( x ) = y ( x k ) + ( x − x k ) y '( x k ) + O ( ∆ x 2 ) y ( x k +1 ) = y ( x k ) + hy '( x k ) + O ( h 2 ) h = x k +1 − x k
注：解微分方程组时，如果所给的输出个数与方程个数相同， 解微分方程组时，如果所给的输出个数与方程个数相同， 按词典顺序输出 则方程组的解按词典顺序输出；如果只给一个输出， 则方程组的解按词典顺序输出；如果只给一个输出，则输出 的是一个包含解的结构 结构（ 的是一个包含解的结构（structure）类型的数据。 ）类型的数据。
solver 为Matlab的ODE求解器（可以是 ode45、ode23、 求解器（ 的 求解器
ode113、ode15s、ode23s、ode23t、ode23tb）
问题， 没有一种算法可以有效地解决所有的 ODE 问题，因此 MATLAB 提供了多种ODE求解器，对于不同的ODE， 提供了多种 求解器，对于不同的 ， 求解器 可以调用不同的求解器 求解器。 可以调用不同的求解器。
Runge-Kutta 方法
Euler 法与 R-K法误差比较 法误差比较
Matlab 解初值问题
用 Maltab自带函数 解初值问题 自带函数 求解析解： 求解析解：dsolve 求数值解： 求数值解：
ode45、ode23、 、 、 ode113、ode23t、ode15s、 、 、 、 ode23s、ode23tb 、
k k
Βιβλιοθήκη Baidu
k = 0 , 1, 2 , K , n − 1
Euler 折线法举例
例：用 Euler 法解初值问题 2x dy = y+ 2 y x ∈ [0, 2 ] dx y(0) = 1 解：取步长 h = (2 - 0)/n = 2/n，得差分方程 ，
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 = x + h k k +1
y ( x k +1 ) − y ( x k ) y ( x k +1 ) − y ( x k ) dy = + O ( h) ≈ dx x h h k
初值问题的Euler折线法 折线法 初值问题的
具体步骤：分割求解区间，差商代替微商， 具体步骤：分割求解区间，差商代替微商，解代数方程
分割求解区间
{ x k }k = 0 为分割点
n
等距剖分： 等距剖分： a = x 0 < x1 < x 2 < L < x n −1 < x n = b 步长：h = x k +1 − x k = ( b − a ) / n , 步长： 差商代替微商
y ( x k +1 ) − y ( x k ) dy y ( x k +1 ) ≈ y ( x k ) + h y '( x k ) ≈ dx x h k k = 0, 1, 2, ..., n-1 y0 = y ( x 0 ) 得方程组： 得方程组： y k +1 = y k + h f ( x k , y k ) x = x + h y 是 y (x ) 的近似 k k +1
其中
L1 = L2 = L3 = L4 =
f ( x k , yk ) f ( x k + h / 2 , yk + hL1 / 2 ) f ( x k + h / 2 , yk + hL2 / 2 ) f ( x k + h, yk + hL3 )
方法源程序 四阶 R-K 方法源程序
当 h=0.4，即 n=5 时，Matlab 源程序见 fulu1.m ，
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-')