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

引用 happy 2006-6-20 11:06

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 2006-6-20 11:07

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 2006-8-8 21:57

这程序是不是可以直接调用
来计算微分方程组阿
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 2006-10-23 15:28

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

引用 suffer 2006-11-29 09:24

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

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

引用 faith 2006-12-3 20:41

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

引用 suffer 2006-12-9 09:47

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

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

引用 shenyongjun 2006-12-12 14:32

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

引用 stephenhope 2006-12-15 14:48

楼上你 help ode23
你就知道了

引用 faithdy 2006-12-15 17:37

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

引用 insects 2007-5-16 10:13

本帖最后由 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

引用 insects 2007-5-16 10:16

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

引用 insects 2007-5-16 13:42

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

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

引用 eight 2007-5-16 15:25

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

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

引用 chenvy 2009-2-26 23:39

??? Error using ==> ode5
Too many output arguments.

Error in ==> Lorenzresponse at 19
    [T,x]=ode5(@Lorenz,tspan,yinit,[],yo1,yo2,yo3,yo4);
麻烦happy教授帮忙解释一下

引用 ansys10 2009-2-27 09:39

谢谢分享；
前段时间自己编的不行；解算例可以，实际一用就发散了，有空试试这个

引用 天涯万礼 2009-5-10 20:19

请问一下该程序可以用到电力系统吗？谢谢
3个发电机9个节点的  龙格库塔法:lol
连接  http://forum.vibunion.com/forum/viewthread.php?tid=81369
谢谢大哥帮帮忙 小弟非常感谢

[ 本帖最后由 ChaChing 于 2009-5-10 23:47 编辑 ]

引用 ldlnihao 2009-10-17 21:22

楼主，可以在定步长里加点东西让他变步长吗？期待中

引用 happy 2010-10-28 21:19

ldlnihao 发表于 2009-10-17 21:22

                               
登录/注册后可看大图

楼主，可以在定步长里加点东西让他变步长吗？期待中

现在的人的需求怎么都这么怪异
matlab提供变步长的，有人就非要定步长的

现在给出定步长的又想要改回变步长的

直接用matlab自带的ode函数不就完了？

查看全部评论(23)