当前位置:文档之家› 几种数值积分方法

几种数值积分方法


Runge-Kutta 方法
Euler 法与 R-K法误差比较 法误差比较
Matlaห้องสมุดไป่ตู้ 解初值问题
用 Maltab自带函数 解初值问题 自带函数 求解析解: 求解析解:dsolve 求数值解: 求数值解:
ode45、ode23、 、 、 ode113、ode23t、ode15s、 、 、 、 ode23s、ode23tb 、
y ( x k +1 ) − y ( x k ) y ( x k +1 ) − y ( x k ) dy = + O ( h) ≈ dx x h h k
初值问题的Euler折线法 折线法 初值问题的
具体步骤:分割求解区间,差商代替微商, 具体步骤:分割求解区间,差商代替微商,解代数方程
分割求解区间
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);
注:解微分方程组时,如果所给的输出个数与方程个数相同, 解微分方程组时,如果所给的输出个数与方程个数相同, 按词典顺序输出 则方程组的解按词典顺序输出;如果只给一个输出, 则方程组的解按词典顺序输出;如果只给一个输出,则输出 的是一个包含解的结构 结构( 的是一个包含解的结构(structure)类型的数据。 )类型的数据。
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-')
如果省略初值条件,则表示求通解; 如果省略初值条件,则表示求通解; 如果省略自变量, 如果省略自变量,则默认自变量为
dsolve('Dy=2*x','x'); dsolve('Dy=2*x');
t
% dy/dx = 2x % dy/dt = 2x
若找不到解析解,则返回其积分形式。 若找不到解析解,则返回其积分形式。
某时刻人口数为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 为自变量。 为初值条件, 为自变量。
solver 为Matlab的ODE求解器(可以是 ode45、ode23、 求解器( 的 求解器
ode113、ode15s、ode23s、ode23t、ode23tb)
问题, 没有一种算法可以有效地解决所有的 ODE 问题,因此 MATLAB 提供了多种ODE求解器,对于不同的ODE, 提供了多种 求解器,对于不同的 , 求解器 可以调用不同的求解器 求解器。 可以调用不同的求解器。
求微分方程的解
研究具有实际应用背景的数学函数 物理量在空间的分布情况其及随时间变化的规律。 物理量在空间的分布情况其及随时间变化的规律。
物理现象
物理定律
微分方程
数值解
建模型: 物理学家、化学家、生物学家、 建模型 物理学家、化学家、生物学家、社会学家
问题背景和实验目的
自牛顿发明微积分以来, 自牛顿发明微积分以来,微分方程在描述事物运 动规律上已发挥了重要的作用。实际应用问题通过 动规律上已发挥了重要的作用。 数学建模所得到的方程,绝大多数是微分方程。 数学建模所得到的方程,绝大多数是微分方程。 由于实际应用的需要,人们必须求解微分方程。 由于实际应用的需要,人们必须求解微分方程。 然而能够求得解析解的微分方程十分有限, 然而能够求得解析解的微分方程十分有限,绝大多 数微分方程需要利用数值方法来近似求解。 数微分方程需要利用数值方法来近似求解。 首先重点介绍一个求解微分方程的基本数值解法 --Euler折线法。 折线法。 -- 折线法
其中
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 方法源程序
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'''
k k
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
当 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-')
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 举例
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]);
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
相关主题