下载文档原格式
龙哥库塔法or欧拉法求解微分⽅程matlab实现举例：分别⽤欧拉法和龙哥库塔法求解下⾯的微分⽅程我们知道的欧拉法（Euler）"思想是⽤先前的差商近似代替倒数"，直⽩⼀些的编程说法即：f(i+1)=f(i)+h*f(x,y)其中h是设定的迭代步长，若精度要求不⾼，⼀般可取0.01。
在定义区间内迭代求解即可。
龙哥库塔法⼀般⽤于⾼精度的求解，即⾼阶精度的改进欧拉法，常⽤的是四阶龙哥库塔，编程语⾔如下：y(i+1)=y(i)+h*(k1+2*K2+2*k3+k4)/6;k1=f(xi,yi)k2=f(xi+h/2,yi+h*k1/2);k3=f(xi+h/2,yi+h*k2/2);k4=f(xi+h,yi+h*k3);设置终⽌条件迭代求解。
matlab实现程序如下：%% 龙哥库塔or欧拉法求解微分⽅程t=0:0.01:3; %⾃变量范围f = [];g = [];f(1) = 0.1; %f初值g(1) = 1; %g初值h=0.001;for i=1:length(t)%% 欧拉法% f(i+1) =f(i)+h*(exp(f(i)*sin(t(i)))+g(i));% g(i+1) =g(i)+h*(exp(g(i)*cos(t(i)))+f(i));%% 龙哥库塔kf1 = exp(f(i)*sin(t(i)))+g(i);%g(i)相当于常值kf2 = exp((f(i)+kf1*h/2)*sin(t(i)+h/2))+g(i);kf3 = exp((f(i)+kf2*h/2)*sin(t(i)+h/2))+g(i);kf4 = exp((f(i)+kf3*h)*sin(t(i)+h))+g(i);f(i+1) = f(i) + h*(kf1+2*kf2+2*kf3+kf4)/6;kg1 = exp(f(i)*cos(t(i)))+f(i);%f(i)相当于常值kg2 = exp((f(i)+kg1*h/2)*cos(t(i)+h/2))+f(i);kg3 = exp((f(i)+kg2*h/2)*cos(t(i)+h/2))+f(i);kg4 = exp((f(i)+kg3*h)*cos(t(i)+h))+f(i);g(i+1) = g(i) + h*(kg1+2*kg2+2*kg3+kg4)/6;endfigure(1)plot(t,f(1:length(t)),t,g(1:length(t)));legend('f函数','g函数')。
欧拉方法matlab
欧拉方法是一种数值解微分方程的方法，适用于一阶常微分方程和某些高阶微分方程。
本文将介绍如何使用MATLAB实现欧拉方法。
首先，我们需要定义微分方程和初始条件。
例如，对于一阶常微分方程dy/dx = f(x,y)，初始条件为y(x0) = y0，我们可以定义函数f和初始值x0和y0：
function dydx = f(x,y)
dydx = ... % 定义f(x,y)的具体表达式
end
x0 = ... % 初始值x0
y0 = ... % 初始值y0
接下来，我们需要设置步长和计算区间。
步长越小，计算结果越精确，但计算量也会增加。
计算区间可以根据需要设置。
h = ... % 步长
xspan = ... % 计算区间
然后，我们可以使用欧拉方法计算微分方程的数值解。
欧拉方法的基本思想是在每个步长上使用当前值和导数的近似值（如斜率）来估计下一个值。
首先，我们可以定义一个数组来存储计算结果，设初始值为y0： y = [y0];
然后，我们可以使用for循环来计算每个步长上的结果。
在每个步长上，我们需要计算当前值的导数，以及使用欧拉方法计算下一个
值。
龙格库塔法（Runge-Kutta method）和欧拉法（Euler's method）是两种常用的数值求解常微分方程的方法。
这里分别给出它们的MATLAB实现：1. 龙格库塔法（Runge-Kutta method）：```matlabfunction [y, t] = runge_kutta(f, y0, t0, tf, h)% f: 微分方程函数，输入为[y, t]，输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Nk1 = h * f(y(:, i), t(i));k2 = h * f(y(:, i) + k1 / 2, t(i) + h / 2);k3 = h * f(y(:, i) + k2 / 2, t(i) + h / 2);k4 = h * f(y(:, i) + k3, t(i) + h);y(:, i + 1) = y(:, i) + (k1 + 2 * k2 + 2 * k3 + k4) / 6;t(i + 1) = t(i) + h;endend```2. 欧拉法（Euler's method）：```matlabfunction [y, t] = euler_method(f, y0, t0, tf, h)% f: 微分方程函数，输入为[y, t]，输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Ny(:, i + 1) = y(:, i) + h * f(y(:, i), t(i));t(i + 1) = t(i) + h;endend```使用这两个函数时，需要定义一个表示微分方程的函数`f`，例如：```matlabfunction dydt = my_ode(y, t)dydt = -y; % 一个简单的一阶线性微分方程：dy/dt = -yend```然后调用相应的求解函数，例如：```matlaby0 = 1; % 初始值t0 = 0; % 初始时间tf = 5; % 结束时间h = 0.1; % 步长[y_rk, t_rk] = runge_kutta(@my_ode, y0, t0, tf, h);[y_euler, t_euler] = euler_method(@my_ode, y0, t0, tf, h);```。
欧拉法是数值分析中常用的一种方法，用于求解常微分方程的数值解。
在MATLAB中，可以通过编写相应的代码来实现欧拉法求解微分方程。
下面我们将通过具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。
我们要了解欧拉法的基本原理。
欧拉法是一种通过迭代逼近微分方程解的方法，它基于微分方程的定义，通过离散化的方法逼近微分方程的解。
其基本思想是利用微分方程的导数定义，将微分方程以差分形式进行逼近。
具体而言，欧拉法通过将微分方程转化为差分方程的形式，然后通过迭代逼近得到微分方程的数值解。
接下来，我们通过一个具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。
假设我们要求解以下的一阶常微分方程：(1) dy/dx = x + y(2) y(0) = 1现在我们来编写MATLAB代码来实现欧拉法求解这个微分方程。
我们需要确定微分方程的迭代步长和迭代范围。
假设我们将x的范围取为0到10，步长为0.1。
接下来，我们可以编写MATLAB代码如下：```matlab欧拉法求解微分方程 dy/dx = x + y定义迭代步长和范围h = 0.1;x = 0:h:10;初始化y值y = zeros(1,length(x));y(1) = 1;使用欧拉法迭代求解for i = 1:(length(x)-1)y(i+1) = y(i) + h * (x(i) + y(i));end绘制图像plot(x,y,'-o');xlabel('x');ylabel('y');title('欧拉法求解微分方程 dy/dx = x + y');```在这段MATLAB代码中，我们首先定义了迭代的步长和范围，并初始化了微分方程的初始值y(0) = 1。
然后通过for循环使用欧拉法进行迭代求解微分方程，最后绘制出了微分方程的数值解的图像。
通过以上的实例讲解，我们可以看到，在MATLAB中使用欧拉法求解微分方程是非常简单而直观的。
解微分方程的欧拉法，龙格-库塔法及其MATLAB简单实例欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前进EULER法、后退EULER法、改进的EULER法。
缺点：欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。
因此欧拉格式一般不用于实际计算。
改进欧拉格式：为提高精度，需要在欧拉格式的基础上进行改进。
采用区间两端的斜率的平均值作为直线方程的斜率。
改进欧拉法的精度为二阶。
算法为：微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。
对于常微分方程：?x∈[a,b]y(a) = y0可以将区间[a,b]分成n段，那么方程在第xi点有y'(xi) = f(xi,y(xi))，再用向前差商近似代替导数则为：在这里，h是步长，即相邻两个结点间的距离。
因此可以根据xi点和yi点的数值计算出yi+1来：i=0,1,2,L这就是向前欧拉格式。
改进的欧拉公式：将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均，即上式便是梯形的欧拉公式。
可见，上式是隐式格式，需要迭代求解。
为了便于求解，使用改进的欧拉公式：数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。
龙格-库塔方法的基本思想：在区间[xn,xn+1]内多取几个点，将他们的斜率加权平均，作为导数的近似。
龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。
令初值问题表述如下。
则，对于该问题的RK4由如下方程给出：其中这样，下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：k1是时间段开始时的斜率；k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值；k4是时间段终点的斜率，其y值用k3决定。
解微分方程的欧拉法，龙格-库塔法及其MATLAB简单实例欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前进EULER法、后退EULER法、改进的EULER法。
缺点：欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。
因此欧拉格式一般不用于实际计算。
改进欧拉格式：为提高精度，需要在欧拉格式的基础上进行改进。
采用区间两端的斜率的平均值作为直线方程的斜率。
改进欧拉法的精度为二阶。
算法为：微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。
对于常微分方程：?x∈[a,b]y(a) = y0可以将区间[a,b]分成n段，那么方程在第xi点有y'(xi) = f(xi,y(xi))，再用向前差商近似代替导数则为：在这里，h是步长，即相邻两个结点间的距离。
因此可以根据xi点和yi点的数值计算出yi+1来：i=0,1,2,L这就是向前欧拉格式。
改进的欧拉公式：将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均，即上式便是梯形的欧拉公式。
可见，上式是隐式格式，需要迭代求解。
为了便于求解，使用改进的欧拉公式：数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。
龙格-库塔方法的基本思想：在区间[xn,xn+1]内多取几个点，将他们的斜率加权平均，作为导数的近似。
龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。
令初值问题表述如下。
则，对于该问题的RK4由如下方程给出：其中这样，下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：k1是时间段开始时的斜率；k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值；k4是时间段终点的斜率，其y值用k3决定。
第8章常微分方程初值问题数值解法8.1 引言8.2 欧拉方法8.3 龙格-库塔方法8.4 单步法的收敛性与稳定性8.5 线性多步法8.1 引 言考虑一阶常微分方程的初值问题00(,),[,],().y f x y x a b y x y '=∈=（1.1） （1.2）如果存在实数 ，使得121212(,)(,).,R f x y f x y L y y y y -≤-∀∈（1.3）则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点)(x y<<<<<+121n n x x x x 上的近似值 .,,,,,121+n n y y y y 相邻两个节点的间距 称为步长.n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ),2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”.即求解过程顺着节点排列的次序一步一步地向前推进.00(,),[,],().y f x y x a b y x y '=∈=描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式.1+n y 一类是计算时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值， 1+n y k 11,,,+--k n n n y y y 称为 步法.k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题.n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式.),(y x f y ='8.2 简单的数值方法8.2.1 欧拉法与后退欧拉法积分曲线上一点 的切线斜率等于函数 的 值.),(y x ),(y x f 如果按函数 在平面上建立一个方向场，那么， ),(y x f xy 积分曲线上每一点的切线方向均与方向场在该点的方向相一致.在 平面上，微分方程的解 称 作它的积分曲线.xy )(x y y =),(y x f y ='基于上述几何解释，从初始点 出发， ),(000y x P 先依切线 在该点的方向推进到 上一点 ，然后再 1x x =1P 从 依切线 的方向推进到 上一点 ，循此前进 1P 2x x =2P 做出一条折线 （图8-1).210P P P 图8-179一般地，设已做出该折线的顶点 ，过 依 切线 的方向再推进到 ，显然两个顶点的坐标有关系 ),(n n n y x P nP ),(111+++n n n y x P 1,+n n P P 11n n n n y y x x ++-=-反解，得).,(1n n n n y x hf y y +=+（2.1）这就是著名的欧拉（Euler ）方法. 欧拉法实际上是对常微分方程中的导数用差商近似，即1()()n n y x y x h+-≈直接得到的.(,)d d n n x y y x =(,)n n f x y 00(,),[,],().y f x y x a b y x y '=∈=1n n y y h +-=()(,())n n n y x f x y x '=111(,)n n n n y y hf x y +=+例欧拉方法13若初值已知，则依公式（2.1）可逐步算出 0y ),,(0001y x hf y y +=),,(1112y x hf y y +=).,(1n n n n y x hf y y +=+（2.1）例1 求解初值问题⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y x y y （2.2）解 欧拉公式的具体形式为).2(1nnn n n y x y h y y -+=+取步长 ，计算结果见表8-1.1.0=h 初值问题（2.2）的解为 ，按这个解析式 子算出的准确值 同近似值 一起列在表8-1中，两者相比较可以看出欧拉方法的精度不高.x y 21+=)(n x y n y 1(,)n n n n y y hf x y +=+12y x=+81()()0.1 1.1000 1.09540.6 1.5090 1.48320.2 1.1918 1.18320.7 1.5803 1.54920.3 1.2774 1.26490.8 1.6498 1.61250.4 1.3582 1.34160.9 1.7178 1.67330.51.43511.41421.01.78481.7321n nn n nn x y y x x y y x -表计算结果对比 还可以通过几何直观来考察欧拉方法的精度.假设 ，即顶点 落在积分曲线 上，那么，按欧拉方法做出的折线 便是 过点 的切线（图8-2）. )(n n x y y =nP )(x y y =1+n n P P )(x y y =nP图8-2从图形上看，这样定出的顶点 显著地偏离了原来 的积分曲线，可见欧拉方法是相当粗糙的.1+n P 为了分析计算公式的精度，通常可用泰勒展开将)(1+n x y在处展开，则有n x)()(1h x y x y n n +=+在 的前提下， )(n n x y y =)())(,(),(n n n n n x y x y x f y x f '==11()n n y x y ++-=称为此方法的局部截断误差.()()n n y x y x h '=++于是可得欧拉法（2.1）的误差 （2.3）),(22n x y h ''≈如果对方程 从 到积分，得 n x 1+n x ),(y x f y =').,(1n n n n y x hf y y +=+（2.1） 2()2n h y ξ''1(,)n n n x x ξ+∈1(,)n n n n y y hf x y +=+(上一值精确相等时，估计下一值误差!) 2()2n h y ξ''1()n y x +=（2.4）1()(,())n nx n xy x f t y t dt++⎰18右端积分用左矩形公式 近似.))(,(n n x y x f h 再以 代替 n y ),(n x y 如果在（2.4）中右端积分用右矩形公式 ))(,(11++n n x y x f h 111(,)n n n n y y hf x y +++=+（2.5）称为后退的欧拉法.代替也得到（2.1）， )(1+n x y 1+n y 局部截断误差也是（2.3）. 近似，则得另一个公式它也可以通过利用均差近似导数 ，即11()()n n n ny x y x x x ++-≈-)(1+'n x y 1(,)n n n n y y hf x y +=+2()2n hy ξ''1()(,())n nx n x y x f t y t dt++⎰111()(,())n n n y x f x y x +++'=欧拉公式是关于 的一个直接的计算公式，这类公式称作是显式的；1+n y 后退欧拉公式的右端含有未知的 ，它是关于 的一个函数方程，这类公式称作是隐式的.1+n y 1+n y 1(,)n n n n y y hf x y +=+111(,)n n n n y y hf x y +++=+隐式方程通常用迭代法求解，而迭代过程的实质是逐步显示化. 设用欧拉公式),()0(1n n n n y x f h y y+=+给出迭代初值 ，用它代入（2.5）式的右端，使之转化 为显式，直接计算得)0(1+n y ),,()0(11)1(1++++=n n n n yx f h y y然后再用 代入（2.5）式，又有)1(1+n y ).,()1(11)2(1++++=n n n n yx f h y y111(,)n n n n y y hf x y +++=+（2.5）21如此反复进行，得).,1,0(),,()(11)1(1=+=++++k yx f h y yk n n n k n （2.6）由于 对 满足李普希茨条件(1.3). ),(y x f y 由(2.6)减(2.5)得(1)11k n n yy +++-=.1)(1++-≤n k n y y hL 由此可知，只要迭代法(2.6)就收敛到解 . 1<hL 1+n y 从积分公式可以看到后退欧拉方法的公式误差与欧拉法是相似的.111(,)n n n n y y hf x y +++=+（2.5） ()1111(,)(,)k n n n n h f x y f x y ++++-121212(,)(,).,Rf x y f x y L y y y y -≤-∀∈8.2.2 梯形方法 若用梯形求积公式近似等式（2.4）右端的积分并分别用 代替 则可得到比欧拉法 精度高的计算公式1,+n n y y ),(),(1+n n x y x y 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） 称为梯形方法. 梯形方法是隐式单步法，可用迭代法求解.⎰+1))(,(n nx x dtt y t f .))(,()()(11⎰++=+n nx xn n dt t y t f x y x y （2.4） 1()(,())n nx n xy x f t y t dt++⎰111(,())[(,)(,)]2n n x n n n n x hf t y t dt f x y f x y +++≈+⎰23⎪⎪⎩⎪⎪⎨⎧+=+);,()0(1n n n ny x f h y y 为了分析迭代过程的收敛性，将（2.7）与（2.8）式相减，得)],(),([2)(11)1(1k nn n n n k n y x f y x f h y y ++++++=（2.8） ).,2,1,0( =k 同后退的欧拉方法一样，仍用欧拉方法提供迭代初值， 则梯形法的迭代公式为111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） (1)()111111[(,)(,)]2k k n n n n n n h y yf x y f x y +++++++-=-如果选取 充分小，使得h 12hL<则当 时有 ， ∞→k 1)(1++→n k n y y 此时迭代过程是收敛的.于是有(1)11k n n y y+++-≤式中 为关于 的李普希茨常数. ),(y x f y L (1)()111111[(,)(,)]2k k n n n n n n h y yf x y f x y +++++++-=-121212(,)(,).,Rf x y f x y L y y y y -≤-∀∈()112k n n hL y y ++-计算方法及MATLAB 实现8.2.3 改进欧拉公式梯形方法虽然提高了精度，但其算法复杂. 在应用迭代公式（2.8）进行实际计算时，每迭代一次，都要重新计算函数 的值.),(y x f 为了控制计算量，通常只迭代一两次就转入下一步的 计算，这就简化了算法.具体地，先用欧拉公式求得一个初步的近似值 , 1n y + 而迭代又要反复进行若干次，计算量很大，而且往往难以预测.称之为预测值，)],(),([2)(11)1(1k n n n n n k n y x f y x f h y y ++++++=计算方法及MATLAB 实现这样建立的预测-校正系统通常称为改进的欧拉公式： 预测值 的精度可能很差，再用梯形公式（2.7）将它校正一次，即按（2.8）式迭代一次得 ，这个结果称校正值. 1+n y 1+n y 预测 ),,(1n n n n y x f h y y +=+校正)].,(),([2111+++++=n n n n n n y x f y x f hy y （2.9）也可以表为下列平均化形式),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(1y y y +=111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） )],(),([2)(11)1(1k n n n n n k n y x f y x f h y y ++++++=(参阅教材392页!)计算方法及MATLAB 实现例2 用改进的欧拉方法求解初值问题（2.2）.解 这里 改进的欧拉公式为⎪⎪⎪⎩⎪⎪⎪⎨⎧-+=),2(nnn n p y x y h y y ⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y x y y （2.2）2(,)(),xf x y y y=-),2(1p n p n c y x y h y y +-+=).(211c p n y y y +=+),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(211c p n y y y +=+计算方法及MATLAB 实现282(,)xf x y y y=-),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(211c p n y y y +=+取 ，计算结果见表8-2.1.0=h 同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.()()0.11.0959 1.09540.6 1.4860 1.48320.2 1.1841 1.18320.7 1.5525 1.54920.3 1.2662 1.26490.8 1.6153 1.61250.4 1.3434 1.34160.9 1.6782 1.67330.51.41641.41421.01.73791.7321n n n n n n x y y x x y y x -表8212y x=+9.2.4 单步法的局部截断误差与阶 初值问题(1.1)，(1.2)的单步法可用一般形式表示为),,,,(11h y y x h y y n n n n n +++=ϕ（2.10）其中多元函数 与 有关，ϕ),(y x f 当 含有 时，方法是隐式的，若不含 则为显式方法，ϕ1+n y 1+n y ),,,(1h y x h y y n n n n ϕ+=+（2.11） 称为增量函数， ),,(h y x ϕ).,(),,(y x f h y x =ϕ所以显式单步法可表示为 例如对欧拉法（2.1）有 它的局部截断误差已由(2.3)给出.⎩⎨⎧=='.)(),,(00y x y y x f y （1.1） （1.2） ).,(1n n n n y x hf y y +=+（2.1） )(2)(211n n n y h y x y ξ''=-++（2.3） ),(22n x y h ''≈对一般显式单步法则可如下定义. 定义1 设 是初值问题(1.1)，(1.2)的准确解， 称 )(x y y =)),(,()()(11h x y x h x y x y T n n n n n ϕ--=++（2.12） 为显式单步法（2.11）的局部截断误差. 之所以称为局部的，是假设在 前各步没有误差.1+n T n x 当 时，计算一步，则有 )(n n x y y =11()n n y x y ++-=1()()(,(),)n n n n y x y x h x y x h ϕ+=--⎩⎨⎧=='.)(),,(00y x y y x f y （1.1）（1.2） ),,,(1h y x h y y n n n n ϕ+=+（2.11）1()[(,,)]n n n n y x y h x y h ϕ+-+1.n T +=1(,)n n n n y y hf x y +=+).,(),,(y x f h y x =ϕ在前一步精确的情况下用显式单步公式计算产生的公式误差. 根据定义，欧拉法的局部截断误差 ))(,()()(11n n n n n x y x hf x y x y T --=++即为（2.3）的结果. 这里 称为局部截断误差主项. )(22n x y h '')(2h O T =局部截断误差可理解为用显式单步方法计算一步的误差，即 ()()n n y x h y x =+-22()()2n h y x o h ''=+显然 ),,,(1h y x h y y n n n n ϕ+=+（2.11）)(2)(211n n n y h y x y ξ''=-++（2.3） ),(22n x y h ''≈()n hy x '-(,)y f x y '=1()()n n y x y x h +=+(泰勒公式) (小o 高阶，大O 同阶)计算方法及MATLAB 实现 定义2 设 是初值问题(1.1)，(1.2)的准确解， 若存在最大整数使显式单步法(2.11)的局部截断误差满足 )(x y p ),(),,()()(11++=--+=p n h O h y x h x y h x y T ϕ（2.13） 则称显式单步法具有 阶精度. p 若将（2.13）展开式写成),())(,(211++++=p p n n n h O h x y x T ψ则 称为局部截断误差主项. 1))(,(+p n n hx y x ψ 以上定义对隐式单步法（2.10）也是适用的.⎩⎨⎧=='.)(),,(00y x y y x f y （1.1） （1.2） ),,,,(11h y y x h y y n n n n n +++=ϕ（2.10） 1(,,)n n n n y y h x y h ϕ+=+(小o 高阶，大O 同阶)对后退欧拉法（2.5）其局部截断误差为))(,()()(1111++++--=n n n n n x y x hf x y x y T 这里 ，是1阶方法，局部截断误差主项为 . 1=p )(22n x y h ''-23()()()2n n h hy x y x O h '''=++).()(232h O x y h n +''-=),,(111++++=n n n n y x hf y y （2.5）2[()()()]n n h y x hy x O h '''-++对梯形法（2.7）有)]()([2)()(111+++'+'--=n n n n n x y x y h x y x y T 所以梯形方法是二阶的，其局部误差主项为 )(123n x y h '''-23()()()23!n n n h h hy x y x y x ''''''=++).()(1243h O x y h n +'''-=)],,(),([2111+++++=n n n n n n y x f y x f h y y （2.7） 24[()()()()]()22n n n n h h y x y x hy x y x O h '''''''-++++3364h h -8.3龙格-库塔方法8.3.1 显式龙格-库塔法的一般形式 上节给出了显式单步法的表达式),,,(1h y x h y y n n n n ϕ+=+其局部截断误差为),(),,()()(11++=--+=p n h O h y x h x y h x y T ϕ对欧拉法 ，即方法为阶， )(21h O T n =+1=p ))].,(,(),([21n n n n n n n n y x hf y h x f y x f h y y ++++=+（3.1）若用改进欧拉法，它可表为此时增量函数 ))].,(,(),([21),,(n n n n n n n n y x hf y h x f y x f h y x +++=ϕ（3.2）与欧拉法的 相比，增加了计算一个 右函数 的值，可望 .),(),,(n n n n y x f h y x =ϕf 2=p 若要使得到的公式阶数 更大， 就必须包含更多的值.p ϕf ,))(,()()(11⎰+=-+n n x x n n dx x y x f x y x y （3.3）从方程 等价的积分形式（2.4），即),(y x f y ='若要使公式阶数提高，就必须使右端积分的数值求积公式 精度提高，必然要增加求积节点.为此可将（3.3）的右端用求积公式表示为线性组合1(,())n n x x f x y x dx +≈⎰点数 越多，精度越高， r 上式右端相当于增量函数 ， ),,(h y x ϕ为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为),,,(1h y x h y y n n n n ϕ+=+（3.4）其中1(,())r i n i n i i h c f x h y x h λλ=++∑,))(,()()(11⎰+=-+n n x x n n dx x y x f x y x y （3.3）(可变步长因子!) (组合系数!),),,(1∑==r i i i n n K c h y x ϕ),,(1n n y x f K =),(11∑-=++=i j j ij n i n i K h y h x f K μλ（3.5）,,,2r i =这里 均为常数. ij i i c μλ,, (3.4)和(3.5)称为 级显式龙格-库塔（Runge-Kutta ）法， r 简称R-K 方法.当 时，就是欧拉法，),(),,(,1n n n n y x f h y x r ==ϕ此时方法的阶为. 1=p 当 时，改进欧拉法(3.1)，(3.2)也是其中的一种，2=r ),,,(1h y x h y y n n n n ϕ+=+（3.4） ),,(1h y x h y y n n n n ϕ+=+))].,(,(),([21),,(n n n n n n n n y x hf y h x f y x f h y x +++=ϕ).,(1n n n n y x hf y y+=+1(,,)n n n n y y h x y h ϕ+=+(步长调节因子) (组合系数) (斜率组合系数)8.3.2 二阶显式R-K 方法 1(,,)rn n i ii x y h c K ϕ==∑1(,,)n n n n y y h x y h ϕ+=+ 对 的R-K 方法，计算公式如下2=r 11122122211()(,)(,)n n n n n n y y h c K c K K f x y K f x h y hK λμ+=++⎧⎪=⎨⎪=++⎩（3.6）这里 均为待定常数.21221,,,μλc c ),(11∑-=++=i j j ij n i n i K h y h x f K μλ若取 得计算公式,2/1,1,021221====μλc c 12n n y y hK +=+称为中点公式，相当于数值积分的中矩形公式.上式也可表示为1(,(,))22n n n n n n h hy y hf x y f x y +=+++1(,)n n K f x y =（3.10）21(,)22n n h hK f x y K =++)).,(2,2(1n n n n n n y x f hy h x hf y y +++=+11122122211()(,)(,)n n n n n n y y h c K c K K f x y K f x h y hK λμ+=++⎧⎪=⎨⎪=++⎩继续上述过程，经过较复杂的数学演算，可以导出各 种四阶龙格-库塔公式，下列经典公式是其中常用的一个：可以证明其截断误差为 . )(5h O 四阶龙格-库塔方法的每一步需要计算四次函数值 ，f ⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+).,(),2,2(),2,2(),,(),22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n n n n n n n （3.13） (凸线性组合!) (中点公式!) (欧拉公式!)例 设取步长 ，从 到 用四阶龙格-库塔方法求解 初值问题 ⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y xy y 解 这里 ，经典的四阶龙格-库塔公式为 yx y y x f 2),(-=),22(643211K K K K hy y n n ++++=+,2nx y K -=2.0=h 0=x 1=x 112341213243(22),6(,),(,),22(,),22(,).n n n nn n n n n n h y y K K K K K f x y h h K f x y K h hK f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪⎪=++⎪⎩计算方法及MATLAB 实现,222223K h y h x K h y K n n n ++-+=,222112K h y h x K hy K n n n ++-+=.)(2334hK y h x hK y K n n n ++-+= 表8-3列出了计算结果，同时了出了相应的精确解. 龙格-库塔方法的精度较高.112341213243(22)6(,),(,),22(,),22(,).n n n nn n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪⎪=++⎪⎩yxy y x f 2),(-=,21nnn y x y K -=83()0.2 1.1832 1.18320.4 1.3417 1.34160.6 1.4833 1.48320.8 1.6125 1.61251.01.73211.7321n n n x y y x 表计算结果。
