龙格库塔积分算法
龙格-库塔方法基本原理3
2020/4/7
15
令 y(xi1) yi1 对应项的系数相等，得到
c1 c2 1 ,
a2c2
1 2
,
b21c2
1 2
这里有 4 个未知 数，3 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶 龙格 - 库塔格式。
2020/4/7
称为P阶龙格－库塔方法。
其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作 Tailor展开，通过相同项的系数确定参数。
2020/4/7
9
Runge-Kutta方法的推导思想
对于常微分方程的初值问题
y f (x, y) a x b
y(a)
y0
的解y=y(x)，在区间[xi, xi+1]上使用微分中值定理，有
hf (xi 2
,
yi
1 2
K1 )
K3 hf (xi h, yi K1 2K 2 )
2020/4/7
22
四阶(经典)龙格—库塔法
如果需要再提高精度，用类似上述的处理 方法，只需在区间[xi,xi+1]上用四个点处的斜 率加权平均作为平均斜率K*的近似值，构成一 系列四阶龙格—库塔公式。具有四阶精度，即 局部截断误差是O(h5)。
f x(xi , yi ) f (xi , yi ) f y (xi , yi ) O(h3 )
类似地，若取前P+1项作为y(xi+1)的近似值，便得到
yi1
yi
hyi
h2 2!
yi
hP P!
yi(P)
P阶泰勒方法
其中 yi f , yi f (xi , yi )x f x ff y
龙格库塔
数值分析中，龙格－库塔法（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阶方法控制误差。
第三部分龙格-库塔方法
内江师范学院数学与信息科学学院 吴开腾 制作
于是有
其中
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 )
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。
如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。
一阶常微分方程可以写作：y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理，存在0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th))这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。
利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：K1＝f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6);所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。
仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。
想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。
编写的定步长的龙格库塔计算函数：function [x,y]=runge_kutta1(ufunc,y0,h,a,b,Vg)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点,发酵罐体积（参数形式参考了ode45函数）n=floor((b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii),Vg);k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2,Vg);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2,Vg);k4=ufunc(x(ii)+h,y(:,ii)+h*k3,Vg);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end调用的子函数以及其调用语句：function dy=test_fun(x,y)dy = zeros(3,1);%初始化列向量dy(1) = y(2) * y(3);dy(2) = -y(1) + y(3);dy(3) = -0.51 * y(1) * y(2);对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：[T,F] = ode45(@test_fun,[0 15],[1 1 3]);subplot(121)plot(T,F)%Matlab自带的ode45函数效果title('ode45函数效果')[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数subplot(122)plot(T1,F1)%自编的龙格库塔函数效果title('自编的龙格库塔函数')运行结果如下：。
龙格库塔公式
龙格-库塔（Runge-Kutta）公式是一种常用于数值解求微分方程的方法，其中包括了多种不同类型的龙格-库塔公式.
其中最常见的是四阶龙格-库塔公式（Runge-Kutta 4th order method），这种方法通过在当前时刻的函数值和当前时刻的导函数的近似值计算下一时刻的函数值。
它的计算公式如下：
k1 = hf(xn, yn)
k2 = hf(xn + h/2, yn + k1/2)
k3 = hf(xn + h/2, yn + k2/2)
k4 = hf(xn + h, yn + k3)
yn+1 = yn + (k1 + 2k2 + 2k3 + k4)/6
xn+1 = xn + h
其中k1,k2,k3,k4 是近似值，f(x,y)是微分方程，xn, yn 是当前时刻的函数值，xn+1, yn+1 是下一时刻的函数值, h是时间步长.
这种方法具有较高的精度,并且易于实现,但在求解某些类型的方程时,会存在稳定性问题.
除此之外,还有其他类型的龙格-库塔公式，比如二阶、三阶等等,它们也有各自的适用条件.。
龙格-库塔方法-文档资料
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法：
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法：N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时，微分方程是不稳定的； 而 f y 0 时，微分方程是稳定的。
上面讨论的稳定性，与数值方法和方程中 f 有关
•21
实验方程： y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程，
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的：h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例，来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n，
其中 y n 为精确值， n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值，使得公式具有 2阶精度！！
•3
由泰勒公式展开，要使公式具有 2 阶精度，只需
数值分析9-3 龙格-库塔方法
最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ：
y i 1 y i h ( K1 2K 2 2K 3 K 4 ) 6 K1 K2 K3 K4 f ( xi , yi )
h f ( xi h , y K1 ) i 2 2 h f ( xi h , y K2 ) i 2 2
一、再看Taylor方法
h y ( x h) y ( x) hy' ( x) y" ( x) 2 h3 h 4 ( 4) y ( x) y ( x) 6 24 一般可取公式为如下形式
2
y n 1
h h ( p) y n hy' n y"n yn 2! p!
yi 1 K1 K2 yi h [1 K 1 2 K 2 ] f ( x i , yi ) f ( xi ph, yi phK 1 )
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O( h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
考察改进的欧拉法，可以将其改写为： 斜率 一定取K1 K2 的平均值吗？
y i 1 K1 K2
1 1 yi h K 1 K 2 2 2 f ( xi , yi ) f ( x i h, y i hK 1 )
步长一定是一个h 吗？
将改进欧拉法推广为：
龙格库塔方法理解及应用
龙格库塔方法理解及应用龙格库塔方法是一种常用的数值解微分方程的方法，也是许多科学、工程和经济领域中常用的算法之一。
本文将介绍龙格库塔方法的原理及其应用。
一、龙格库塔方法原理在计算微分方程时，往往需要对方程进行离散化，采用数值方法处理。
龙格库塔方法（Runge-Kutta method）就是一种离散化的数值方法，其原理可以概括为：通过相应的递推公式，将微分方程在离散时间点上进行逼近，从而得到近似的解。
具体来说，假设要求解如下形式的一阶常微分方程：$$ y'=f(t,y) $$其中，$f(t,y)$是已知的函数，$y(t)$是未知函数，并且已知初值$y(t_0)=y_0$。
为了离散化这个方程，我们可以采用以下的递推公式：$$ \begin{aligned} y_1 &=y_0 + h\varphi_1 \\ y_2 &=y_0 +h\varphi_2 \\ \cdots &=\cdots \\ y_n &=y_0 + h\varphi_n \end{aligned} $$其中，$h$是离散时间点的时间步长，$t_n=t_0+nh$，$\varphi_i$是与$t_i$有关的递推公式。
根据龙格库塔方法的不同级别，$\varphi_i$也有不同的形式。
二、龙格库塔方法的应用由于龙格库塔方法的较高精度和鲁棒性，以及易于实现等特点，它在各个领域都有着广泛的应用。
1. 数学领域在数学领域，龙格库塔方法可以用于求解常微分方程、偏微分方程、常微分方程组等等，特别是对于复杂的高阶微分方程，龙格库塔方法更是可以发挥其优势。
2. 物理学领域在物理学领域，各种微分方程是研究物理过程的基础。
龙格库塔方法在求解各种物理问题时也得到了广泛的应用，如天体力学、流体力学、电磁场问题等等。
3. 经济学领域在经济学领域，许多经济问题可以通过微分方程的形式进行建模，并采用龙格库塔方法进行数值求解。
滤波算法 龙格库塔算法-概述说明以及解释
滤波算法龙格库塔算法-概述说明以及解释1.引言1.1 概述概述：滤波算法和龙格库塔算法是计算机科学领域中常用的算法之一，它们在数据处理和数值计算中有着重要的应用价值。
滤波算法被广泛应用于信号处理、图像处理、通信系统等领域，用于消除信号中的噪声和提高数据的质量。
而龙格库塔算法则是一种常用的数值求解微分方程的方法，能够有效地对复杂的数学模型进行数值求解，具有较高的准确性和稳定性。
本文将分别介绍滤波算法和龙格库塔算法的原理、优缺点以及应用领域，希望读者通过本文能够对这两种算法有更深入的了解，并在实际应用中能够灵活运用。
1.2 文章结构本文将分为四个部分来探讨滤波算法和龙格库塔算法。
首先在引言部分，对滤波算法和龙格库塔算法进行简要介绍，并说明本文的结构和目的。
接着在第二部分，详细介绍滤波算法的概念、常见算法和应用场景，以便读者对滤波算法有个全面的了解。
然后在第三部分，深入探讨龙格库塔算法的简介、原理和优缺点，帮助读者更好地理解这一种数值计算方法。
最后，在结论部分对两种算法进行总结，并展望未来可能的发展方向，以及得出结论。
通过以上四个部分的内容，读者能够全面了解和掌握滤波算法和龙格库塔算法的相关知识。
1.3 目的本文的主要目的是介绍和探讨滤波算法和龙格库塔算法这两种在计算机科学和工程领域中广泛应用的算法。
通过对这两种算法的概述、原理和应用进行详细分析，能够帮助读者全面了解它们的工作原理和特点。
同时，通过对这两种算法的比较和讨论，可以帮助读者更好地理解它们在不同应用场景下的适用性和优缺点。
此外，本文还旨在为读者提供一个深入学习和掌握这两种算法的基础知识和入门指南。
通过本文的学习，读者可以加深对滤波算法和龙格库塔算法的理解，为进一步的研究和实践打下坚实的基础。
同时，希望本文能够激发读者对算法领域的兴趣，促使他们深入研究和探索更多先进的算法及其应用。
2.滤波算法2.1 滤波算法概述滤波算法是一种用于处理信号或数据的技术，其主要目的是通过去除噪声或不需要的信息，从而提取出所需的信号或数据。
(完整版)第二节龙格-库塔方法
因为单步法是 p阶的：h0 , 0 h h0 满足| Tn1 | Chp1
| en1 || en | hL | en | Ch p1 | en |
其中 1 hL, Ch p1
en O(hp )
即取
K * 1K1 2 K2 yi1 yi h(1K1 2 K2 )
其中1 和 2 是待定常数。若取 K1 f ( xi , yi ) ，则
问题在于如何确定 xi p 处的斜率 K2 和常数 1 和 2 。
仿照改进的欧拉方法，用欧拉方法预测 y( xi p ) 的值，
yi p yi phK1
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x (0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
E(h) 1 h
绝对稳定域： 1 h 1
当 R 时， 1 h 1 2 h 0
绝对稳定区间：(2, 0)
❖经典的R-K公式：yn1
yn
h 6
(k1
2k2
2k3
k4 )
k1 f ( xn , yn ) yn
k2
f
(
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
龙格-库塔方法
h 按照数学学 习一般方法， yn 1 yn （k1 2k 2 2k3 +k 4） 6 继续上面的 做法，经过 k1 f ( xn , yn ) 较复杂的数 h 学推理和计 k （16） 2 f ( xn 1 , yn k1 ) 算，可以导 2 2 出四阶龙格 h －库塔格式， k f ( x 1 , yn k 2 ) 这里有一个 3 n 2 2 较有用的四 k f ( x , y hk ) 阶格式。 n 1 n 3 4
？
3、欧拉格式
yn 1 yn hf ( xn , yn ),
取点 xn 的斜率 k
*
(3)
n 0,1, 2,
f ( xn , yn )作为区间 [ xn , xn1 ] 上的平均斜率，精度低。
4、改进的欧拉格式
于是：
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn hf ( xn , yn ))] 2 h yn 1 yn [k1 k2 ] 2
龙格－库塔方法
考虑一维经典初值问题
dy f ( x , y ) , y( x0 ) y0 , x [a, b] dx
(1)
龙格－库塔方法
一、设计思想
加权平均斜率。
1、微分中值定理
如果 f (x) 在 [a, b]上连续，在 ( a, b) 内可导，则在( a, b) 内至 少存在一点 ，有 f (b) f (a) f ' ( ) ba
xnq xn qh, 0 p q 1
令
yn1 yn h[(1 )k1 k2 k3 ] (14)
k1 , k 2
第二节_龙格_库塔方法
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
如果 |1 hf y | 1，则误差是不增的，故可认为是稳定的
例如：对于初值问题
y y
y(
x0
)
a
精确解为
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x
(0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
h
y
ae x x0
y y
而实际求解的初值问题为
y(
x0
)
a
a
精确解为 y (a a)e x x0在 xn 处的误差为ae xn x0
可见误差随着 xn的增加呈指数函数增长
y y
如果初值问题为
y( x0 ) a
精确解为 y ae x0 x
y y
实际求解的初值问题为
y( x0 )
2C( h)5 2
y( xn1)
y(h/ 2) n1
y( xn1)
y(h) n1
1 16
记
|
y(h/ 2) n1
y(h) n1
|
16( y( xn1)
y(h/ 2) n1
)
y( xn1)
y(h) n1
y( xn1)
y(h/ 2) n1
四阶龙格库塔法原理
四阶龙格库塔法原理
四阶龙格库塔法是一种常见的数值计算方法，用于求解常微分方程的初值问题。
该方法通过逐步逼近准确解来得到数值解。
在四阶龙格库塔法中，我们将求解区间 [a, b] 平均分成 n 个子
区间，每个子区间的长度为 h = (b - a) / n，其中 a 是初始点，b 是终点。
首先，我们需要给出初始条件 y(a)，即方程在 a 点的值。
然后，我们利用以下公式依次计算 y(a + kh) 的近似值，其中 k 为当
前步数（从 0 开始）：
\[
\begin{align*}
k_1 &= hf(a, y(a)) \\
k_2 &= hf(a + \frac{h}{2}, y(a) + \frac{k_1}{2}) \\
k_3 &= hf(a + \frac{h}{2}, y(a) + \frac{k_2}{2}) \\
k_4 &= hf(a + h, y(a) + k_3) \\
y(a + kh) &= y(a) + \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}
\end{align*}
\]
其中，f(a, y(a)) 为方程的右端函数，描述了方程随时间的变化
规律。
在每一步中，我们利用当前步数的近似值 k1 到 k4 来计算下一步的近似值。
最后，通过不断迭代计算，直到达到指定的终点b，我们可以得到方程的数值解。
四阶龙格库塔法的主要原理是通过将步长 h 下的误差控制在O(h^5) 的级别，从而在保证计算效率的同时，提高数值解的精度。
该方法具有较好的稳定性和精确性，常被应用于各种科学计算和工程问题的数值求解中。
4.2 龙格库塔法
重复积分过程，直到相邻两次的积分值R2n 与Rn 满足 下列关系：
| R2n Rn | | R2n | 1
R2 n Rn | | R2 n
其中 为允许的误差限。
| R2n | 1
例1：利用龙贝格法计算单摆的振荡周期 T l /2 d T 4 1 g 0 2 0 [1 sin ( )sin 2 ]2 2 g 9.8m/s2 ， 其中 l 5m ， 0 20/ ， 105 program main external f call romb(0.,1.5708,1e-5,r2,f) write(*,"(3x,'B=',e12.6)") r2 End function f(x) f=4*sqrt(5/9.8)/(1-sin(10/3.14159)**2*sin(x)**2)**(1/2) end
subroutine romb(a,b,ep,r2,f) h=b-a t1=h/2.*(f(a)+f(b)) n=1 5 s=0.0 do k=0,n-1 s=s+f(a+(k+0.5)*h) end do t2=t1/2.+h/2.*s s2=t2+(t2-t1)/3. if(n/=1)goto 20 15 n=n+n h=h/2 t1=t2 s1=s2 goto 5 20 c2=s2+(s2-s1)/15. if(n/=2)goto 40
共有 n 1个分点 xk a kh, h b a , k 0,1,2,, n 用Tn 表示简单梯形法所求得的积分值。 1、考察一个子段 [ xk , xk 1 ], 其中点为 xk 1
2
n
1 ( xk xk 1 ), 2
变步长的龙格库塔法
每步须算Ki 的个数 2 3
4
5
6
可达到的最高精度 O (h2 ) O (h3 ) O(h4 ) O(h4 ) O(h5 )
7
O(h6 )
n8
O(hn2 )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
§2 Runge-Kutta Method
y y h[ K K ... K ]
i 1
i
11
22
mm
K f(x , y )
1
ii
K f ( x h, y hK )
2
i
2
i
21
1
K f ( x h, y hK hK )
3
i
3
i
31
1
32
2
... ...
K f ( x h, y hK hK ... hK )
从节点由于局部截断误差为故有然后将步长折半即以为步长从节点x出发跨两步到节点xi1再求得一个近似值截断误差是因此有这样16由此可得这表明以作为的近似值其误差可用步长折半前后两次计算结果的偏差当要求的数值精度为时
变步长的龙格—库塔方法
以 经 典 四 阶 龙 格 — 库 塔 公 式 为 例 。 从 节 点 x n出 发 ， 以 h为
115 (yi( h 21)
yi( h1))
h
这表明以
()
y2 i1
作为 y(xi1)的近似值，其误差可用步
长折半前后两次计算结果的偏差
来判断所选步长是否适当
y y ( h ) 2 i 1
龙格库塔积分算法
龙格库塔法龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。
龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。
其中4 阶龙格库塔法是最常用的一种方法。
因为它相当精确、稳定、容易编程。
在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。
如果需要较高的精度, 可采取减小步长的方法即可。
4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。
1、初值问题对于一阶常微分方程的初值问题根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。
2、离散化取步长h=(b-a)/n，将区间[a , b]分成n个子区间：a=<=b在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲线上至少有一点，它的斜率与整段曲线的平均斜率相同，得=y’() (0<<1)其中，=可以将上式改写成y()=y()+h*K (2.1)其中K为平均斜率，K=f()公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。
欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。
3、Euler法欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙格库塔法的基础。
首先，令、为y() 及y()的近似值，并且令平均斜率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1）4、改进的欧拉法此种方法是取、两点的斜率的平均值作为平均斜率K，即K= ,其中、均为y()以及y()的近似值，就得到改进后的欧拉公式（4.1）其中、分别为、两点的斜率值，即= ,=在上面的（4.1）式中，k2是未知的，采用一种叫预报法的方法来求解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
龙格库塔法
龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。
龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。
其中4 阶龙格库塔法是最常用的一种方法。
因为它相当精确、稳定、容易编程。
在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。
如果需要较高的精度, 可采取减小步长的方法即可。
4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。
1、初值问题
对于一阶常微分方程的初值问题
根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。
2、离散化
取步长h=(b-a)/n，将区间[a , b]分成n个子区间：
a=<=b
在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲
线上至少有一点，它的斜率与整段曲线的平均斜率相同，
得=y’() (0<<1)
其中，=
可以将上式改写成y()=y()+h*K (2.1)
其中K为平均斜率，K=f()
公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。
欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。
3、Euler法
欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙
格库塔法的基础。
首先，令、为y() 及y()的近似值，并且令平均斜
率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1）
4、改进的欧拉法
此种方法是取、两点的斜率的平均值作为平均斜率K，
即K= ,其中、均为y()以及y()的近似值，就得到
改进后的欧拉公式（4.1）
其中、分别为、两点的斜率值，即
= ,=
在上面的（4.1）式中，k2是未知的，采用一种叫预报法的方
法来求解。
即先用欧拉法求得一个初步的近似值，记作 ,称为预报值，预报值的精度不高，但是可以用它代替k2中的再直接计算，便可得到校正后的。
欧拉法求得预报值=，
校正=（）
将带入校正公式可得
改进后的欧拉公式
=（4.2）通常写作
（）
由此可见，为了提高精度，可以多取几点的斜率值作加权平均
来得到平均斜率K，这就是龙格库塔法的基本思想。
5、局部截断误差
初值问题
的单步法可用一般形式表示：
其中，多元函数与有关，如果中含有，那么方法是隐式的，若不含则为显式单步法，所以，显式单步法可以表示为
其中称为增量函数，对欧拉公式有
=f(x,y)
从开始计算，如果考虑每一步产生的误差，直到，则误差
，称该方法为在点的整体截断误差，分析求解整体截断误差是复杂的，为此，仅考虑从到的局部情况，并假设之前的计算没有误差，即，则一般显式单步法的局部截断误差可以定义为
一般，为了求解局部截断误差，我们会将对在
x=点做泰勒展开，进行化简。
6、龙格库塔法
龙格库塔法的基本思想就是用区间上的若干个点的导数
，将它们做加权平均得到平均斜率K，以提高方法的阶数，即提高方法的精度。
一般的显式龙格库塔法可以写成
其中，=
=f
（）
i=2,3......L
其中，，
它的局部截断误差是
其中，与的区别在于，用微分方程精确解代替中所
含的就得到，参数，，待定，确定这些常数的原则：使局
部截断误差误差关于h的阶数尽可能高。
二阶龙格库塔公式
令L=2，则
其局部截断误差为
（）（6.1）利用泰勒公式，将中各项作泰勒展开，特别是在
（，（））点做二元泰勒展开，并利用 ’，，=+y’()=y’’()
则有
+h*++O()
，
将上面各式带入（6.1）式中，整理得
选取、、使局部截断误差尽量小，就是使h和项的系
数为0，于是，得到方程组
这里是含有三个未知量的2个方程，它可以有无穷组解，都统
称为二阶龙格库塔公式。
取=0.5，则，=1，得到中点公式：
以上给出了构造2级龙格库塔公式的方法，它可以推广到L=3，4，5......的情况。
经典的龙格库塔法是一种4阶方法，计算公式是：
在具体计算中，利用已知的h，，，求得、、、
的值，代入公式中，就得到了的值。
7、变步长的龙格库塔法
单从每一步看，步长越小，截断误差就越小，但随着步长的减小，不但引起计算量的增大，而且可能导致舍入误差的严重积累，
因此，同积分的数值计算一样，微分方程的数值解法也要注意合理
的选择步长。
设从出发，以h为步长，经过一步计算，得到的近
似值，记作，若计算方法是p阶的，则有
（7.1）然后，将步长减半，即取为步长，从出发经两步计算，求得的近似值，记作
（7.2）当h很小时，以上两式中的c可看作同一常数，用乘式（7.2）的两端后与式（7.1）相减，得
因此，取=，显然要比
或的精度高。
当p=4时，可取
=
由此，可得下列近似式，
≈ =
上式表明，以作为的近似值，其误差可用先后两
次计算结果的差来表示。
可以用来判断所选取的步长是否合适，具
体可分为以下两种情况：
(1) 对于给定的精度，如果，就将反复将步长减半计算，直至为止，这时取步长折半后的“新值” 作为结果。
(2) 如果，就反复将步长加倍，直至为止，这时取上一次的步长的计算所得到的值作为。
这种通过步长减半或加倍来计算的方法就叫变步长法，表面上看，为了选择步长，每一步的计算量似乎增加了，但从总体上考虑往往是合理的。
8、龙格-库塔法在MATLAB中的用法
ode23:使用二阶龙格-库塔法求微分方程，调用格式为：
[t， ] od 3 ‘ n m ’ tsp n
ode45:使用四阶龙格-库塔法求微分方程，调用格式为：
[t， ] od 4 ‘ n m ’ tsp n
9、隐式龙格库塔公式
（
i=1,2,3......L
隐式龙格库塔法适用于解刚性问题，对于基于数值求积公式的隐式龙格库塔法，可以利用3种方法建立隐式龙格库塔方法，分别是基于高斯求积公式的隐式龙格库塔法，基于拉道求积公式的隐式龙格库塔法，基于洛巴托求积公式的隐式龙格库塔法。
而且一般的隐式龙格库塔法计算比较复杂，为了简化计算，可以考虑一种对角隐式龙格库塔法。
隐式龙格库塔法的稳定性函数
d t
d t ，这是z的有理函数。