Download<< Back
function conv_analysis(method)
% COMV_ANALYSIS Performs convergence analysis of the specified Implicit method
%
% Argument: method - name of the Implicit method
% $$ 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)';
% Solve using the specified solver with $h=(\frac12)^p$, $p=1,\dots,10$,
% get the numerical solution at the time point T, and compute the error
% $\mathcal{E}_N = \Vert u(T)=u_N\Vert$:
error = zeros(no_approx,1);
for i=1:no_approx
un = solve(method, @logistic, 0, 2, 2, h(i), timepoint);
error(i) = norm(exact-un);
end
% We now plot the error vs. the time step $h$ on a log-log plot:
figure
loglog(h,error,'-ro');
title('Order')
xlabel('h');
ylabel('error');
% We use linear regression to estimate the convergence rate in the form
%
% $$ \log_{10}\mathcal{E}_N = q+p\log_{10} h $$,
%
% and, hence, estimate the order $p$ of the method.
[p,q]=regrese(h,error);
fprintf('log_{10}E_N = %f + %f x log10(error)\n', q, p);
end
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
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)];
r=A\data(:,2);
p=r(1);
q=r(2);
end
