计算机应用基础-5-常微分方程求解

2 三阶龙格-库塔
h y y ( k1 4 k 2 k 3 ) n 1 n 6 k1 f ( x n , y n ) k f ( x h , y h k ) n n 1 2 2 2 k f ( x h, y hk 2hk ) n n 1 2 3
h y n 1 y n k1 3k 2 3k 3 k 4 8 k1 f x n , y n 1 1 k f x h , y hk 2 n n 1 3 3 2 1 k f x h , y hk hk n 3 n 1 2 3 3 k 4 f x n h, y n hk1 hk2 hk3
y '' x F ( x, y, y ' )
在区间[a;b]上求解该方程的 解，且在这两个边界点上满足：
a ya a y ' a a b yb b y ' b b
上式即为边界条件
5.3 常微分方程的边值问题求解
对于线性微分方程
求得。类似地，由
y ( xn 1 ) y ( xn ) y' ( xn ) 以及 h y ' ( xn ) f ( xn , y ( xn ))
得到计算近似值的向前欧拉公式
yn1 yn hf ( xn , yn )
5.2 常微分方程的初值问题求解
由yn直接算出yn+1值的计算格式称为显式格式， 向前欧拉公式是显式格式。而欧拉方法的几何意义是 以y1作为斜率，通过点(x0, y0)做一条直线，它与直线 x=x1的交点就是y1。依此类推，是以yn+1作为斜率，经 过点的直线与直线x=xn+1的交点。故欧拉法也称为欧 拉折线法，如下图所示。
y '' x pxy ' x qxyx f x
5.2 常微分方程的初值问题求解
3、四阶龙格—库塔公式
h y y n1 n k1 2k 2 2k 3 k 4 6 k1 f x n , y n 1 1 ( 2) k 2 f x n h, y n hk1 2 2 1 1 k f x h , y hk2 3 n n 2 2 k 4 f xn h, y n hk3
y: 状态变量，求解
5.2 常微分方程的初值问题求解
例5-1 求解初值问题：
2x ' y y y y0 1
0 x 1
求解程序”xODE.m”
计算结果： ode45 x=1
ode23
精确值
1.7320508
1.73205081676653 1.73215487941780
(1)
(2)
(3)
h 2k1 3k 2 4k 3 y y n n 1 9 k1 f x n , y n k f x 1 h, y 1 hk n n 1 2 2 2 3 3 k f x h , y hk n n 2 3 4 4
图4.1
欧拉折线法
5.2 常微分方程的初值问题求解 二 龙格-库塔方法
龙格-库塔法是求解常微分方程较常用的一种方法，它通 过巧妙的线性组合，在显式格式的情况下获得理想的计算精度， 大大提高了计算速度。该方法的推导过程不是本课程要研究的 对象，本课程主要研究其实际应用,，故直接给出各类龙格-库 塔公式。 1、 二阶龙格-库塔
5.2 常微分方程的初值问题求解
x dy y dx y (0) 1
下面用一个例子说明：
用标准四阶龙格-库塔方法计算得
步长 h=1 h=0.5 h=0.25 h=0.125 y(2) 2.23940 2.23620 2.23607 2.23607 计算次数 2 4 8 16
把它与上式相比，可得
( h / 2) y ( xn 1 ) y n 1 1 ( h) 16 y ( xn 1 ) y n 1
经整理可得
( h / 2) ( h) 1 y( xn1 ) yn [ y ( x ) y 1 n1 n1 ] 15
的近似值，其误差可用先后两次计算 ( h / 2) y ( x ) y n 1 n 1 结果之差来表示，因而，只需考察 是否 ( h / 2) 成立。若成立，则可将 y n1 作为y(xn+1)的近似值；若不成立， 则将步长再次对折进行计算，直到不等式成立为止，并取最后 ( h / 2) 的 yn 1 作为计算结果。以上方法就是计算过程中自动选择步长 的方法，也称为变步长方法。从表面上看，为了选择适当的步 长，每一步的计算量增加了，但从总体考虑，工作量往往还是 减少的。
第五章 常微分方程的求解
5.1 常微分方程基本特征 5.2 常微分方程的初值问题求解 5.3 常微分方程的边值问题求解 5.4 高阶常微分方程的求解
5.1 常微分方程的数值解
一般微分方程式描述系统内部变量的变化率如何受系 统内部变量和外部激励，如输入的影响。 当常微分方程式能够解析求解时，可用 MATLAB 的 符号工具箱中的功能找到精确解。 在微分方程难以获得解析解的情况下，可以方便地在 数值上求解。 微分方程数值解是数值计算的基本内容，由于微分方 程的多样性，有不同的解法。MATLAB给出的7种解 法分别由7个不同的函数来完成。
5.2 常微分方程的初值问题求解 龙格-库塔方法——步长的选择 怎样选取合适的步长，这在实际计算中是很重 要的。由于步长愈小，每步计算的截断误差就愈小； 但在一定的求解范围内，需要完成的步数就愈多， 这不但引起计算量的增加，而且计算步数的增加往 往造成舍入误差的严重积累，反而降低了计算精度。 上面介绍的龙格 - 库塔方法是对定步长（即步长 h 为 常数）而言的，但为了保证精确度，一种有效的措 施是在计算过程中自动进行步长调整，即变步长技 巧。
5.2 常微分方程的初值问题求解
二阶龙格-库塔公式的局部截断误差为O(h3), 是二阶 精度的计算公式。也可建立高阶的龙格-库塔公式， 如三阶龙格-库塔公式、四阶龙格-库塔公式。较常用 的是四阶龙格-库塔公式，四阶龙格-库塔公式的局部 截断误差界为O(h5), 是四阶精度的计算公式。
5.2 常微分方程的初值问题求解
由积分法算得y(2)=2.23607,当h=0.5时相差0.00013，而h=0.25误 差小于0.000001，但当 h=0.125时虽然足够精确但计算次数却比 h=0.25多了一倍。因此合理选择步长既能保证精度又能减少计算量。
5.3 常微分方程的边值问题求解 边值问题求解之一-打靶算法
二阶微分方程的边值问题的 求解：
(1)
5.2 常微分方程的初值问题求解
龙格-库塔方法的调用格式：
[t, y]=ode45(@ODEfun, tspan, y0, options, p1…)
求解器 自定义的函数 积分限或离 散的向量 初值
Options：积分参数值，不选为默认值 p1: 传递给ODEfun的参数 输出： t: 时间变量，因变量，列向量
5.1 常微分方程的数值解 一、刚性问题(Stiff Problem)
1、刚性方程
对于一个常微分方程组，如果其 Jacobian 矩阵的特征 值相差十分悬殊，那么这个方程组就称为刚性方程组。
对于刚性方程组，为保持解法的稳定，步长选取很 困难。有些解法不能用来解刚性方程组，而有些解 法对稳定性的要求不严格，可以用来解决刚性问题。
h y y ( k1 k 2 ) n n 1 2 k1 f ( x n , y n ) k f ( x h, y hk ) n n 1 2
y n 1 y n hk2 k1 f ( x n , y n ) h h k 2 f ( x n , y n k1 ) 2 2
( h / 2) y 这表明以 n1 作为y(x
n+1)
5.2 常微分方程的初值问题求解
龙格-库塔方法的主要优点是计算精确度较高，能满足通常 的计算要求；容易编制程序；每次计算y(xn+1) ，只用到前 一步的计算结果yn，因此，在已知初始值y0的条件下，就可 自动地进行计算，是单步法，而且计算过程中可随时改变 步长。它的缺点是每前进一步需要计算多次f(x,y)的值，因 此，计算工作量较大，且其截断误差难以估计。 在实际应用上，一般当要求更高精确度时，采用的办 法是缩小步长，而不是采用更高阶的公式，因为高阶公式 的计算太复杂，一般选用标准四阶龙格-库塔方法即可。
一 欧拉公式
y' ( x0 ) f ( x0 , y( x0 ))
y ( x1 ) y ( x0 ) f ( x0 , y ( x0 )) h
于是得到：
y ( x1 ) y ( x 0 ) y ' ( x0 ) h
5.2 常微分方程的初值问题求解
故y(x1)的近似值y1可按
y1 y 0 f ( x0 , y 0 )或y1 y 0 hf ( x0 , y 0 ) h
dy f x, y , a x b dx x b, x a, ga, y 0;
边值问题求解：
gb, y 0;
打靶法：将边值问题转化成初值问题
松弛法：包括有限差分法和配置法
5.2 常微分方程的初值问题求解
对于常微分方程初值问题,在求解区间[a,b]上作等距分割，步长 ba h m 记xn = xn-1+h , n =1,2…,m。用差商近似导数计算常微分方程。 y ' ( x) 做 的在x = x0处的一阶向前差商得：
常微分方程（组）包括初值问题和边值问题两种： （1） 初值问题： 初始条件已知，形式如下