MATLAB EC

profilewaleed39
mat_lab_ec_codes.docx

function [t,y] = euler(f,tspan,y0,N)

% Solves the IVP y' = f(t,y), y(t0) = y0 in the time interval tspan = [t0,tf]

% using Euler's method with N time steps.

% Input:

% f = name of inline function or function M-file that evaluates the ODE

% (if not an inline function, use: euler(@f,tspan,y0,N))

% For a system, the f must be given as column vector.

% tspan = [t0, tf] where t0 = initial time value and tf = final time value

% y0 = initial value of the dependent variable. If solving a system,

% initial conditions must be given as a vector.

% N = number of steps used.

% Output:

% t = vector of time values where the solution was computed

% y = vector of computed solution values.

m = length(y0);

t0 = tspan(1);

tf = tspan(2);

h = (tf-t0)/N; % evaluate the time step size

t = linspace(t0,tf,N+1); % create the vector of t values

y = zeros(m,N+1); % allocate memory for the output y

y(:,1) = y0'; % set initial condition

for n=1:N

y(:,n+1) = y(:,n) + h*f(t(n),y(:,n)); % implement Euler's method

end

t = t'; y = y'; % change t and y from row to column vectors

end

function [t,y] = impeuler(f,tspan,y0,N)

% Solves the IVP y? = f(t,y), y(t0) = y0 in the time interval tspan=[t0,tf]

% using improved Euler's method with N time steps.

% Input:

% f = name of inline function or function M-file that evaluates the ODE

% (if not an inline function, use: euler(@f,tspan,y0,N))

% For a system, the f must be given as column vector.

% tspan = [t0, tf] where t0 = initial time value and tf = final time value

% y0 = initial value of the dependent variable. If solving a system,

% initial conditions must be given as a vector.

% N = number of steps used. % Output:

% t = vector of time values where the solution was computed

% y = vector of computed solution values.

m = length(y0);

t0 = tspan(1);

tf = tspan(2);

h = (tf-t0)/N; % evaluate the time step size

t = linspace(t0,tf,N+1); % create the vector of t values

y = zeros(m,N+1); % allocate memory for the output y

y(:,1) = y0'; % set initial condition

for n=1:N

f1=f(t(n),y(:,n));

f2=f(t(n+1),y(:,n)+h*f1);

y(:,n+1) = y(:,n) + h/2*(f1+f2); % implement improved Euler's method

end

t = t'; y = y'; % change t and y from row to column vectors

end