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 < 5, KC = 1; end %with modifier by default
  7. if nargin < 4 | N <= 0, N = 100; end %default maximum number of iterations
  8. if nargin < 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. OK, I Understand
 
Top