impeuler.m

ode routines in matlab

function [x, y] = ImpEuler( x0, x_end, y0, h, rhs_fctn )
%
%  Purpose:
%  Solve the initial value problem
%               
%                      y'(x) = f(x,y(x))
%                      y(x0) = y0
%
%  on the interval [ x0, x_end ] using the implicit Euler method
%  The implicit Euler method is given by
%
%    y(:,i+1) - h * f(x(i+1),y(:,i+1)) - y(:,i) = 0
%
%
%
%  Usage:
%
%       [x, y] = ImpEuler( x0, x_end, y0, h, rhs_fct, rhs_fct_y )
%
%  Parameters:
%
%  Input:
%
%  x0          initial x
%
%  x_end       final x
%
%  y0          initial values
%
%  h           step size
%
%  rhs_fctn    function evaluating f(x,y) for given x and y
%              rhs_fctn(x, y, '') returns  f(x,y)
%              rhs_fctn(x, y, 'jacobian') returns Jacobian df(x,y)/dy
%
%
%  Output:
%  
%  x        vector of x values x(i) = x0 + (i-1)*h
%
%  y        approximation of solution at x
%


% find number of time steps needed to reach final time
mx = ceil((x_end-x0)/ h);


% allocate workspace
n = size(y0(:),1);
x = (x0:h:x0+mx*h);
y = zeros(n,mx+1);

% tolerance for Newton's method
Newt_tol = 0.01*h^2;


% Set initial values
y(:,1) = y0(:);

for i = 1:mx

    % solve the nonlinear system for y(:,i+1) using Newton's 
    % method with starting value  y(:,i)

    y(:,i+1) = y(:,i);  % starting value for Newton's method
    Newt_it  = 0;       % counter for Newton iterations
    Newt_fct = y(:,i+1) - h * feval(rhs_fctn,x(i+1),y(:,i+1),'') - y(:,i);
    while( norm(Newt_fct) > Newt_tol*norm(y(:,i+1)) )
        % Compute Jacobian 
        Newt_Jac  = eye(n) - h * feval(rhs_fctn,x(i+1),y(:,i+1),'jacobian');
        Newt_step = - (Newt_Jac \ Newt_fct);   % Newton step
        y(:,i+1)  = y(:,i+1) + Newt_step;      % new Newton iterate
        Newt_it   = Newt_it + 1;      
        Newt_fct  = y(:,i+1) - h * feval(rhs_fctn,x(i+1),y(:,i+1),'') - y(:,i);
    end
    %disp([ x(i+1), Newt_it ])
end