微分方程的数值解法matlab(四阶龙格—库塔法).
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
说明：
Solver解算指令的使用格式
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol); ode45
输出宗量形式
y1 (t0பைடு நூலகம்) y (t ) Y 1 1 y1 (t 2 )
ode23 ode113
ode15s
ode23tb
ode23s
一.解ODE的基本机理:
1. 列出微分方程
f ( y, y , t) y
(2.1)
令
y1 y, y2 y
(0) y 0 y(0) y0 , y
初始条件
2. 把高阶方程转换成一阶微分方程组 1 y f1 (t , Y ) Y f (t , Y ) y f ( t , Y ) 2 2
subplot (122), plot (Y( :, 1), Y( :, 2))
Van der Pol方程 % 子程序 （程序名： dYdt.m ）
function Ydot = dYdt (t, Y) Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];
或写为
function Ydot = dYdt (t, Y) Ydot=zeros(size(Y)); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];
%M function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y)); Yd (1) Y (2); Yd (2) (Y (1).^2 1) *Y (2) Y (1);
4. 使编写好的ODE函数文件和初值 Y0 供微分 方程解算指令（solver）调用
y2 (t0 ) y2 (t1 ) y2 (t 2 )
Y (:,1) y1 (t ) Y (:,2) y2 (t )
例题1：著名的Van der Pol方程
2 y0 y ( y 1) y
解法1：采用 ODE命令
% 主程序 （程序名：VanderPol _ex1.m） t0 = 0; tN = 20; tol = 1e-6; Y0 = [0.25; 0.0]; [t, Y]=ode45 (‘dYdt’, t0, tN, Y0, tol); subplot (121), plot (t, Y)
(2.2)
y1 (0) Y0 y ( 0 ) 2
(2.3)
例：著名的Van der Pol方程
2 y0 y ( y 1) y
令
y1 y, y2 y
y1 Y y2
降为一阶
1 y y 2 Y y 2 2 ( y1 1) y2 y1
微分方程的数值解法
四阶龙格—库塔法 （The Fourth-Order Runge－Kutta Method）
常微分方程（Ordinary differential equations, ODE)
初值问题---给出初始值 边值问题---给出边界条件
ode45 ode23t
与初值常微分方程解算有关的指令
单步法-Runge-Kutta 方法
多步法-Admas方法
计算 y(tn1 ) 的近似值 yn1 时只用到 t n , yn ，是自开始 方法
Runge-Kutta法是常微分方程的一种经典解法
MATLAB 对应命令：ode45
四阶Runge-Kutta公式
h yn 1 yn (k1 2k 2 2k3 k 4 ) 6 k1 f (t n , yn ) 1 h k 2 f (t n h, yn k1 ) 2 2 1 h k3 f (t n h, yn k 2 ) 2 2 k 4 f (t n h, yn hk3 )
各种solver 解算指令的特点
解法指令 ode45 ode23 ode113 ode23t ode15s ode23s 解题类 型 非刚性 非刚性 非刚性 适度刚 性 刚性 刚性 特 点 采用4、5阶Runge－Kutta法 采用Adams算法 多步法；采用Adams算法；高 低精度均可（10－3～10－6） 采用梯形法则算法 多步法；采用2阶Rosenbrock 算式，精度中等 一步法；采用2阶Rosenbrock 算式，低精度 采用梯形法则－反向数值微分 两阶段算法，低精度 适合场合 大多数场合的首选算法 较低精度（10－3）场合 ode45计算时间太长时 取代ode45 适度刚性 当ode45失败时使用； 或存在质量矩阵时 低精度时，比ode15s有 效；或存在质量矩阵时 低精度时，比ode15s有 效；或存在质量矩阵时
y1 (0) y10 Y0 y 2 (0) y 20
初始条件
3. 根据式(2.2)编写计算导数的M函数文件ODE文件
作为输出宗量 把t,Y作为输入宗量，把 Y
%M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展开式 例Van der Pol方程
ode23tb 刚性
二. 四 阶 Runge-Kutta 法
dy f (t , y) a t b dt y(t0 ) y0
对 I=[a,b]作分割 步长
a t0 t1 t N 1 b
hi hi 1 hi , i 0,1,, N 1
初值问题的数值 解法分为两大类
四 阶 Runge-Kutta 法计算流程图
开始 初始条件： t 0；y0 迭代次数： N 积分步长：
h
t n t0
for i = 1 : N
Solver解算指令的使用格式
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol); ode45
输出宗量形式
y1 (t0பைடு நூலகம்) y (t ) Y 1 1 y1 (t 2 )
ode23 ode113
ode15s
ode23tb
ode23s
一.解ODE的基本机理:
1. 列出微分方程
f ( y, y , t) y
(2.1)
令
y1 y, y2 y
(0) y 0 y(0) y0 , y
初始条件
2. 把高阶方程转换成一阶微分方程组 1 y f1 (t , Y ) Y f (t , Y ) y f ( t , Y ) 2 2
subplot (122), plot (Y( :, 1), Y( :, 2))
Van der Pol方程 % 子程序 （程序名： dYdt.m ）
function Ydot = dYdt (t, Y) Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];
或写为
function Ydot = dYdt (t, Y) Ydot=zeros(size(Y)); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];
%M function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y)); Yd (1) Y (2); Yd (2) (Y (1).^2 1) *Y (2) Y (1);
4. 使编写好的ODE函数文件和初值 Y0 供微分 方程解算指令（solver）调用
y2 (t0 ) y2 (t1 ) y2 (t 2 )
Y (:,1) y1 (t ) Y (:,2) y2 (t )
例题1：著名的Van der Pol方程
2 y0 y ( y 1) y
解法1：采用 ODE命令
% 主程序 （程序名：VanderPol _ex1.m） t0 = 0; tN = 20; tol = 1e-6; Y0 = [0.25; 0.0]; [t, Y]=ode45 (‘dYdt’, t0, tN, Y0, tol); subplot (121), plot (t, Y)
(2.2)
y1 (0) Y0 y ( 0 ) 2
(2.3)
例：著名的Van der Pol方程
2 y0 y ( y 1) y
令
y1 y, y2 y
y1 Y y2
降为一阶
1 y y 2 Y y 2 2 ( y1 1) y2 y1
微分方程的数值解法
四阶龙格—库塔法 （The Fourth-Order Runge－Kutta Method）
常微分方程（Ordinary differential equations, ODE)
初值问题---给出初始值 边值问题---给出边界条件
ode45 ode23t
与初值常微分方程解算有关的指令
单步法-Runge-Kutta 方法
多步法-Admas方法
计算 y(tn1 ) 的近似值 yn1 时只用到 t n , yn ，是自开始 方法
Runge-Kutta法是常微分方程的一种经典解法
MATLAB 对应命令：ode45
四阶Runge-Kutta公式
h yn 1 yn (k1 2k 2 2k3 k 4 ) 6 k1 f (t n , yn ) 1 h k 2 f (t n h, yn k1 ) 2 2 1 h k3 f (t n h, yn k 2 ) 2 2 k 4 f (t n h, yn hk3 )
各种solver 解算指令的特点
解法指令 ode45 ode23 ode113 ode23t ode15s ode23s 解题类 型 非刚性 非刚性 非刚性 适度刚 性 刚性 刚性 特 点 采用4、5阶Runge－Kutta法 采用Adams算法 多步法；采用Adams算法；高 低精度均可（10－3～10－6） 采用梯形法则算法 多步法；采用2阶Rosenbrock 算式，精度中等 一步法；采用2阶Rosenbrock 算式，低精度 采用梯形法则－反向数值微分 两阶段算法，低精度 适合场合 大多数场合的首选算法 较低精度（10－3）场合 ode45计算时间太长时 取代ode45 适度刚性 当ode45失败时使用； 或存在质量矩阵时 低精度时，比ode15s有 效；或存在质量矩阵时 低精度时，比ode15s有 效；或存在质量矩阵时
y1 (0) y10 Y0 y 2 (0) y 20
初始条件
3. 根据式(2.2)编写计算导数的M函数文件ODE文件
作为输出宗量 把t,Y作为输入宗量，把 Y
%M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展开式 例Van der Pol方程
ode23tb 刚性
二. 四 阶 Runge-Kutta 法
dy f (t , y) a t b dt y(t0 ) y0
对 I=[a,b]作分割 步长
a t0 t1 t N 1 b
hi hi 1 hi , i 0,1,, N 1
初值问题的数值 解法分为两大类
四 阶 Runge-Kutta 法计算流程图
开始 初始条件： t 0；y0 迭代次数： N 积分步长：
h
t n t0
for i = 1 : N