Download
<< Back
OutputCode
%% Question 1 - Order of Implicit & Explicit methods
%% Exact solution
function q1_solution
%%
% We solve the logistic equation
% \[\begin{array}{rcll}
% x'(t) &=& (a-bx(t))x(t), &\qquad t\in[0,2], \\
% x(0) &=& x_0,
% \end{array}\]
% with $a=b=1$, $x_0=2$, which has the known exact solution
% \[
% x(t) = \frac{x_0 e^t}{1-x_0(1-e^t)}.
% \]
x0 = 2;
timepoint = 2;
exact = x0*exp(timepoint)/(1-x0*(1-exp(timepoint)));
no_approx = 10;
h = 0.5.^(1:no_approx)';
%% Euler
error = zeros(no_approx,1);
for ei=1:no_approx
x = solve(@eul, @logistic, 0, 2, timepoint, h(ei), timepoint);
error(ei) = norm(exact-x);
end
[p,q]=regrese(h,error);
figure;
plot(h,error,'-ro');
title('Euler');
xlabel('\tau');
ylabel('error');
figure;
loglog(h,error,'-ro');
title('Euler')
xlabel('\tau');
ylabel('error');
fprintf('Euler = %f + %f x log10(tau)\n', q, p);
%%
% We can see that for Euler $p\approx 1$.
%% Runge
error = zeros(no_approx,1);
for ei=1:no_approx
x = solve(@runge, @logistic, 0, 2, timepoint, h(ei), timepoint);
error(ei) = norm(exact-x);
end
[p,q]=regrese(h,error);
figure;
plot(h,error,'-ro');
title('Runge');
xlabel('\tau');
ylabel('error');
figure;
loglog(h,error,'-ro');
title('Runge')
xlabel('\tau');
ylabel('error');
fprintf('Runge = %f + %f x log10(tau)\n', q, p);
%%
% We can see that for Runge $p\approx 2$.
%% Runge-Kutta
error = zeros(no_approx,1);
for ei=1:no_approx
x = solve(@rk_classical, @logistic, 0, 2, timepoint, h(ei), timepoint);
error(ei) = norm(exact-x);
end
[p,q]=regrese(h,error);
figure;
plot(h,error,'-ro');
title('Runge-Kutta');
xlabel('\tau');
ylabel('error');
figure;
loglog(h,error,'-ro');
title('Runge-Kutta')
xlabel('\tau');
ylabel('error');
fprintf('Runge-Kutta = %f + %f x log10(tau)\n', q, p);
%%
% We can see that for Runge-Kutta $p\approx 4$.
%% Heun
error = zeros(no_approx,1);
for ei=1:no_approx
x = solve(@heun, @logistic, 0, 2, timepoint, h(ei), timepoint);
error(ei) = norm(exact-x);
end
[p,q]=regrese(h,error);
figure;
plot(h,error,'-ro');
title('Heun');
xlabel('\tau');
ylabel('error');
figure;
loglog(h,error,'-ro');
title('Heun')
xlabel('\tau');
ylabel('error');
fprintf('Heun = %f + %f x log10(tau)\n', q, p);
%%
% We can see that for Heun $p\approx 2$.
%% Implicit Euler
% Note, we start from h=1/4, as Implicit Euler fails to
% converge for h=1/2
error = zeros(no_approx,1);
h = 0.5.^(1+(1:no_approx))'; % start from 1/4 (.5^2)
for ei=1:no_approx
x = solve(@ieuler, @logistic, 0, 2, timepoint, h(ei), timepoint);
error(ei) = norm(exact-x);
end
[p,q]=regrese(h,error);
figure;
plot(h,error,'-ro');
title('Implicit Euler');
xlabel('\tau');
ylabel('error');
figure;
loglog(h,error,'-ro');
title('Implicit Euler')
xlabel('\tau');
ylabel('error');
fprintf('Implicit Euler = %f + %f x log10(tau)\n', q, p);
%%
% We can see that for Implicit Euler $p\approx 1$.
%% Crank-Nicholson
error = zeros(no_approx,1);
for ei=1:no_approx
x = solve(@crank_nicholson, @logistic, 0, 2, timepoint, h(ei), timepoint);
error(ei) = norm(exact-x);
end
[p,q]=regrese(h,error);
figure;
plot(h,error,'-ro');
title('Crank-Nicholson');
xlabel('\tau');
ylabel('error');
figure;
loglog(h,error,'-ro');
title('Crank-Nicholson')
xlabel('\tau');
ylabel('error');
fprintf('Crank-Nicholson = %f + %f x log10(tau)\n', q, p);
%%
% We can see that for Crank-Nicholson $p\approx 2$.
%% Solve function
% Function to run the specified solver and extract the value at the
% specified timepoint
function y = solve(solver, field, t0, T, x0, h, timepoint)
%SOLVE Use the specified ODE solver to solve the ODE and get value at
% specified time point
%
% Parameters:
% solver -- ODE solve to use
% 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)
% timepoint -- Timepoint to get numerical solution at
%
% Returns: Computed value at timepoint
[t,x] = feval(solver, field, t0, T, x0, h);
n1 = floor((timepoint-t0)/h);
n2 = ceil((timepoint-t0)/h);
if n1 == n2
y = x(n1+1,:);
else
y1 = x(n1+1,:);
y2 = x(n2+1,:);
t1 = t(n1+1);
t2 = t(n2+1);
y = y1 + (y2-y1)*(timepoint-t1)/(t2-t1);
end
end
%% Linear Regression
% Function to perform linear regression in log-log coordinates
function [p,q]=regrese(x,y)
% REGRESE Compute linear regression
%
% Fits the line log10(y) = p log10(x) + q in "loglog" coordinates
%
% Returns: p, q
% p = estimated order of the method
data=[log10(x),log10(y)];
m=length(x);
A=[data(:,1), ones(m,1)];
coeffs=A\data(:,2);
p=coeffs(1);
q=coeffs(2);
end
end
