[分享]定步长龙格库塔法程序（三阶、四阶、五阶）

happy

经常看到很多朋友问定步长的龙格库塔法设置问题，下面吧定步长三阶、四阶、五阶龙格库塔程序贴出来，有需要的可以看看


ODE3 三阶龙格－库塔法

1. 
   
2. function Y = ode3(odefun,tspan,y0,varargin)
   
3. %ODE3 Solve differential equations with a non-adaptive method of order 3.
   
4. % Y = ODE3(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates
   
5. % the system of differential equations y' = f(t,y) by stepping from T0 to
   
6. % T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
   
7. % The vector Y0 is the initial conditions at T0. Each row in the solution
   
8. % array Y corresponds to a time specified in TSPAN.
   
9. %
   
10. % Y = ODE3(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters
    
11. % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).
    
12. %
    
13. % This is a non-adaptive solver. The step sequence is determined by TSPAN
    
14. % but the derivative function ODEFUN is evaluated multiple times per step.
    
15. % The solver implements the Bogacki-Shampine Runge-Kutta method of order 3.
    
16. %
    
17. % Example
    
18. % tspan = 0:0.1:20;
    
19. % y = ode3(@vdp1,tspan,[2 0]);
    
20. % plot(tspan,y(:,1));
    
21. % solves the system y' = vdp1(t,y) with a constant step size of 0.1,
    
22. % and plots the first component of the solution.
    
23. %
    
24. 
    
25. if ~isnumeric(tspan)
    
26. error('TSPAN should be a vector of integration steps.');
    
27. end
    
28. 
    
29. if ~isnumeric(y0)
    
30. error('Y0 should be a vector of initial conditions.');
    
31. end
    
32. 
    
33. h = diff(tspan);
    
34. if any(sign(h(1))*h <= 0)
    
35. error('Entries of TSPAN are not in order.')
    
36. end
    
37. 
    
38. try
    
39. f0 = feval(odefun,tspan(1),y0,varargin{:});
    
40. catch
    
41. msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
    
42. error(msg);
    
43. end
    
44. 
    
45. y0 = y0(:); % Make a column vector.
    
46. if ~isequal(size(y0),size(f0))
    
47. error('Inconsistent sizes of Y0 and f(t0,y0).');
    
48. end
    
49. 
    
50. neq = length(y0);
    
51. N = length(tspan);
    
52. Y = zeros(neq,N);
    
53. F = zeros(neq,3);
    
54. 
    
55. Y(:,1) = y0;
    
56. for i = 2:N
    
57. ti = tspan(i-1);
    
58. hi = h(i-1);
    
59. yi = Y(:,i-1);
    
60. F(:,1) = feval(odefun,ti,yi,varargin{:});
    
61. F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
    
62. F(:,3) = feval(odefun,ti+0.75*hi,yi+0.75*hi*F(:,2),varargin{:});
    
63. Y(:,i) = yi + (hi/9)*(2*F(:,1) + 3*F(:,2) + 4*F(:,3));
    
64. end
    
65. Y = Y.';
    

复制代码

[ 本帖最后由 suffer 于 2006-10-9 19:17 编辑 ]

happy

ODE4 四阶龙格－库塔法

1. function Y = ode4(odefun,tspan,y0,varargin)
   
2. %ODE4 Solve differential equations with a non-adaptive method of order 4.
   
3. % Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates
   
4. % the system of differential equations y' = f(t,y) by stepping from T0 to
   
5. % T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
   
6. % The vector Y0 is the initial conditions at T0. Each row in the solution
   
7. % array Y corresponds to a time specified in TSPAN.
   
8. %
   
9. % Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters
   
10. % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).
    
11. %
    
12. % This is a non-adaptive solver. The step sequence is determined by TSPAN
    
13. % but the derivative function ODEFUN is evaluated multiple times per step.
    
14. % The solver implements the classical Runge-Kutta method of order 4.
    
15. %
    
16. % Example
    
17. % tspan = 0:0.1:20;
    
18. % y = ode4(@vdp1,tspan,[2 0]);
    
19. % plot(tspan,y(:,1));
    
20. % solves the system y' = vdp1(t,y) with a constant step size of 0.1,
    
21. % and plots the first component of the solution.
    
22. %
    
23. 
    
24. if ~isnumeric(tspan)
    
25. error('TSPAN should be a vector of integration steps.');
    
26. end
    
27. 
    
28. if ~isnumeric(y0)
    
29. error('Y0 should be a vector of initial conditions.');
    
30. end
    
31. 
    
32. h = diff(tspan);
    
33. if any(sign(h(1))*h <= 0)
    
34. error('Entries of TSPAN are not in order.')
    
35. end
    
36. 
    
37. try
    
38. f0 = feval(odefun,tspan(1),y0,varargin{:});
    
39. catch
    
40. msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
    
41. error(msg);
    
42. end
    
43. 
    
44. y0 = y0(:); % Make a column vector.
    
45. if ~isequal(size(y0),size(f0))
    
46. error('Inconsistent sizes of Y0 and f(t0,y0).');
    
47. end
    
48. 
    
49. neq = length(y0);
    
50. N = length(tspan);
    
51. Y = zeros(neq,N);
    
52. F = zeros(neq,4);
    
53. 
    
54. Y(:,1) = y0;
    
55. for i = 2:N
    
56. ti = tspan(i-1);
    
57. hi = h(i-1);
    
58. yi = Y(:,i-1);
    
59. F(:,1) = feval(odefun,ti,yi,varargin{:});
    
60. F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
    
61. F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:});
    
62. F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});
    
63. Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));
    
64. end
    
65. Y = Y.';

复制代码

