Download<< Back
function [t,x] = ieuler_newton(ode, t0, T, x0, h)
%IEULER_NEWTON Implements the Implicit Euler one-step ODE solver using Newton
%
% Parameters:
% ode -- Function which returns two function handles:
% - Right hand side function of ODE system: x'=f(t,x)
% - Derivative of right hand side w.r.t second (space) argument f_x(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
% Tolerance to use for solving the Newton iteration
tol = 0.01;
n = ceil((T-t0)/h);
t = t0+h*(0:n).';
x = ones(n+1,length(x0));
x(1,:) = x0;
I = eye(length(x0));
[f, df] = feval(ode);
for i=1:n
k1_old = feval(f, t(i), x(i,:).'); % k1_{(0)}
F = k1_old - feval(f, t(i)+h, x(i,:).'+h*k1_old);
Fx = I - h*feval(df, t(i)+h, x(i,:).'+h*k1_old);
k1_new = k1_old - Fx\F; % k1_{(1)}
while norm(k1_old-k1_new, 'inf') >= tol
k1_old = k1_new; % k1_{(n)}
F = k1_old - feval(f, t(i)+h, x(i,:).'+h*k1_old);
Fx = I - h*feval(df, t(i)+h, x(i,:).'+h*k1_old);
k1_new = k1_old - Fx\F; % k1_{(n+1)}
end
x(i+1,:) = x(i,:).' + h*k1_new;
end
