[分享]定步长龙格库塔法程序（三阶、四阶、五阶）
经常看到很多朋友问定步长的龙格库塔法设置问题，下面吧定步长三阶、四阶、五阶龙格库塔程序贴出来，有需要的可以看看ODE3 三阶龙格－库塔法
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 = 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,);
% 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.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['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:N
ti = 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));
end
Y = Y.';
[ 本帖最后由 suffer 于 2006-10-9 19:17 编辑 ]
回复：（happy）[分享]定步长龙格库塔法程序（三阶、四...
ODE4 四阶龙格－库塔法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 = 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,);
% 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.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['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:N
ti = 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));
end
Y = Y.';
[ 本帖最后由 suffer 于 2006-10-9 19:18 编辑 ]
回复：（happy）[分享]定步长龙格库塔法程序（三阶、四...
ODE5 五阶龙格－库塔法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 = 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,);
% 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.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['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 -- Butcher's tableau
%
% C | A
% --+---
% | B
C = ;
A = [ 1/5, 0, 0, 0, 0
3/40, 9/40, 0, 0, 0
44/45 -56/15, 32/9, 0, 0
19372/6561, -25360/2187, 64448/6561, -212/729, 0
9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];
B = ;
% More convenient storage
A = A.';
B = B(:);
nstages = length(B);
F = zeros(neq,nstages);
Y(:,1) = y0;
for i = 2:N
ti = tspan(i-1);
hi = h(i-1);
yi = Y(:,i-1);
% General explicit Runge-Kutta framework
F(:,1) = feval(odefun,ti,yi,varargin{:});
for stage = 2:nstages
tstage = 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);
end
Y = Y.';
[ 本帖最后由 suffer 于 2006-10-9 19:18 编辑 ] 这程序是不是可以直接调用
来计算微分方程组阿
function =eulerpro(fun,x0,xf,y0,h)
n=fix((xf-x0)/h);x(1)=x0;y(1)=y0
for i=1:(n-1)
x(i+1)=x0+i*h; y1=y(i)+h*feval(fun,x(i),y(i))
y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2
end
function xEulerpro
clear all;clc
=eulerpro(@fun,0,1,1,0.1)
plot(x,y)
function f=fun(x,y)
f=-y+x+1
我编了这样一个程序但是现实这样的错误提示
??? Input argument "xf" is undefined.
Error in ==> haowan2 at 2
n=fix((xf-x0)/h);x(1)=x0;y(1)=y0
不知道什么意思
程序是想用欧拉法解微分方程
谢谢大虾指点阿 xf是没有赋值啊。
那怎么计算？
这个函数可以当成子函数，进行调用使用。
当然要在主程序中调用了。 原帖由 ziding1763 于 2006-10-23 15:28 发表
xf是没有赋值啊。
那怎么计算？
这个函数可以当成子函数，进行调用使用。
当然要在主程序中调用了。
同意，实际上这几个程序函数的调用是和ode45一样的
求教
这个程序我在M文件中运行总是出现错误，请问斑竹怎么回事啊？ 原帖由 faith 于 2006-12-3 20:41 发表这个程序我在M文件中运行总是出现错误，请问斑竹怎么回事啊？
是你调用的问题，最好把你的具体情况说一下，当然最理想的就是给代码了 楼主的这几个程序是那里搞到的？我的Matlab7中没有这几个函数啊？
另问：Matlab7中的定步长数值积分函数？ 楼上你 help ode23
你就知道了
编写四阶龙格库塔数值积分程序
楼主，我有个难题想请教；和这个ode45类似，请赐教！万分感谢！结合自己的课题或学过的课程中的微分方程动态模型，编写四阶龙格库塔数值积分程序求解
救命啊！
本帖最后由 wdhd 于 2016-5-31 09:27 编辑有没有谁会编个二阶用龙格库塔法解非线性微分方程的程序给我看看，用matlab语言的，方程是
d2θe(t)/dt2+(1/τ1+K*(τ2/τ1)*cosθe(t))* dθe(t)/dt+k/τ1*sinθe(t)-△ω0/τ1=0
二阶微分方程
我原来就是用ode45编的程序，但是老师要我用这个龙格库塔法再编一个，这个方法还不会用，请指教。有错误怎么办？
运行了楼主的ode4的程序就出现下面这个了：??? Input argument "tspan" is undefined.
Error in ==> ode4 at 24
if ~isnumeric(tspan)
怎么改啊！
回复 14楼 insects 的帖子
建议先好好看看基础书，欲速则不达[ 本帖最后由 ChaChing 于 2010-5-14 21:18 编辑 ]
页:
[1]
2