UtamaDonny

Ode_Hamming

Dec 9th, 2019
133
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