Code snippet: Runge-Kutta-Gil (Matlab)
From CFD-Wiki
(Difference between revisions)
m (Runge-Kutta-Gil (Matlab) moved to Code snippet: Runge-Kutta-Gil (Matlab)) |
Latest revision as of 08:58, 19 May 2006
function [tout,xout] = myrkg(fname,tspan,x0) % % Based on code written by: % Marc Compere % CompereM@asme.org % created : 06 October 1999 % modified: 17 January 2001 % % Copyright (C) 2001, 2000 Marc Compere % % rkg.m is free software; you can redistribute it and/or modify it % under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2, or (at your option) % any later version. % % <krivilli at unberwoot.net> - Modified to use Runge-Kutta-Gil method (06/03/2006): % without adaptive step, possible stiffness related issues. % % Kunge-Kutta-Gil coefficients w(1)=1/6; w(2)=(2-sqrt(2))/6; w(3)=(2+sqrt(2))/6; w(4)=1/6; c(1)=0; c(2)=1/2; c(3)=1/2; c(4)=1; a(2,1)=1/2; a(3,1)=-(sqrt(2)-1)/2; a(3,2)=(2-sqrt(2))/2; a(4,1)=0; a(4,2)=-sqrt(2)/2; a(4,3)=(2+sqrt(2))/2; t0 = tspan(1); tfinal = tspan(2); t = t0; h = (tfinal - t)/500.0; % '500' is our step. Change it as needed. x = x0(:); % column vector tout = t; xout = x.'; while (t < tfinal) if t + h > tfinal, h = tfinal - t; end ydot=feval(fname,t,x); k(:,1)=h*ydot; ydot=feval(fname,t+c(2)*h,x+a(2,1)*k(:,1)); k(:,2)=h*ydot; ydot=feval(fname,t+c(3)*h,x+a(3,1)*k(:,1)+a(3,2)*k(:,2)); k(:,3)=h*ydot; ydot=feval(fname,t+c(4)*h,x+a(4,2)*k(:,2)+a(4,3)*k(:,3)); k(:,4)=h*ydot; x2=x+w(1)*k(:,1)+w(2)*k(:,2)+w(3)*k(:,3)+w(4)*k(:,4); t = t + h; x = x2; tout = [tout; t]; xout = [xout; x.']; end;