第7章 常微分方程初值问题的数值解法 - 图文

第7章 常微分方程初值问题的数值解法

7.1 引 言

科学研究和工程技术中的许多问题在数学上往往归结为微分方程的求解问题。为了确定微分方程的解，一般要加上定解条件，根据不同的情况，这些定解条件主要有初始条件(initial condition)和边界条件(boundary condition). 只含初始条件作为定解条件的微分方程求解问题称为初值问题(initial-value problem); 例如天文学中研究星体运动,空间技术中研究物体飞行等，都需要求解常微分方程初值问题(initial-value problem for ordinary differential equations). 只含边界条件作为定解条件的微分方程求解问题称为边值问题(boundary-value problem).

除特殊情形外，微分方程一般求不出解析解，即使有的能求出解析解，其函数表示式也比较复杂，计算量比较大，而且实际问题往往只要求在某一时刻解的函数值. 为了解决这个问题，有两种方法可以逼近原方程的解。第一种方法是：将原微分方程化简为可以准确求解的微分方程，然后使用化简后的方程的解近似原方程的解；第二种方法是：将求原微分方程的解析解转化为求原方程的数值解，这是实际中最常用的方法。

本章将介绍求解常微分方程初值问题的常用的数值方法。第8章将介绍常微分方程的边值问题的常用的数值方法。

为简明起见，本章主要介绍形如

?y?(t)?f(t,y), a?t?b, ? (7.1.1)

y(a)???的初值问题的数值解法. 在介绍这些方法之前，需要了解常微分方程的一些相关定义和结果.

2定义7.1.1 函数f(t,y)称为在集合D?R上关于变量y满足Lipschitz条件

(Lipschitz condition)，简称Lip条件，如果存在常数L?0，使得

f(t,y1)?f(t,y2)?Ly1?y2 (7.1.2)

对所有(t,y1),(t,y2)?D都成立. 常数L称为Lipschitz常数(Lipschitz constant)，简称Lip常数.

定义7.1.2 如果对所有(t1,y1),(t2,y2)?D，都有

?(1??)t1??y1,(1??)t2??y2)??D， (7.1.3)

2其中0???1，那么称集合D?R为凸集(convex set).

注 如果没有特别说明，本章中的集合D?R为D?{(t,y)a?t?b,???y??}，显然它是凸集。 定理7.1.1 设函数f(t,y)在D?{(t,y)a?t?b,???y??}上连续. 如果2f(t,y)在D上关于变量y满足Lipschitz条件，那么初值问题(7.1.1)对a?t?b有唯一解y(t).

1

2定理7.1.2 设函数f(t,y)定义在凸集D?R上. 如果存在常数L?0，使得

fy?(t,y)?L (7.1.4)

对一切(t,y)?D成立，那么f(t,y)在D上关于变量y满足Lipschitz条件，且L为Lipschitz常数. 注 因为Lipschitz条件一般不易验证，所以往往用条件(7.1.4)代替Lipschitz条件(7.1.2)。 由定理7.1.1立得： 定理7.1.3 设函数f(t,y)在D?{(t,y)a?t?b,???y??}上连续. 如果存在常数L?0，使得fy?(t,y)?L对一切(t,y)?D成立，，那么初值问题(7.1.1)对a?t?b有唯一解y(t). 理论上，初始值是准确的，微分方程的计算过程也是准确的；而实际情况并非如此. 事实上，初始数据可能有误差，计算过程中的数字的舍入也会产生误差。对初值问题(7.1.1)，若初值?有误差?0，即实际得到的初值是???0，再考虑到计算过程中的数字的舍入也会产生误差，则得到的初值问题变为 ?z?(t)?f(t,z)??(t), a?t?b, (7.1.5) ?z(a)????,0?这种误差的扰动在计算过程中是否会增长很快，以至影响结果，这就是初值问题的适定性问题. 定义7.1.3 设初值问题(7.1.1)和(7.1.5)的均有唯一解，分别是y?y(t)和z?z(t)，且?(t)连续。若对任意??0，存在常数K?K(?)?0，当?0??,?(t)??时，对任意t?[a,b]，恒有y(t)?z(t)?K?，则称初值问题(7.1.1)是适定的(well-posed)。 定理7.1.6 设函数f(t,y)在D?{(t,y)a?t?b,???y??}上连续. 如果f(t,y)在D上关于变量y满足Lipschitz条件，那么初值问题(7.1.1)是适定的。 注 初值问题只有是适定的才有意义，因此本章考察的所有初值问题都是适定的。 求微分方程数值解的主要问题：

