计算方法龙格库塔方法概要
龙格-库塔(Runge-Kutta)法
龙格-库塔（Runge-Kutta）法 1.1 龙格-库塔(Runge-Kutta)法的基本思想
Euler公式可改写成
yi1 yi hK1 K1 f ( xi , yi )
则yi+1的表达式y(xi+1)与的Taylor展开式的前两项 完全相同,即局部截断误差为 O(h 2 ) 。
为了进一步提高精度，设除 xi p 外再增加一点
xiq xi qh ( p q 1)
并用三个点 xi ,xi p , xiq 的斜率k1,k2,k3加权平均
得出平均斜率k*的近似值，这时计算格式具有形式:
yi1 yi h(1 )k1 k2 k3
k1 f (xi , yi ) k2 f (xi ph, yi phk1 )
格式。
若取 1 0 ,则 2 法的计算公式为
1,
p
1 2
，此时二阶龙格-库塔
ky1i
1
f
yi hk2 ( xi , yi )
k
2
h
f
(
x
i
1
,
yi
2
2 k1 )
i 0,1,2, n 1
此计算公式称为变形的二阶龙格—库塔法。式中
x 1 i 2
为区间
xi , xi1
的中点。
1.3 三阶龙格-库塔法
拉法,将 xi p 视为 xi1，即可得
k2 f (xi ph, yi phk1 ) 对常微分方程初值问题(7.1)式的解 y=y(x),根据微 分中值定理，存在点 (xi , xi1 ) ，使得
也即
y(xi1 ) y(xi ) y( )( xi1 xi )
y( xi1 ) y( xi ) hK
龙格库塔法
一、高阶泰勒法
假设初值问题
龙格—库塔法 龙格 库塔法
dy = f (t , y ) dt y (a) = α 的解y (t)及f (t , y )足够光滑.
将y (ti +1 )在ti处作n阶泰勒展开, 得
a≤t ≤b
(1)
y′′(ti ) 2 y ( n ) (ti ) n y ( n +1) (ξ i ) n +1 y (ti +1 ) = y (ti ) + y′(ti )h + h +L+ h + h n! 2! (n + 1)! 其中, ti < ξ i < ti +1.
2
i
i
1
3
i
i
2
4
i
i
3
i +1
i
6123 Nhomakorabea4
作业 教材P198 习题3
(2)
(3)
首先将y (ti +1 )在ti处展成幂级数 h2 y (ti +1 ) = y (ti ) + hy′(ti ) + y′′(ti ) + O(h 3 ) 2 将 y′(t ) = f (t , y (t )) y′′(t ) = f t′(t , y (t )) + f y (t , y (t )) f (t , y (t )) 代入上式, 得 h2 y (ti +1 ) = y (ti ) + hf + ( f t + ff y ) + O(h 3 ) (3) 2 其中f , f t , f y′分别表示相应函数在点(ti , y (ti ))处的函数值.
龙格库塔
数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
经典四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。
令初值问题表述如下。
则，对于该问题的RK4由如下方程给出：其中这样，下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：∙k1是时间段开始时的斜率；∙k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点t n + h/2的值；∙k3也是中点的斜率，但是这次采用斜率k2决定y值；∙k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值：RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
显式龙格库塔法显示龙格－库塔法是上述RK4法的一个推广。
它由下式给出其中如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(h p+1)时的条件。
这些可以从舍入误差本身的定义中导出。
例如，一个2阶精度的2段方法要求b1 + b2 = 1, b2c2 = 1/2, 以及b2a21 = 1/2。
在Matlab下输入：edit，然后将下面两行百分号之间的内容，复制进去，保存%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%function dxdt=ode_Miss_ghost(t,x)%分别用x(1),x(2),x(3),x(4)代替N1,P1,N2,P2N1=x(1);P1=x(2);N2=x(3);P2=x(4);K=2;tau_c=3e-9;tan_p=6e-12;beta =5e-5;delta=0.692;eta =0.0001;fm =8e6;Ith =26e-3;Ib =1.5*Ith;Im =0.3*Ith;I1=Ib+Im*sin(2*pi*fm*t)+K*P2;I2=Ib+Im*sin(2*pi*fm*t)+K*P1;dxdt=[(I1/Ith-N1-(N1-delta)/(1-delta)*P1)/tau_e;((N1-delta)/(1-delta)*(1-eta*P1)*P1-P1+beta*N1)/tau_p;(I2/Ith-N2-(N2-delta)/(1-delta)*P2)/tau_e;((N2-delta)/(1-delta)*(1-eta*P2)*P2-P2+beta*N2)/tau_p;]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在Matlab下面输入：t_start=0;t_end=2e-9;y0=[1e-3;1e-4;0;0]; %初值[x,y]=ode15s('ode_Miss_ghost',[0,t_end],y0);plot(x,y);legend('N1','P1','N2','P2');xlabel('x');Mathlab定步长龙格－库塔法[345阶]、自适应步长rkf45经常看到很多朋友问定步长的龙格库塔法设置问题，下面吧定步长三阶、四阶、五阶龙格库塔程序贴出来，有需要的可以看看ODE3 三阶龙格－库塔法CODE:function Y = ode3(odefun,tspan,y0,varargin)%ODE3 Solve differential equations with a non-adaptive method of order 3.% Y = ODE3(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates% the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE3(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN% but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Bogacki-Shampine Runge-Kutta method of order 3.%% Example% tspan = 0:0.1:20;% y = ode3(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1, % and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,3);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});F(:,3) = feval(odefun,ti+0.75*hi,yi+0.75*hi*F(:,2),varargin{:});Y(:,i) = yi + (hi/9)*(2*F(:,1) + 3*F(:,2) + 4*F(:,3));endY = Y.';ODE4 四阶龙格－库塔法CODE:function Y = ode4(odefun,tspan,y0,varargin)%ODE4 Solve differential equations with a non-adaptive method of order 4.% Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the classical Runge-Kutta method of order 4.%% Example% tspan = 0:0.1:20;% y = ode4(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,4);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:}); F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:}); F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));endY = Y.';定步长RK4(自编):CODE:function Y=RungeKutta4(f,xn,y0)% xn=0:.1:1;% y0=1;% y_n=[];% f=@(X1,Y1) Y1-2*X1/Y1;y_n=[];h=diff(xn(1:2));for i=1:length(xn)-1K1=f(xn(i),y0);K2=f(xn(i)+h/2,y0+h*K1/2);K3=f(xn(i)+h/2,y0+h*K2/2);K4=f(xn(i)+h,y0+h*K3);y_n=[y_n;y0+h/6*(K1+2*K2+2*K3+K4)];y0=y_n(end);endY=y_n;ODE5 五阶龙格－库塔法CODE:function Y = ode5(odefun,tspan,y0,varargin)%ODE5 Solve differential equations with a non-adaptive method of order 5.% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Dormand-Prince method of order 5 in a general % framework of explicit Runge-Kutta methods.%% Example% tspan = 0:0.1:20;% y = ode5(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);% Method coefficients -- Butcher's tableau%% C | A% --+---% | BC = [1/5; 3/10; 4/5; 8/9; 1];A = [ 1/5, 0, 0, 0, 03/40, 9/40, 0, 0, 044/45 -56/15, 32/9, 0, 019372/6561, -25360/2187, 64448/6561, -212/729, 0 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]; % More convenient storageA = A.';B = B(:);nstages = length(B);F = zeros(neq,nstages);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);% General explicit Runge-Kutta frameworkF(:,1) = feval(odefun,ti,yi,varargin{:});for stage = 2:nstageststage = ti + C(stage-1)*hi;ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));F(:,stage) = feval(odefun,tstage,ystage,varargin{:});endY(:,i) = yi + F*(hi*B);endY = Y.';-------------------------------------------------------------------------------------------------------------------自适应步长RKF45(相当于ode45)ODE45 是4阶方法提供候选解，5阶方法控制误差。
matlab龙格库塔法程序,给出实例
一、介绍龙格库塔法龙格库塔法（Runge-Kutta method）是一种数值计算方法，用于求解常微分方程的数值解。
它通过多步迭代的方式逼近微分方程的解，并且具有较高的精度和稳定性。
二、龙格库塔法的原理龙格库塔法采用迭代的方式来逼近微分方程的解。
在每一步迭代中，计算出当前时刻的斜率，然后根据这个斜率来求解下一个时刻的值。
通过多步迭代，可以得到微分方程的数值解。
三、龙格库塔法的公式龙格库塔法可以表示为以下形式：k1 = f(tn, yn)k2 = f(tn + h/2, yn + h/2 * k1)k3 = f(tn + h/2, yn + h/2 * k2)k4 = f(tn + h, yn + h * k3)yn+1 = yn + h/6 * (k1 + 2k2 + 2k3 + k4)其中，k1、k2、k3、k4为斜率，h为步长，tn为当前时刻，yn为当前时刻的解，yn+1为下一个时刻的解。
四、使用matlab实现龙格库塔法在MATLAB中，可以通过编写函数来实现龙格库塔法。
下面是一个用MATLAB实现龙格库塔法的简单例子：```matlabfunction [t, y] = runge_kutta(f, tspan, y0, h)t0 = tspan(1);tf = tspan(2);t = t0:h:tf;n = length(t);y = zeros(1, n);y(1) = y0;for i = 1:n-1k1 = f(t(i), y(i));k2 = f(t(i) + h/2, y(i) + h/2 * k1);k3 = f(t(i) + h/2, y(i) + h/2 * k2);k4 = f(t(i) + h, y(i) + h * k3);y(i+1) = y(i) + h/6 * (k1 + 2*k2 + 2*k3 + k4);endend```以上就是一个简单的MATLAB函数，可以利用该函数求解给定的微分方程。
四阶龙格—库塔法的原理及其应用
《四阶龙格—库塔法的原理及其应用》
龙格—库塔法(又称龙格库塔法)是由一系列有限的、独立的可能解组成的无穷序列,这些解中每个都与原来的数列相差一个常数。
它是20世纪30年代由匈牙利著名数学家龙格和库塔提出的,故得此名。
1.它的基本思想是:在n 阶方阵M 上定义一个函数,使得当n 趋于无穷时,它在m 中所表示的数值为M 的某种特征值,从而构造出一族具有某种特性的可计算函数f (x)= Mx+ C (其中C 为任意正整数)。
例如,若f (x)=(a-1) x+ C,则称之为(a-1) x 的龙格—库塔法。
2.它的应用很广泛,可以求解各类问题,且能将大量的未知数变换成少数几个已知数,因此它是近似计算的一种重要工具。
3.
它的优点主要有:(1)可以将多项式或不等式化成比较简单的形式;(2)对于同一问题可以用不同的方法来解决,并取得同样的结果;(3)适合处理高次多项式或者不等式,尤其适合处理多元函数的二次型。
第三部分龙格-库塔方法
内江师范学院数学与信息科学学院 吴开腾 制作
于是有
其中
y ( xn +1 ) − y ( xn ) = y '(ξ ), ξ ∈ ( xn , xn +1 ) h y ( xn +1 ) = y ( xn ) + hf (ξ , y (ξ ))
k * = f (ξ , y (ξ )) 称作区间 [ xn , xn +1 ] 上的平均斜率。 上的平均斜率 平均斜率。 问题：计算近似值y ( xn +1 ) 的关键是如何选择算法确定平均斜率 k *
（15）
f ( xn +1 , yn + h ( − k1 + 2 k 2 ))
内江师范学院数学与信息科学学院 吴开腾 制作
注释1 可以用Taylor展示证明格式（14） 注释1：可以用Taylor展示证明格式（14）具有三阶精 展示证明格式
度，并且还可以用类似的方法得到四阶及其以上的更高 阶精度的Runge-Kutta格式 阶精度的Runge-Kutta格式。 Runge 格式。
内江师范学院数学与信息科学学院 吴开腾 制作
h yn + （ k1 + 2 k 2 + 2 k3 +k 4） 6 f ( xn , y n ) h f ( x 1 , yn + k1 ) n+ 2 2 h f ( x 1 , yn + k 2 ) n+ 2 2 f ( xn +1 , yn + hk3 ) （16）
四阶龙格－ 四阶龙格－库塔格式计算结果
xn
0.1 0.2 0.3 0.4 0.5
yn
欧拉格式计算结果 xn yn y ( xn )
龙格库塔法介绍
yn
hf
(xn, yn ))]，
(x, y,h) 1[ f (x, y) f (x h, y hf (x, y))],
2
|
( x,
y1,
h)
(x,
y2 ,
h)
|
[L
2
L 2
(1
hL)]
|
y1
y2
|,
L
L(1
h0L)，h 2
h0.
类似地,不难验证其他龙格 库塔方法的收敛性.
这里c1,c2,c3,2,3, 21, 31, 32均为待定参数.
Tn1 y(xn1) yn1 O(h4 )
(3.11)
c1 c2 c3 1
2
21
3 31 32
c22
c33
1 2
cc232223c2332
将步长折半，从xn用两步求xn1处的近似值，则有
y(xn1)
h
yn21
2c
h 2
5
.
从而
h
y ( xn 1) y ( xn 1)
yn21 ynh1
1， 16
得到事后估计式：
y ( xn 1)
h
yn21
1 15
(
h
yn21
ynh1).
通过检查步长折半前后计算结果的偏差，
y(x) (x, y(x),0) 0 p 1 单步法(4.1)收敛. 定义4 若单步法(4.1)增量函数(x, y,h)是否满足
龙格-库塔法
四阶龙格－库塔法求解常微分方程的初值问题1.算法原理对于一阶常微分方程组的初值问题⎪⎪⎪⎩⎪⎪⎪⎨⎧=⋯⋯==⋯⋯=⋯⋯⋯⋯=⋯⋯=0020********'212'2211'1)(,,)(,)())(,),(),(,()())(,),(),(,()())(,),(),(,()(n n n n n n n y x y y x y y x y x y x y x y x f x y x y x y x y x f x y x y x y x y x f x y ， 其中b x a ≤≤。
若记Tn Tn Tn y x f y x f y x f y x f y y y y x y x y x y y x y )),(,),,(),,((),(),,,())(),(),(()(2102010021⋯⋯=⋯⋯=⋯⋯=，，则可将微分方程组写成向量形式⎩⎨⎧=≤≤=0')()),(,()(y a y b x a x y x f x y微分方程组初值问题在形式上和单个微分方程处置问题完全相同，只是数量函数在此变成了向量函数。
因此建立的单个一阶微分方程初值问题的数值解法，可以完全平移到求解一阶微分方程组的初值问题中，只不过是将单个方程中的函数转向向量函数即可。
标准4阶R-K 法的向量形式如下：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()21,2()21,2(),()22(61342312143211K y h x hf K K y h x hf K K y h x hf K y x hf K K K K K y y n n n n n n n n n n 其分量形式为n j K y K y K y h x hf K K y K y K y h x hf K K y K y K y h x hf K y y y x hf K K K K K y y n ni i i i j j n nii i i j j n nii i i j j ni i i i j j j j j j i j i j ,,2,1).,,,;(),2,2,2;2(),2,2,2;2(),,,,;(),22(6132321314222212131212111221143211,1,⋯⋯=⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+⋯⋯+++=+⋯⋯+++=+⋯⋯+++=⋯⋯=++++=++，，2.程序框图3.源代码%该函数为四阶龙格-库塔法function [x,y]=method(df,xspan,y0,h)%df为常微分方程，xspan为取值区间，y0为初值向量，h为步长x=xspan(1):h:xspan(2);m=length(y0);n=length(x);y=zeros(m,n);y(:,1)=y0(:);for i=1:n-1k1=feval(df,x(i),y(:,i));k2=feval(df,x(i)+h/2,y(:,i)+h*k1/2);k3=feval(df,x(i)+h/2,y(:,i)+h*k2/2);k4=feval(df,x(i)+h,y(:,i)+h*k3);y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;end%习题9.2clear;xspan=[0,1];%取值区间h=0.05;%步长y0=[-1,3,2];%初值df=@(x,y)[y(2);y(3);y(3)+y(2)-y(1)+2*x-3];[xt,y]=method(df,xspan,y0,h)syms t;yp=t*exp(t)+2*t-1;%微分方程的解析解yp1=xt.*exp(xt)+2*xt-1%计算区间内取值点上的精确解[xt',y(1,:)',yp1']%y(1,:)为数值解，yp1为精确解ezplot(yp,[0,1]);%画出解析解的图像hold on;plot(xt,y(1,:),'r');%画出数值解的图像4.计算结果。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
同理，改进Euler公式可改写成
1 1 yi 1 yi 2 K1 2 K 2 K1 hf ( xi , yi ) K hf ( x h, y K ) i i 1 2
局部截断误差为O(h3)
上述两组公式在形式上共同点:都是用f(x,y)在某 些点上值的线性组合得出y(xi+1)的近似值yi+1, 且增 加计算的次数f(x,y)的次数,可提高截断误差的阶。如 欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进 欧拉法需计算两次f(x,y)的值，为二阶方法。
K hy(i ) hf i , y (i )
yi 1 yi K
K可以认为是y y( x)在区间 [ xi , xi 1 ]上的平均斜率
y
只要使用适当的方法求 出y ( x)在区 间[ xi , xi 1 ]上平均斜率的近似值K
y y( x)
K
就可得到相应的Runge-Kutta方法
2018/11/9 5
于是可考虑用函数f(x,y)在若干点上的函数值的 线性组合来构造近似公式，构造时要求近似公式在 (xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式 的前面几项重合，从而使近似公式达到所需要的阶数。 既避免求高阶导数,又提高了计算方法精度的阶数。 或者说,在[xi,xi+1]这一步内多计算几个点的斜率值， 然后将其进行加权平均作为平均斜率，则可构造出更 高精度的计算格式，这就是龙格—库塔（RungeKutta）法的基本思想。
称为P阶龙格－库塔方法。 其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作 Tailor展开，通过相同项的系数确定参数。
2018/11/9 7
Runge-Kutta方法的推导思想 对于常微分方程的初值问题
y f ( x , y ) y( a ) y0 a xb
的解y=y(x)，在区间[xi, xi+1]上使用微分中值定理，有
y( xi 1 ) y( xi ) y ( i )( xi 1 xi )
其中 i ( xi , xi 1 )
即
2018/11/9
y( xi 1 ) y( xi ) hy( i )
8
引入记号
y ( xi 1) y ( xi ) K
2018/11/9
xi
xi 1
x
9
如果以y( x)在xi处的斜率作为y( x)在[ xi , xi 1 ]上的平均斜率K
即
K hy( xi ) hf [ xi , y ( xi )]
hf ( xi , yi )
如下图
y
则上式化为
yi 1 yi hf ( xi , yi )
2018/11/9
6
一般龙格－库塔方法的形式为
yi 1 yi c1K1 c2 K 2 c p K p K1 hf ( xi , yi ) K 2 hf ( xi a2h, yi b21K1 ) • • • • • • • • • • • K p hf ( xi a p h, yi b p1K1 b p , p 1K p 1 )
h2 y( xi 1 ) y( xi ) hy( xi ) y( xi ) O(h3 ) 2 yi hf ( xi , yi )
h2 ( xi , yi ) O(h 3 ) f x ( xi , yi ) f ( xi , yi ) f y 2
类似地，若取前P+1项作为y(xi+1)的近似值，便得到
2018/11/9 2
显然p=1时,
y i+1=y i+hf(xi,y i)
它即为我们熟悉的Euler方法。 当p≥2时,要利用泰勒方法就需要计算f(x,y)的高 阶微商。这个计算量是很大的，尤其当f(x,y)较复 杂时，其高阶导数会很复杂。因此，利用泰勒公 式构造高阶公式是不实用的。但是泰勒级数展开 法的基本思想是许多数值方法的基础。 R-K方法不是直接使用Taylor级数,而是利用它的思想
§
9.4 龙格－库塔方法
得到高精度方法的一个直接想法是利用Taylor展开
假设式 y' =f(x,y) (a≤x≤b) 中的 f(x,y) 充分光滑,将y(xi+1)在x i点作Taylor展开，若 取右端不同的有限项作为y(xi+1)的近似值，就可得到 计算y(xi+1)的各种不同截断误差的数值公式。
例如：取前两项可得到
y( xi 1 ) y( xi ) hy( xi ) O(h 2 )
y( xi ) hf ( xi , y( xi )) O(h 2 ) yi hf ( xi , yi ) O(h 2 )
2018/11/9 1
若取前三项，可得到截断误 P ( P) yi hyi yi yi 2! P!
P阶泰勒方法
其中 yi f , yi f ( xi , yi ) x f x ff y )x f xx 2 f xy f f yy f 2 f x f y ( fy )2 f yi ( f x ff y
2018/11/9 3
Runge-Kutta 方法是一种高精度的单步法,简称R-K法
9.4.1 龙格-库塔(R-K)法的基本思想
Euler公式可改写成
yi 1 yi K K hf ( xi , yi )
则 yi+1 的表达式与 y(xi+1) 的 Taylor 展开式的前两项 完全相同,即局部截断误差为O(h2)。
y y( x)
K
K
即Euler方法
xi
xi 1
x
Euler方法也称为一阶Runge-Kutta方法
2018/11/9 10
9.4.2