The Three eighths rule in Matlab
The Three Eighths rule is a Runge-Kutta method with order 4. This is the code for a program written in Matlab for the initial value problem
y'=y
y(0)=1
We want to know the y value at t = 1. The obvious solution is
y(t)=et, and therefore y(1) = e
function [y,f_calls]=rk3_8(fun,t0,t1,y0,h,rpar)
t=t0;
y=y0;
fc=0;
while t < t1
if t+h>t1;h=t1-t;end
k1=feval(fun,t,y,rpar);
k2=feval(fun,t+h/3,y+h*k1/3,rpar);
k3=feval(fun,t+2*h/3,y+h*(k2-k1/3),rpar);
k4=feval(fun,t+h,y+h*(k3-k2+k1),rpar);
y=y+h*(k1+3*(k2+k3)+k4)/8;
fc=fc+4;
t=t+h;
end
if nargout > 1 fcalls=fc; end
function[y] = fun(x, y, rpar)
u=x;
Save both files at same directory, and form matlab console use the cd command to go this directory
Run the following command from console
>>y = rk3_8('fun', 0, 1, y0, 0.0001,0);
Now we have the result stored in the variable y, run this to see it
>> y
and that gives us the result
y = 2.718
Execute here the Three eighths rule method in our Runge kutta calculator
y'=y
y(0)=1
We want to know the y value at t = 1. The obvious solution is
y(t)=et, and therefore y(1) = e
Matlab source code for TheThree Eighths rule: file: rk3_8.m
function [y,f_calls]=rk3_8(fun,t0,t1,y0,h,rpar)
t=t0;
y=y0;
fc=0;
while t < t1
if t+h>t1;h=t1-t;end
k1=feval(fun,t,y,rpar);
k2=feval(fun,t+h/3,y+h*k1/3,rpar);
k3=feval(fun,t+2*h/3,y+h*(k2-k1/3),rpar);
k4=feval(fun,t+h,y+h*(k3-k2+k1),rpar);
y=y+h*(k1+3*(k2+k3)+k4)/8;
fc=fc+4;
t=t+h;
end
if nargout > 1 fcalls=fc; end
File fun.m
This is a file containing the function to evaluatefunction[y] = fun(x, y, rpar)
u=x;
Running
Save both files at same directory, and form matlab console use the cd command to go this directory
Run the following command from console
>>y = rk3_8('fun', 0, 1, y0, 0.0001,0);
Now we have the result stored in the variable y, run this to see it
>> y
and that gives us the result
y = 2.718
Execute here the Three eighths rule method in our Runge kutta calculator
