Personal Website:
Scott Congreve

Download<< Back

function one_step_order
%ONE_STEP_ORDER Experimental calculation of order of ODE methods

% problem: x'=(a-b*x)*x ,  x(0)=2;   a=1, b=1
% exact solution
x0 = 2;
timepoint = 2;
exact = x0*exp(timepoint)/(1-x0*(1-exp(timepoint)));

%--------------------------------------------------------------------------
% the Euler method: p ?
%--------------------------------------------------------------------------
% a sequence of approximations
no_approx = 10;
error = zeros(no_approx,1);
h = 0.5.^(1:no_approx)';
for i=1:no_approx
  x = solve(@eul, @logistic, 0, 2, timepoint, h(i), timepoint);
  error(i) = norm(exact-x);
end
figure
plot(h,error,'-ro');
title('order');
xlabel('\tau');
ylabel('error');

figure
loglog(h,error,'-ro');
title('order')
xlabel('\tau');
ylabel('error');
[p,q]=regrese(h,error);
fprintf('Euler = %f + %f x log10(tau)\n', q, p);

%--------------------------------------------------------------------------
% the Runge method: p ?
%--------------------------------------------------------------------------
% TODO


%--------------------------------------------------------------------------
% the Runge-Kutta method: p ?
%--------------------------------------------------------------------------
% TODO


%--------------------------------------------------------------------------
% the Heun method: p ?
%--------------------------------------------------------------------------
% TODO


%--------------------------------------------------------------------------
% the Implicit-Euler method: p ?
%--------------------------------------------------------------------------
% TODO
no_approx = 10;
error = zeros(no_approx,1);
h = 0.5.^(1+(1:no_approx))'; % start from 1/4 (.5^2)


%--------------------------------------------------------------------------
% the Crank-Nicholson method: p ?
%--------------------------------------------------------------------------
% TODO


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)];
coeffs=A\data(:,2);
p=coeffs(1);
q=coeffs(2);
end