龙格库塔.doc

数值分析中，龙格库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔龙格和马丁威尔海姆库塔于1900年左右发明。经典四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。令初值问题表述如下。则，对于该问题的RK4由如下方程给出：其中这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均： k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值； k3也是中点的斜率，但是这次采用斜率k2决定y值； k4是时间段终点的斜率，其y值用k3决定。 当四个斜率取平均时，中点的斜率有更大的权值：RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。注意上述公式对于标量或者向量函数(y可以是向量)都适用。显式龙格库塔法显示龙格库塔法是上述RK4法的一个推广。它由下式给出其中(注意：上述方程在不同著述中由不同但却等价的定义)。要给定一个特定的方法，必须提供整数s (阶段数)，以及系数 aij (对于1 j i s), bi (对于i = 1, 2, ., s)和ci (对于i = 2, 3, ., s)。这些数据通常排列在一个助记工具中，称为龙格库塔表：0c2a21c3a31a32csas1as2as,s 1b1b2bs 1bs龙格库塔法是自洽的，如果如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(hp+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.) end tryf0 = feval(odefun,tspan(1),y0,varargin:);catchmsg = Unable to evaluate the ODEFUN at t0,y0. ,lasterr;error(msg); end y0 = y0(:); % Make a column vector.if isequal(size(y0),size(f0)error(Inconsistent sizes of Y0 and f(t0,y0).);end neq = 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.) end tryf0 = feval(odefun,tspan(1),y0,varargin:);catchmsg = Unable to evaluate the ODEFUN at t0,y0. ,lasterr;error(msg); end y0 = y0(:); % Make a column vector.if isequal(size(y0),size(f0)error(Inconsistent sizes of Y0 and f(t0,y0).);end neq = 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.) end tryf0 = feval(odefun,tspan(1),y0,varargin:);catchmsg = Unable to evaluate the ODEFUN at t0,y0. ,lasterr;error(msg); end y0 = y0(:); % Make a column vector.if isequal(size(y0),size(f0)error(Inconsistent sizes of Y0 and f(t0,y0).);end neq = length(y0);N = length(tspan);Y = zeros(neq,N);% Method coefficients - Butchers 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, 09017/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:);end Y(:,i) = yi + F*(hi*B);endY = Y.;-自适应步长RKF45(相当于ode45)ODE45 是 4阶方法提供候选解，5阶方法控制误差。自适应步长RKF45(From J.K.Mathews，我做了小幅度修改，去掉了我认为不必要的测试语句，加了些注释):CODE:function R=RungeKutta45(f,a,b,ya,m,tol)% Input - f function% - a left endpoint of a,b% -b right endpoint of a,b% - ya initial value% -m initial guess for number of steps% Output - T solution: vector of abscissas% - Y solution: vector of ordinates% NUMERICAL METHODS: Matlab Programs% (c) 2004 by John H. Mathews and Kurtis D. Fink% Complementary Software to accompany the textbook:% NUMERICAL METHODS: Using Matlab, Fourth Editiona2 = 1/4; b2 = 1/4; a3 = 3/8; b3 = 3/32; c3 = 9/32; a4 = 12/13;b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197; a5 = 1;b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104; a6 = 1/2;b6 = -8/27; c6 = 2; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40;r1 = 1/360; r3 = -128/4275; r4 = -2197/75240; r5 = 1/50;r6 = 2/55; n1 = 25/216; n3 = 1408/2565; n4 = 2197/4104; n5 = -1/5;big = 1e15;h = (b-a)/m;hmin = h/64;% 步长自适应范围下限hmax = 64*h;% 步长自适应范围上限max1 = 200;% 迭代次数上限Y(1) = ya;T(1) = a;j = 1;tj = T(1);br = b - 0.00001*abs(b);while (T(j)br), h = b - T(j); endtj = T(j);yj = Y(j);y1 = yj;k1 = h*f(tj,y1);y2 = yj+b2*k1; if bigabs(y2) return, endk2 = h*f(tj+a2*h,y2);y3 = yj+b3*k1+c3*k2; if bigabs(y3) return, endk3 = h*f(tj+a3*h,y3);y4 = yj+b4*k1+c4*k2+d4*k3; if bigabs(y4) return, endk4 = h*f(tj+a4*h,y4);y5 = yj+b5*k1+c5*k2+d5*k3+e5*k4; if bigabs(y5) return, endk5 = h*f(tj+a5*h,y5);y6 = yj+b6*k1+c6*k2+d6*k3+e6*k4+f6*k5; if bigabs(y6) return, endk6 = h*f(tj+a6*h,y6);err = abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);ynew = yj+n1*k1+n3*k3+n4*k4+n5*k5;if (errtol)|(hbr),T(j+1) = b;elseT(j+1) = tj + h;endj = j+1;tj = T(j); endif (err=0),s = 0;elses = 0.84*(tol*h/err)(0.25);% 最佳步长值endif (s2*hmin), h = h/2; endif (s1.50)&(2*hhmax), h = 2*h; endif (bigabs(Y(j)|(max1=j), return, endendR=T Y;经我测试，结果均准确。后一个程序的精度较高。同步长大致相差一个数量级。-龙格-库塔法求解微分方程function yy=weifen_lgkt(A,B,x0,y0,fx)% yy=weifen_lgkt(A,B,x0,y0,fx)% 系统状态方程% dy/dx-Ay=B*fx 微分方程,A,B为参数% fx x的函数,可选参数,如不输入该参数, 则要在该函数体内编辑该函数% x0,y0 状态变量初始值%注：x0,y0可以为列向量,求解微分方程组N=1000; %计算步数h=0.002; %计算步长x=x0;y=y0;%-for k=1:N%求解微分方程%fx(:,k)=x0+(k-1)*h; %例子1，在此选f(x)=x k0=A*y+B*fx(:,k); %四阶龙格-库塔法则求解微分方程k1=A*(y+h*k0/2)+B*fx(:,k);k2=A*(y+h*k1/2)+B*fx(:,k);k3=A*(y+h*k2)+B*fx(:,k);y=y+(k0+2*k1+2*k2+k3)*h/6;%-%yy(k)=y; %输出xx(k)=x0+k*h;%y_jiexi(k)=2*fx(:,k)-2+4*exp(-fx(:,k); %例子1的解析解%y_e(k)=y_jiexi(k)-yy(k); %例子1解析解与数值解的误差endplot(x0 xx,y0 yy)程序完毕。说明：如求解方程dy/dx+y=2x 初值x0=0 y0=2 则选A=-1 B=2 x0=0 y0=2 fx=x 见程序内例1实现方法 ，当然也可把fx当作参数输进去。调用函数即可（例1需把%fx(:,k)=x0+(k-1)*h;前的%去掉）A=-1;B=2;x0=0;y0=2;yy=weifen_lgkt(A,B,x0,y0)就可以求得该微分方程龙格库塔法RKF45的Matlab实现 4阶5级龙格库塔法用于解一阶微分方程（组），对于高阶微分方程，可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写（数值方法matlab版），但只能解微分方程，不能解微分方程组。由LiuLiuuestc修改，使之能够解微分方程组。该程序精度比matlab自带的ode45更高。 rkf45.m:function Rt Rx=rkf45(f,tspan,ya,m,tol)% Input:% - f function column vector% - tspana,b left & right point of a,b% - ya initial value column vector% -m initial guess for number of steps% -tol tolerance% Output:% - Rt solution: vector of abscissas% - Rx solution: vector of ordinates% Program by John.Mathews, improved by liuliuuestcif length(tspan)=2error(length of vector tspan must be 2.);endif isnumeric(tspan)error(TSPAN should be a vector of integration steps.);endif isnumeric(ya)error(Ya 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.) ;end a=tspan(1);b=tspan(2);ya=ya(:);a2 = 1/4; b2 = 1/4; a3 = 3/8; b3 = 3/32; c3 = 9/32; a4 = 12/13;b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197; a5 = 1;b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104; a6 = 1/2;b6 = -8/27; c6 = 2; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40;r1 = 1/360; r3 = -128/4275; r4 = -2197/75240; r5 = 1/50;r6 = 2/55; n1 = 25/216; n3 = 1408/2565; n4 = 2197/4104; n5 = -1/5;big = 1e15;h = (b-a)/m;hmin = h/64;% 步长自适应范围下限hmax = 64*h;% 步长自适应范围上限max1 = 200;% 迭代次数上限Y(1,:) = ya;T(1) = a;j = 1;% tj = T(1);br = b - 0.00001*abs(b);while (T(j)br), h = b - T(j); end%caculate values of k1.k6,y1.y6tj = T(j);yj = Y(j,:);y1 = yj;k1 = h*feval(f,tj,y1);y2 = yj+b2*k1; if bigabs(max(y2) return, endk2 = h*feval(f,tj+a2*h,y2);y3 = yj+b3*k1+c3*k2; if bigabs(max(y3) return, endk3 = h*feval(f,tj+a3*h,y3);y4 = yj+b4*k1+c4*k2+d4*k3; if bigabs(max(y4) return, endk4 = h*feval(f,tj+a4*h,y4);y5 = yj+b5*k1+c5.*k2+d5*k3+e5*k4; if bigabs(max(y5) return, endk5 = h*feval(f,tj+a5*h,y5);y6 = yj+b6*k1+c6.*k2+d6*k3+e6*k4+f6*k5; if bigabs(max(y6) return, endk6 = h*feval(f,tj+a6*h,y6);err = abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);ynew = yj+n1*k1+n3*k3+n4*k4+n5*k5;% error and step size controlif ( (errtol) | (hbr),T(j+1) = b;elseT(j+1) = tj + h;endj = j+1;tj = T(j); endif (max(err)=0),s = 0;elses1 = 0.84*(tol.*h./err).(0.25);% 最佳步长值s=min(s1);endif (s2*hmin), h = h/2; endif (s1.50)&(2*hhmax), h = 2*h; endif ( (bigabs(Y(j,:) | (max1=j) ), return, endend% Rt Rx=T Y;Rt=T;Rx=Y; 使用方法：首先编写方程（组）文件（注意与ode45不同，这儿方程组为1Xn数组：function dx= fun(t,x)dx=zeros(1,2);dx(1)=x(1)+x(2)*2+0*t;dx(2)=3*x(1)+x(2)*2+0*t;然后使用：Rt,Rx=rkf45(fun,0,0.2,6;4,100,1e-7)龙格库塔法的c+编程龙格库塔法的c+编程 CODE: #include #include /*n表示几等分，n+1表示他输出的个数*/ int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double,double) double h=(b-a)/n,k1,k2,k3,k4; int i; / x=(double*)malloc(n+1)*sizeof(double); / y=(double*)malloc(n+1)*sizeof(double); x0=a; y0=y0; switch(style) case 2: for(i=0;in;i+) xi+1=xi+h; k1=function(xi,yi); k2=function(xi+h/2,yi+h*k1/2); yi+1=yi+h*k2; break; case 3: for(i=0;in;i+) xi+1=xi+h; k1=function(xi,yi); k2=function(xi+h/2,yi+h*k1/2); k3=function(xi+h,yi-h*k1+2*h*k2); yi+1=yi+h*(k1+4*k2+k3)/6; break; case 4: for(i=0;in;i+) xi+1=xi+h; k1=function(xi,yi); k2=function(xi+h/2,yi+h*k1/2); k3=function(xi+h/2,yi+h*k2/2); k4=function(xi+h,yi+h*k3); yi+1=yi+h*(k1+2*k2+2*k3+k4)/6; break; default: return 0; return 1; double function(double x,double y) return y-2*x/y; /例子求y=y-2*x/y(0x1);y0=1; /* int main() double x6,y6; printf(用二阶龙格-库塔方法n); RungeKutta(1,0,1,5,x,y,2,function); for(int i=0;i6;i+) printf(x%d=%f,y%d=%fn,i,xi,i,yi); printf(用三阶龙格-库塔方法n); RungeKutta(1,0,1,5,x,y,3,function); for(i=0;i6;i+) printf(x%d=%f,y%d=%fn,i,xi,i,yi); printf(用四阶龙格-库塔方法n); RungeKutta(1,0,1,5,x,y,4,function); for(i=0;i6;i+) printf(x%d=%f,y%d=%fn,i,xi,i,yi); return 1; 龙格库塔求解微分方程数值解工程中很多的地方用到龙格库塔求解微分方程的数值解， 龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。 此代码只是演示求一个微分方程 /*y=y-2x/y,x0,0.6 /*y(0)=1 的解，要求解其它的微分方程，可以自己定义借口函数，退换程序里面的函数： float f(float , float)即可； CODE: /* /*编程者： /*完成日期：XXXX，XX，XX /*e.g:y=y-2x/y,x0,0.6 /*y(0)=1 /*使用经典四阶龙格-库塔算法进行高精度求解 /* y(i+1)=yi+( K1+ 2*K2 +2*K3+ K4)/6 /* K1=h*f(xi,yi) /* K2=h*f(xi+h/2,yi+K1/2) /* K3=h*f(xi+h/2,yi+K2/2) /* K4=h*f(xi+h,yi+K3) */ #include #include float f(float den,float p0) /要求解的微分方程的右部的函数 e.g: y-2x/y float rus; / den=w/(W0+sl); / rus=k*A*f/Wp*sqrt(RTp)-(k-1)*. rus=p0-2*den/p0; return(rus); /float fden() / / void main() float x0; /范围上限 int x1; /范围下限 float h; /步长 int n; /计算出的点的个数 float k1,k2,k3,k4; float y3; /用于存放计算出的常微分方程数值解 int i=0; int j; /以下为函数的接口 printf(intput x0:); scanf(%f,&x0); printf(input x1:); scanf(%f,&x1); printf(input y0:); scanf(%f,&y0); printf(input h:); scanf(%f,&h); / 以下为核心程序 n=(x1-x0)/h); n=3; for(j=0;jn;j+) k1=h*f(x0,yi); /求K1 k2=h*f(x0+h/2),(yi+k1/2); /求K2 k3=h*f(x0+h/2),(yi+k2/2); /求K3 k4=h*f(x0+h),(yi+k3); /求K4 yi+1=yi+(k1+2*k2+2*k3+k4)/6); /求yi+1 x0+=float(0.2); printf(y%f=%fn,x0,yi+1); +i; 定义 公