• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Ode_Hamming

UtamaDonny Dec 9th, 2019 (edited) 85 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. function [t,y] = ode_Ham(f,tspan,y0,N,KC,varargin)
2. % Hamming method to solve vector d.e. y’(t) = f(t,y(t))
3. % for tspan = [t0,tf] and with the initial value y0 and N time steps
4. % using the modifier based on the error estimate depending on KC = 1/0
5.
6. if nargin &lt; 5, KC = 1; end %with modifier by default
7. if nargin &lt; 4 | N &lt;= 0, N = 100; end %default maximum number of iterations
8. if nargin &lt; 3, y0 = 0; end %default initial value
9. y0 = y0(:)'; %make it a row vector
10. h = (tspan(2)-tspan(1))/N; %step size
11. tspan0 = tspan(1)+[0 3]*h;
12. [t,y] = ode_RK4(f,tspan0,y0,3,varargin{:}); %Initialize by Runge-Kutta
13. t = [t(1:3)' t(4):h:tspan(2)]';
14. for k = 2:4, F(k - 1,:) = feval(f,t(k),y(k,:),varargin{:}); end
15. p = y(4,:); c = y(4,:); h34 = h/3*4; KC11 = KC*112/121; KC91 = KC*9/121;
16. h312 = 3*h*[-1 2 1];
17.
18. for k = 4:N
19. p1 = y(k - 3,:) + h34*(2*(F(1,:) + F(3,:)) - F(2,:)); %Eq.(6.4.9a)
20. m1 = p1 + KC11*(c - p); %Eq.(6.4.9b)
21. c1 = (-y(k - 2,:) + 9*y(k,:) +...
22. h312*[F(2:3,:); feval(f,t(k + 1),m1,varargin{:})])/8; %Eq.(6.4.9c)
23. y(k+1,:) = c1 - KC91*(c1 - p1); %Eq.(6.4.9d)
24. p = p1; c = c1; %update the predicted/corrected values
25. F = [F(2:3,:); feval(f,t(k + 1),y(k + 1,:),varargin{:})];
26. end
27. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top