龙格-库塔方法求解常微分方程（十五）

求解常微分方程—龙格-库塔方法

龙格-库塔方法时求解常微分方程的常用方法之一，其相关的函数命令如表所示。

​ 表1 龙格-库塔命令

命令	说明
ode23	二阶、三阶R-K函数，求解非刚性微分方程的低阶方法
ode45	四阶、五阶R-K函数，求解非刚性微分方程的中阶方法
ode113	求解更高阶或大的标量计算
ode15s	采用多步法求解刚性方程，精度较低
ode23s	采用单步法求解刚性方程，速度较快
ode23t	用于解决难度适中的问题
ode23tb	用于解决难度较大的问题，对于系统的存在常量矩阵的情况有用

注：刚性微分方程与非刚性方程：

在数学领域中，刚性方程（stiffness equation）是指一个微分方程，其数值分析的解只有在时间间隔很小时才会稳定，只要时间间隔略大，其解就会不稳定。目前很难去精确地去定义哪些微分方程是刚性方程，然而粗略而言，若此方程式中包含使其快速变动的项，则其为刚性方程。
在积分微分方程时，若某一区域的解曲线的变化很大，会希望在这个区域的积分间隔密一些，若另一区域的曲线近似直线，且斜率接近零，会希望在这个区域的积分间隔松一些。不过针对一些问题，就算曲线近似直线，仍然需要用非常小的积分间隔来积分，这种现象称为“刚性”。有时可能会出现两个不同问题，一个有“刚性”，另一个没有，但两个问题却有同一个解的情形。因此“刚性”不是解本身的特性，而是微分方程的特性，也可以称为是刚性系统。

特性：刚性系统的特色是该系统所有特征值的实部均为负数，并且其中特征值实部绝对值中，最大和最小的比值远大于1。 |

以上各种函数命令的调用格式如下：

1.[T,Y] = solver(odefun,tspan,y0)

2.[T,Y] = solver(odefun,tspan,y0,options)

3.[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options) options中事件属性要设为on

4.sol = solver(odefun,[t0,tf],y0…)

其中，odefun定义了微分方程的形式，tspan = [t0 tfinal]定义微分方程的积分限，y0是初始条件。

options参数的设置需要利用odeset函数命令，调用格式如下所示。

调用格式	说明
options = odeset(‘name1’,‘value1’,‘name2’,‘value2’,…)	创建一个参数结构，对指定的参数名进行设置，未设置的参数将使用默认值
options = odeset(oldopts,‘name1’,‘value1’,…)	对已有的参数结构进行修改
options = odeset(oldopts,newopts)	将已有的参数结构完整转换为newopts
odeset	显示所以参数的可能值与默认值

例 利用R-K方法求解方程$\{\begin{array}{c}{y}^{\prime }=2t\\ y(0)=0\end{array}(0<x<5)$

tspan = [0 5];
y0 = 0;
[t,y] = ode45(@(t,y) 2*t,tspan,y0)

t =

     0
0.1250
0.2500
0.3750
0.5000
0.6250
0.7500
0.8750
1.0000
1.1250
1.2500
1.3750
1.5000
1.6250
1.7500
1.8750
2.0000
2.1250
2.2500
2.3750
2.5000
2.6250
2.7500
2.8750
3.0000
3.1250
3.2500
3.3750
3.5000
3.6250
3.7500
3.8750
4.0000
4.1250
4.2500
4.3750
4.5000
4.6250
4.7500
4.8750
5.0000

y =

  0
0.0156
0.0625
0.1406
0.2500
0.3906
0.5625
0.7656
1.0000
1.2656
1.5625
1.8906
2.2500
2.6406
3.0625
3.5156
4.0000
4.5156
5.0625
5.6406
6.2500
6.8906
7.5625
8.2656
9.0000
9.7656 
10.5625
11.3906
12.2500
13.1406
14.0625
15.0156
16.0000
17.0156
18.0625
19.1406
20.2500
21.3906
22.5625
23.7656
25.0000

plot(t,y,'-o')

精度求解结果：