Download
<< Back
OutputCode
%% Question 2 - Heun
%% Test Heun
function q2_solution
a=1;
b=1;
logistic = @(t,x)(a-b*x)*x;
x0 = 2;
t0 = 0;
T = 3;
[t,x]=ode23(logistic, [t0, T], x0);
plot(t,x,'-k','DisplayName','ode23');
hold on;
[t,x]=heun(logistic, t0, T, x0, 1/2); plot(t,x,'-ro','DisplayName','Heun (\tau=1/2)');
[t,x]=heun(logistic, t0, T, x0, 1/4); plot(t,x,'-bx','DisplayName','Heun (\tau=1/4)');
[t,x]=heun(logistic, t0, T, x0, 1/8); plot(t,x,'-m+','DisplayName','Heun (\tau=1/8)');
xlabel('time');
ylabel('x');
legend('Location','NorthEast');
%% Heun function
function [t,x] = heun(field, t0, T, x0, h)
%HEUN Implements the Heun one-step ODE solver
%
% Parameters:
% field -- Right hand side function of ODE system: x'=f(t,x)
% t0 -- Initial time
% T -- End time (T > t0)
% x0 -- Initial value
% h -- Size of time step (h <= T-t0)
%
% Outputs:
% t -- [t0; t-0+h, t0+2*h; ...; t0+i*h; ...]
% x -- Matrix containing numerical solution, with each row the value of x
% at each time step
n = ceil((T-t0)/h);
t = t0+h*(0:n).';
x = ones(n+1,length(x0));
x(1,:) = x0;
for i=1:n
k1 = feval(field, t(i), x(i,:).');
k2 = feval(field, t(i)+h, x(i,:).'+h*k1);
x(i+1,:) = x(i,:).' + h*(k1+k2)/2;
end
end
end