[ 本帖最后由 suffer 于 2006-10-9 19:18 编辑 ]

happy

ODE5 五阶龙格－库塔法

1. function Y = ode5(odefun,tspan,y0,varargin)
   
2. %ODE5 Solve differential equations with a non-adaptive method of order 5.
   
3. % Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates
   
4. % the system of differential equations y' = f(t,y) by stepping from T0 to
   
5. % T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
   
6. % The vector Y0 is the initial conditions at T0. Each row in the solution
   
7. % array Y corresponds to a time specified in TSPAN.
   
8. %
   
9. % Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters
   
10. % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).
    
11. %
    
12. % This is a non-adaptive solver. The step sequence is determined by TSPAN
    
13. % but the derivative function ODEFUN is evaluated multiple times per step.
    
14. % The solver implements the Dormand-Prince method of order 5 in a general
    
15. % framework of explicit Runge-Kutta methods.
    
16. %
    
17. % Example
    
18. % tspan = 0:0.1:20;
    
19. % y = ode5(@vdp1,tspan,[2 0]);
    
20. % plot(tspan,y(:,1));
    
21. % solves the system y' = vdp1(t,y) with a constant step size of 0.1,
    
22. % and plots the first component of the solution.
    
23. 
    
24. if ~isnumeric(tspan)
    
25. error('TSPAN should be a vector of integration steps.');
    
26. end
    
27. 
    
28. if ~isnumeric(y0)
    
29. error('Y0 should be a vector of initial conditions.');
    
30. end
    
31. 
    
32. h = diff(tspan);
    
33. if any(sign(h(1))*h <= 0)
    
34. error('Entries of TSPAN are not in order.')
    
35. end
    
36. 
    
37. try
    
38. f0 = feval(odefun,tspan(1),y0,varargin{:});
    
39. catch
    
40. msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
    
41. error(msg);
    
42. end
    
43. 
    
44. y0 = y0(:); % Make a column vector.
    
45. if ~isequal(size(y0),size(f0))
    
46. error('Inconsistent sizes of Y0 and f(t0,y0).');
    
47. end
    
48. 
    
49. neq = length(y0);
    
50. N = length(tspan);
    
51. Y = zeros(neq,N);
    
52. 
    
53. % Method coefficients -- Butcher's tableau
    
54. %
    
55. % C | A
    
56. % --+---
    
57. % | B
    
58. 
    
59. C = [1/5; 3/10; 4/5; 8/9; 1];
    
60. 
    
61. A = [ 1/5, 0, 0, 0, 0
    
62. 3/40, 9/40, 0, 0, 0
    
63. 44/45 -56/15, 32/9, 0, 0
    
64. 19372/6561, -25360/2187, 64448/6561, -212/729, 0
    
65. 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];
    
66. 
    
67. B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
    
68. 
    
69. % More convenient storage
    
70. A = A.';
    
71. B = B(:);
    
72. 
    
73. nstages = length(B);
    
74. F = zeros(neq,nstages);
    
75. 
    
76. Y(:,1) = y0;
    
77. for i = 2:N
    
78. ti = tspan(i-1);
    
79. hi = h(i-1);
    
80. yi = Y(:,i-1);
    
81. 
    
82. % General explicit Runge-Kutta framework
    
83. F(:,1) = feval(odefun,ti,yi,varargin{:});
    
84. for stage = 2:nstages
    
85. tstage = ti + C(stage-1)*hi;
    
86. ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));
    
87. F(:,stage) = feval(odefun,tstage,ystage,varargin{:});
    
88. end
    
89. Y(:,i) = yi + F*(hi*B);
    
90. 
    
91. end
    
92. Y = Y.';

复制代码

[ 本帖最后由 suffer 于 2006-10-9 19:18 编辑 ]

lc622503

这程序是不是可以直接调用
来计算微分方程组阿
function [x,y]=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
[x,y]=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
不知道什么意思
程序是想用欧拉法解微分方程
谢谢大虾指点阿

ziding1763

xf是没有赋值啊。
那怎么计算？
这个函数可以当成子函数，进行调用使用。
当然要在主程序中调用了。

suffer

原帖由 ziding1763 于 2006-10-23 15:28 发表
xf是没有赋值啊。
那怎么计算？
这个函数可以当成子函数，进行调用使用。
当然要在主程序中调用了。

同意，实际上这几个程序函数的调用是和ode45一样的

faith

这个程序我在M文件中运行总是出现错误，请问斑竹怎么回事啊？

suffer

原帖由 faith 于 2006-12-3 20:41 发表
这个程序我在M文件中运行总是出现错误，请问斑竹怎么回事啊？

是你调用的问题，最好把你的具体情况说一下，当然最理想的就是给代码了

shenyongjun

楼主的这几个程序是那里搞到的？我的Matlab7中没有这几个函数啊？
另问：Matlab7中的定步长数值积分函数？

stephenhope

楼上你 help ode23
你就知道了

faithdy

楼主，我有个难题想请教；和这个ode45类似，请赐教！万分感谢！
结合自己的课题或学过的课程中的微分方程动态模型，编写四阶龙格库塔数值积分程序求解

insects

　　有没有谁会编个二阶用龙格库塔法解非线性微分方程的程序给我看看，用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

insects

我原来就是用ode45编的程序，但是老师要我用这个龙格库塔法再编一个，这个方法还不会用，请指教。

insects

运行了楼主的ode4的程序就出现下面这个了：
??? Input argument "tspan" is undefined.

Error in ==> ode4 at 24
if ~isnumeric(tspan)
怎么改啊！

eight

建议先好好看看基础书，欲速则不达

[ 本帖最后由 ChaChing 于 2010-5-14 21:18 编辑 ]