(1) 如何将原微分方程离散化, 并建立求其数值解的递推公式;

(2) 如何求递推公式的局部截断误差, 数值解yn与精确解y(tn)的误差估计； (3) 研究递推公式的稳定性与收敛性.

2

7.2 Euler方法及改进的Euler方法

7.2.1 Euler格式与梯形格式

考虑一阶常微分方程的初值问题

?y?(t)?f(t,y), a?t?b, (7.2.1) ?y(a)??.?设a?t0?t1??tN?1?tN?b，其中tn?t为等距节点，步长?0,?1,,N0?n,hnh?(b?a)N.

在[tn,tn?1](n?0,1,?,N?1)上对y??f(t,y(t))两边积分，得

y(tn?1)?y(tn)??7.2.1.1 Euler格式

tn?1tnf(t,y(t))dt (7.2.2)

用左矩形求积公式计算式(7.2.2)右端积分项，得

?代入式(7.2.2)右端，得

tn?1tnf(t,y(t))dt?hf(tn,y(tn)),

y(tn?1)?y(tn)?hf(tn,y(tn)). (7.2.3)

用y(tn)的近似值yn代入式(7.2.3)右端，记所得结果为yn?1，得

?yn?1?yn?hf(tn,yn),n?0,1,?,N?1, (7.2.4) ?y??,?0并称式(7.2.4)为求解初值问题(7.2.1)的Euler方法(Euler’s method)或Euler格式(Euler scheme)，yn?1?yn?hf(tn,yn)称为差分方程(difference equation).

注 Euler方法是最早的解决一阶常微分方程初值问题的一种数值方法，虽然它的精度不高，很少被采用，但是它反映了微分方程的数值解法的基本思想和特征.

若式(7.2.2)右边的积分由数值积分的右矩形公式近似，并用近似值yn替代y(tn)，近似值yn?1替代y(tn?1)，则可得到

?yn?1?yn?hf(tn?1,yn?1),n?0,1,?,N?1, (7.2.5) ??y0??,并称式(7.2.5)为后退的Euler方法(backward Euler’s method) 或后退的Euler格式(backward Euler’s scheme). yn?1?yn?hf(tn?1,yn?1)是差分方程。

注 在xOy平面上，微分方程y??f(t,y)的解y?y(t)称为积分曲线，积分曲线上一点(t,y)的切线斜率等于函数y?(t)?f(t,y)在x点的值. 如果在

D?{(t,y)?a?tb,???中每一点y?}?(t,y)都画一条以f(t,y)在点(t,y)的值为斜

——————————————————

欧拉(Leonhard Paul Euler 1707年4月15日 –1783年9月18日) 是瑞士杰出的数学家，自然科学家. 他在数学分析，图论，力学，光学和天文学等很多领域都做出了卓越的贡献，被公认为有史以来最伟大的数学家之一.

