四阶龙格库塔算法(ode45)中如何令自变量取特定值

MILAN3

示例程序：
[t,y]=ode45(f,[0 12],y0);
其中， f为函数句柄，自变量为范围为[0 12]，y0为y的初始值，计算时程序在自变量取值范围内自动取值，请教高手，能规
定自变量取一些特殊值么，如2、4、6、8、10、12，如果能，如何编程实现。

messenger

记得matlab faq里有这个问题，答案是不能。如果一定要这样，试试缩小时间区间

MILAN3

自问自答，将程序改为下式即可
[t,y]=ode45(f,[0:0.2:12],y0)
即设定自变量T以0.2的步长在[0 12]的范围内变化，这只是一个特例，因为特殊值刚好有规律，如果没有规律就不知如何解决了。这种设置的另一好处是确定了ODE45的迭代次数，这一设置中程序迭代了(12-0)/0.2=60 次。

inndoor

MILAN3 发表于 2013-8-7 14:25
自问自答，将程序改为下式即可
[t,y]=ode45(f,[0:0.2:12],y0)
即设定自变量T以0.2的步长在[0 12]的范围内变 ...

定步长积分：

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.';

复制代码

MILAN3

inndoor 发表于 2013-8-8 11:44
定步长积分：

非常感谢，计算发现，固定步长后，有时计算不收敛。

inndoor

MILAN3 发表于 2013-9-2 16:29
非常感谢，计算发现，固定步长后，有时计算不收敛。

龙格库塔法收敛与否和步长有点关系。