Matlab Program to Implement Mathematical Models Assignment Solution.

EULER SOLVER

function [tout, yout] = EulerSolver(f, t, y0)

% get number of points

N = length(t);

% Get step

h = t(2)-t(1);

m = length(y0);

% Define the matrix of size mxN where m is the number of variables

yout = zeros(m,N);

tout = t;

% Set the initial condition

yout(:,1) = y0;

% Now, calculate the rest of points

for i = 2:N

Ynew = yout(:,i-1) + h*f(t(i), yout(:,i-1));

yout(:,i) = Ynew;

end

end

IRK4 SOLVER

function [tout, yout] = IRK4Solver(f, t, y0)

% get number of points

N = length(t);

% Get step

h = t(2)-t(1);

m = length(y0);

% Define the matrix of size mxN where m is the number of variables

yout = zeros(m,N);

tout = t;

% Set the initial condition

yout(:,1) = y0;

% Now, calculate the rest of points

options = optimset('Display', 'Off');

for i = 2:N

[x, fval, ex] = fsolve(@(vars)solver(vars, t(i), f, yout(:,i-1), h), yout(:,i-1), options);

% Calculate terms

% k1 = f(t(i), yout(:,i-1));

% k2 = f(t(i) + 0.5*h, yout(:,i-1)+0.5*k1*h);

% k3 = f(t(i) + 0.5*h, yout(:,i-1) + 0.5*k2*h);

% k4 = f(t(i) + h, yout(:,i-1) + k3*h);

% yout(:,i) = yout(:,i-1) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);

yout(:,i) = x;

end

function ff = solver(vars, t, f, varsold, h)

k1 = f(t, vars);

k2 = f(t + 0.5*h, vars + 0.5*k1*h);

k3 = f(t + 0.5*h, vars + 0.5*k2*h);

k4 = f(t + h, vars + k3*h);

ff = vars - varsold - (h/6)*(k1 + 2*k2 + 2*k3 + k4);

end

end

LORENZ

function y = lorenz(t, x, sigma, rho, beta)

dxdt = sigma*(x(2)-x(1));

dydt = x(1)*(rho - x(3)) - x(2);

dzdt = x(1)*x(2) - beta*x(3);

y = [dxdt;dydt;dzdt];

end

RK4 SOLVER

function [tout, yout] = RK4Solver(f, t, y0)

% get number of points

N = length(t);

% Get step

h = t(2)-t(1);

m = length(y0);

% Define the matrix of size mxN where m is the number of variables

yout = zeros(m,N);

tout = t;

% Set the initial condition

yout(:,1) = y0;

% Now, calculate the rest of points

for i = 2:N

% Calculate terms

k1 = f(t(i), yout(:,i-1));

k2 = f(t(i) + 0.5*h, yout(:,i-1)+0.5*k1*h);

k3 = f(t(i) + 0.5*h, yout(:,i-1) + 0.5*k2*h);

k4 = f(t(i) + h, yout(:,i-1) + k3*h);

yout(:,i) = yout(:,i-1) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);

end

end