率并指向t增加方向的有向线段(即：在D上作出了一个由方程y??f(t,y)确定的方向场），

3

那么从几何上看，微分方程y??f(t,y)的解y?y(t)就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向相一致.

?从初始点P0(a,y0)出发，过这点的积分曲线为y?y(t)，斜率为y0?f(a,y0). 设在

切线方程为y?y0?f(a,y0)(t?a). 当t?t1t?a附近y(t)可用过P0点的切线近似表示，

时，y(t1)的近似值为y0?f(a,y0)(t1?a)，并记为y1，这就是得到t?t1时计算y(t1)的近似公式

y1?y0?f(a,y0)(t?a)

当t?t2时，y(x2)的近似值为y1?f(t1,y1)(t2?t1)，并记为y2. 于是就得到当t?t2时计算

y(t2)的近似公式

y2?y1?f(t1,y1)(t2?t1)

重复上面方法，一般可得当t?tn?1的计算y(tn?1)的近似公式

yn?1?yn?f(tn,yn)(tn?1?tn)

如果h?tn?tn?1(n?1,2,?,N)，则上面公式就是式(7.2.4). 将P0,P就得1,?,PN连起来，到一条折线，所以Euler方法又称为折线法(polygon method), 见图7.2.1.

由公式(7.2.4)知，已知y0便可算出y1；已知y1，便可算出y2；如此继续下去，这种只用前一步的值yk便可计算出yk?1的递推公式称为单步法(one-step method).

图7.2.1

7.2.1.2 梯形格式

用梯形求积公式计算式(7.2.2)右端积分项，得 tn?1?tnhf(t,y(t))dt?[f(tn,y(tn))?f(tn?1,y(tn?1))], 2代入式(7.2.2)右端，得

hy(tn?1)?y(tn)?[f(tn,y(tn))?f(tn?1,y(tn?1))]. (7.2.6)

2在式(7.2.6)中，将y(tn)用yn近似替代，所得结果记为yn?1，得

h?y?y?[f(tn,yn)?f(tn?1,yn?1)],n?0,1,?,N?1,?n?1n (7.2.7) 2???y0??,

4

并称(7.2.7)为求解初值问题(7.2.1)的梯形方法(trapezoidal method), 或梯形格式(trapezoidal scheme).yn?1?yn?h[f(tn,yn)?f(tn?1,yn?1)]是差分方程. 27.2.2 改进的Euler方法

因为梯形求积公式的精度比矩形求积公式的精度高，所以求解初值问题(7.2.1)的梯形方法比Euler方法精度高；但梯形方法是隐式方法，不便于计算。为了克服这个困难，先用Euler格式(7.2.4)求得一个初步的近似值yn?1，称之为预测值；然后用公式(7.2.6)作一次迭代得

yn?1，即将yn?1校正一次，得如下方法：

定义7.2.1 对n?0,1,?,N?1，称预测－校正系统(predictor-corrector system)

?预测：yn?1?yn?hf(tn,yn),? (7.2.8) ?h校正：yn?1?yn?[f(tn,yn)?f(tn?1,yn?1)]??2为改进的Euler格式(modified Euler scheme), 称这种方法为求解初值问题(7.2.1)的改进的Euler方法(modified Euler method).

注 改进的Euler格式(7.2.8)还可以表示为下列平均化形式：

?yp?yn?hf(tn,yn),???yc?yn?hf(tn?1,yp), (7.2.9) ?1??yn?1?2(yp?yc).例7.2.1 取步长h?0.1，分别用Euler方法及改进的Euler方法求解初值问题

2?2t?y?(t)?y?te,1?t?2, t???y(1)?0.2t解 这个初值问题的准确解为y(t)?t(e?e). 根据题设知f(t,y)?2y?t2et. tEuler方法的计算格式为 2tnyn?1?yn?0.1?t2y?t. nenn??改进的Euler方法的计算格式为 2tn?yn?1?yn?0.1?2yn?tne,tn? ?tn?1?2tn20.1?22?yn?1?yn?2??tnyn?tne?tyn?1?tn?1e?,n?1?????????或

2tn?yp?yn?0.1?2yn?tne,tn??t22?yc?yn?0.1?tn?1yp?tn?1en?1, ??yn?1?12(yp?yc).????? 5

fun = inline('(2./t).*y+(t.^2).*exp(t)','t','y'); ImprAdams_4(fun, 1, 2, 0.1, 0) 回车，可得结果。

7.7.6 一阶常微分方程组的Runge-Kutta法 %一阶常微分方程组的四阶Runge-Kutta法 %[a, b]是区间, Step是步长, Alpha, Beta是初值

% f, g分别是式(7.6.3)中的f (t, y, z), g (t, y, z)

function [Y Z] = RungeKutta_Sys (f, g, a, b, Step, Alpha, Beta) N = (b-a) / Step; T = a : Step : b; Y(1) = Alpha; Z(1) = Beta; for i = 2 : N+1

K1 = feval(f, T(i-1), Y(i-1), Z(i-1));

M1 = feval(g, T(i-1), Y(i-1), Z(i-1));

K2 = feval(f, T(i-1)+Step/2, Y(i-1)+Step/2*K1, Z(i-1)+Step/2*M1); M2 = feval(g, T(i-1)+Step/2, Y(i-1)+Step/2*K1, Z(i-1)+Step/2*M1); K3 = feval(f, T(i-1)+Step/2, Y(i-1)+Step/2*K2, Z(i-1)+Step/2*M2); M3 = feval(g, T(i-1)+Step/2, Y(i-1)+Step/2*K2, Z(i-1)+Step/2*M2);

K4 = feval(f, T(i-1)+Step, Y(i-1)+Step*K3, Z(i-1)+Step*M3); M4 = feval(g, T(i-1)+Step, Y(i-1)+Step*K3, Z(i-1)+Step*M3); Y(i) = Y(i-1) + Step/6 * (K1+2*K2+2*K3+K4); Z(i) = Z(i-1) + Step/6 * (M1+2*M2+2*M3+M4);

sprintf('t=%3.1f,y=%9.7f,z=%9.7f', T(i), Y(i), Z(i)) end end

例7.7.6 使用微分方程组的四阶Runge-Kutta法求下列一阶微分方程组的近似解??y?(t)?3y1(t)?2y2(t)?(2t2?1)e2t1,0?t?1,?y?(t)?4y1(t)?y2(t)?(t2?2t?4)e2t,0?t?1, ?2?y1(0)?y2(0)?1,取步长h?0.1.

解 在MATLAB命令窗口输入：

f = inline('3.*y+2.*z-(2.*t.^2+1).*exp(2.*t)', 't', 'y', 'z'); g = inline('4.*y+z+(t.^2+2.*t-4).*exp(2.*t)', 't', 'y', 'z'); RungeKutta_Sys (f, g, 0, 1, 0.1, 1, 1) 回车，可得结果。

7.7.7 用四阶Runge-Kutta法解高阶微分方程初值问题 %用四阶Runge-Kutta法解高阶微分方程初值问题 %[a, b]是区间, Step是步长, Alpha, Beta是初值

41

% f分别是式(7.6.25)中的f (t, y, z)

function Y = RungeKutta_HOrder (f, a, b, Step, Alpha, Beta) N = (b-a) / Step; T = a : Step : b; Y(1) = Alpha; Z(1) = Beta; for i = 2 : N+1

M1 = feval(f, T(i-1), Y(i-1), Z(i-1));

M2 = feval(f, T(i-1)+Step/2, Y(i-1)+Step/2*Z(i-1), Z(i-1)+Step/2*M1); M3 = feval(f, T(i-1)+Step/2, Y(i-1) + Step/2*Z(i-1)+ ... Step^2/4*M1, Z(i-1)+Step/2*M2);

M4 = feval(f, T(i-1)+Step, Y(i-1) + Step*Z(i-1)+ ... Step^2/2*M2, Z(i-1)+Step*M3);

Y(i) = Y(i-1) + Step*Z(i-1) + Step^2 / 6 * (M1 + M2 + M3); Z(i) = Z(i-1) + Step / 6 * (M1 + 2*M2 + 2*M3 + M4); sprintf('t=%3.1f,y=%9.7f', T(i), Y(i)) end end

例7.7.7 使用微分方程组的Runge-Kutta法求下列高阶微分方程的近似解

?y???2y??y?tet?t,0?t?1, ???y(0)?y(0)?0,取h?0.1.

解 在MATLAB命令窗口输入： f = inline('2.*z-y+t.*exp(t)-t','t', 'y', 'z'); RungeKutta_HOrder (f, 0, 1, 0.1, 0, 0) 回车，可得结果。

小 结

本章介绍了常微分方程初值问题的常用的数值解法. 利用数值微积分法构造常微分方程初值问题的数值解法: 将微分方程离散化，建立数值解的递推公式. 然后利用Taylor展开式研究这些数值方法的相容性、收敛性和稳定性。

由递推公式的结构来看，可分为单步法和多步法. 也可以分为显式方法和隐式方法. 单步法中，我们主要介绍了Euler方法、梯形方法和Runge-Kutta方法等，其中Euler方法和Runge--Kutta方法是显式方法，而梯形方法是隐式方法. 多步法中，我们主要介绍了Adams方法；其中既有显式Adams方法，也有隐式Adams方法. 一般地，显式方法计算简单，隐式方法计算较复杂, 但稳定性比同阶显式方法好.

Euler方法是最早的解决一阶常微分方程初值问题的数值方法，虽然它的精度不高，很少被采用，但是它反映了解微分方程的数值方法的基本思想和特征，后来的微分方程的数值方法无论从构造，还是理论分析，都借鉴了Euler方法的思想.

在实际问题中，如果要求解的精度高，常用的是四阶Runge-Kutta方法，这是因为四阶Runge-Kutta方法精度高，编程简单，易于调节步长，计算过程稳定；其缺点是计算量较大，

42

计算函数f(x,y)的次数比同阶的线性多步法（四阶Adams方法）多几倍；另外，如果f(x,y)的光滑性较差，Runge-Kutta方法的精度还不如改进的Euler方法高.

一般地，如果f(x,y)比较复杂，那么一般采用四阶Adams方法. 因为Adams方法是多步法，所以通常用四阶Runge-Kutta方法求出启动值，用四阶显式Adams方法作出预测值，再用四阶隐式Adams方法进行校正；其中配合使用误差修正，这就是通常使用的修正的四阶Adams预测校正方法.

不论采用哪种方法，选取的步长h应使?h落在绝对稳定区域内. 一般在保证精度的条件下，步长尽可能大些，这样可以节省计算量.

上述数值方法都能直接用来解一阶常微分方程组初值问题. 而高阶常微分方程初值问题可化为一阶常微分方程组初值问题来解.

习 题 七

1. 取步长h?0.1，用Euler方法求解下列初值问题，并分别与它们的准确解进行比较(小数点后保留5位数字).

?y?(t)??y(t)(1?ty(t)),0?t?1,(1) ? (2)

?y(0)?1,?y?(t)?1?(t?y(t))2,2?t?3, ??y(2)?1.其中(1)的准确解是y(t)?1(2et?t?1)和(2)的准确解是y(t)?t?1(1?t). 2. 取步长h?0.1，分别用改进的Euler方法、中点方法和Heun方法求解下列初值问题，并分别与它们的准确解进行比较(小数点后保留7 位数字).

?y?(t)??y(t)(1?ty(t)),0?t?1,(1) ? (2)

?y(0)?1,?y?(t)?1?(t?y(t))2,2?t?3, ??y(2)?1.其中(1)的准确解是y(t)?1(2et?t?1)和(2)的准确解是y(t)?t?1(1?t). 3. 取步长h?0.1，分别用三阶Runge-Kutta方法(7.3.13)四阶经典Runge-Kutta方法(7.3.14)求解下列初值问题，并分别与它们的准确解进行比较(小数点后保留7位数字).

?y?(t)??y(t)(1?ty(t)),0?t?1,(1) ? (2)

?y(0)?1,?y?(t)?1?(t?y(t))2,2?t?3, ??y(2)?1.4. 取步长h?0.1，用四阶Admas预测-校正方法及修正的Admas预测-校正方法求解下列初值问题，并分别与它们的准确解进行比较(小数点后保留7位数字).

?y?(t)??y(t)(1?ty(t)),0?t?1,(1) ? (2)

y(0)?1,?5. 用改进的Euler方法求初值问题

?y?(t)?1?(t?y(t))2,2?t?3, ??y(2)?1. 43

(1) ??y?(t)?y?0,?5 的解y(t)在t?0.2处的近似值，要求迭代误差不超过10;

?y(0)?1?y?(t)?y(t)?y2(t)sint?0,(2) ?的解y(t)在t?1.2,t?1.4处的近似值，取步长

?y(1)?1h?0.2，要求小数点后保留5位数字。

6. 使用常微分方程组的改进的Euler方法(7.6.4)求下列一阶微分方程组的近似解。

?y?(t)?3t?2z(t),??z?(t)?2y(t)?z(t), ?y(0)?1,z(0)?1.?取h?0.1，计算y(0.1)和z(0.1)的近似值。

7. 使用常微分方程组的Runge-Kutta法求下列一阶微分方程组的近似解，并将结果与精确解进行比较.

?y?(t)??4y(t)?2z(t)?cost?4sint,??z?(t)?3y(t)?z(t)?3sint,?y(0)?0,z(0)??1,?0?t?2,0?t?2,

取h?0.2，精确解是y(t)?2e?t?2e?2t?sint和z(t)??3e?t?2e?2t.

8. 使用常微分方程组的Runge-Kutta法求下列高阶微分方程的近似解，并将结果与精确解进行比较.

?y???2y??2y?e2tsint,0?t?1, ??y(0)??0.4,y?(0)??0.6,取h?0.1，精确解是y(t)?0.2e(sint?2cost).

注：下面的9—13题都是针对下列初值问题

2t?y?(t)?f(t,y),a?t?b (*) ??y(a)??,且数值方法的步长都是h?(b?a)N, 网格点都是tn?a?nh,的初值都是y0?y(a)??.

9. 设求解初值问题(*)的计算格式为

n?0,1,?,N，计算格式?yn?1?yn?h(k2?k3),n?0,1,?N?1,2??k1?f(tn,yn), ??k2?f(tn?qh,yn?qhk1),?k?f(t?(1?q)h,y?(1?q)h).nn?3 44

证明上述计算格式对任意q?(0,1)都是二阶格式。

10. 设求解初值问题(*)的计算格式为

yn?1?yn?h[?f(tn,yn)??f(tn?1,yn?1)],n?1,2,?N?1.

若已知y(tn)?yn,y(tn?1)?yn?1，试确定参数?和?，使上述计算格式的局部截断误差阶为O(h3).

11. 设求解初值问题(*)的计算格式为

yn?1??0yn??1yn?1?h[?0f(tn,yn)??1f(tn?1,yn?1)],n?1,2,?N?1.

若已知y(tn)?yn,y(tn?1)?yn?1，试确定参数?0,式.

12. 设求解初值问题(*)的计算格式为

?1,?0和?1，使上述计算格式为三阶格

yn?1?yn?h[23f(tn,yn)?16f(tn?1,yn?1)?5f(tn?2,yn?2)],n?2,3,?N?1. 12试推导上述计算格式的局部截断误差；并指出它是几阶格式. 13. 设求解初值问题(*)的的预报-校正公式为 ?yn?1?yn?hf(tn,yn),? ?hyn?1?yn?[5f(tn?1,yn?1)?8f(tn,yn)?5f(tn?1,yn?1)],n?1,2,?N?1.??12试推导上述预报-校正公式的局部截断误差；并指出它是几阶格式. 14. 设f(t,y)关于y的偏导数存在，且满足fy?(t,y)?L. (1) 证明迭代格式

(0)?yn?yn?hf(tn,yn),??1 ?(k?1)h(k)f(tn,yn)?f(tn?1,yn?1)?,k?0,1,2,?,?yn?1?yn?????2当

hL?1时收敛。 2(2) 用(1)中的迭代格式求初值问题 ?dy?etsin(ty),0?t?1,?dt ???y(0)?1.如何选取步长h，才能保证上述迭代格式收敛。

15. 证明定理7.4.2. 16. 证明定理7.4.3.

45

本文来源：https://www.bwwdw.com/article/72xf.html