MATLAB lessons N.5

A simple implementation of Euler scheme (vector form): File BVPsolver.m

:open:

%
% Solve the BVP of the form
%
%   y'' + p y' + q y = r
%   y(a) = ya, y(b)=yb
%
% with N subinterval.
%
%   h    = (b-a)/N
%   x(k) = a + h*k,   
%   y(k) ~ y(a+h*k),
%
%  p : reference of a function of one argument p(x)
%  q : reference of a function of one argument q(x)
%  r : reference of a function of one argument r(x)
%
%  BCa : vector of two element [ a , y(a) ] with left BC
%  BCb : vector of two element [ b , y(b) ] with ight BC
%
%  N : number of subdivision of interval [a,b]
%
function [x,y] = BVPsolver( p, q, r, BCa, BCb, N )
  % extract the boundary conditions
  a  = BCa(1) ;
  b  = BCb(1) ;
  ya = BCa(2) ;
  yb = BCb(2) ;
  % fill x
  h  = (b-a)/N ;
  x  = a + h*[0:N]' ;
  % fill the 3 diagonal of the linear system and the RHS
  alpha = zeros( N-1, 1 ) ;
  beta  = zeros( N-1, 1 ) ;
  gamma = zeros( N-1, 1 ) ;
  omega = zeros( N-1, 1 ) ;
  
  alpha = -2/h^2 + feval( q, x(2:end-1) ) ;
  beta  = 1/h^2 - feval( p, x(2:end-1) ) ./(2*h) ;
  gamma = 1/h^2 + feval( p, x(2:end-1) ) ./(2*h) ;
  omega = feval( r, x(2:end-1) ) ;
  omega(1)   = omega(1)   - beta(1)*ya ;
  omega(N-1) = omega(N-1) - gamma(N-1)*yb ;

  % this is equivalent to the following loop
  % for k=1:N-1
  %   alpha(k) = -2/h^2 + feval( q, x(k+1) ) ;
  %   beta(k)  = 1/h^2 - feval( p, x(k+1) ) /(2*h) ;
  %   gamma(k) = 1/h^2 + feval( q, x(k+1) ) /(2*h) ;
  %   omega(k) = feval( r, x(2:end-1) ) ;
  % end
  % omega(1)   = omega(1)   - beta(1)*ya ;
  % omega(N-1) = omega(N-1) - gamma(N-1)*yb ;
  
  % fill the tridiagonal system
  A = zeros( N-1, N-1 ) ;
  A = diag(alpha,0) + diag(gamma(1:N-2),1) + diag( beta(2:N-1),-1) ;

  % this is equivalent to the following loop
  % for k=1:N-1
  %   A(k,k) = alpha(k) ;
  % end  
  % for k=1:N-2
  %   A(k+1,k) = beta(k+1) ; % fill the lower diagonal
  %   A(k,k+1) = gamma(k)  ; % fill the upper diagonal
  % end
  
  % fill the RHS
  b = zeros( N-1, 1 ) ;
  b = omega ;

  % solve the linear system
  sol = A\b ;
  
  % fill the y
  y = [ ya; sol; yb] ;

end

A simple BVP defined by the functions p1, q1, r1 to test the code: File p1.m

:open:

function res = p1( x )
  res = -x ;
end

File q1.m

:open:

function res = q1(x)
  %res = -1 ;
  res = -ones(size(x)) ;
end

File r1.m

:open:

function res = r1(x)
  %res = 1 ;
  res = ones(size(x)) ;
end

Exact solution to check the order: File exact1.m

:open:

function res = exact1(x)
  t1  = exp(0.5);
  t3  = sqrt(2);
  t5  = erf(t3 / 2);
  t10 = erf(x .* t3 / 2);
  t11 = x .^ 2;
  t13 = exp(t11 / 2 - 1 / 2);
  t17 = exp(t11 / 2);
  res = t13 .* t10 ./ t5 .* (3 - t1) + t17 - 1 ;
end

A simple script to check the order: File checkerror.m

:open:

BCa = [0;0] ;
BCb = [1;2] ;

N = [ 10, 20, 40, 80, 160, 320 ] ;

for i=1:length(N)
  [x,y] = BVPsolver( @p1, @q1, @r1, BCa, BCb, N(i) );
  err(i) = norm( y - exact1( x ), inf ) ;
end

for i=2:length(N)
  h1 = 1/N(i-1) ;
  h  = 1/N(i) ;
  p  = log( err(i-1)/err(i) ) / log( h1/h ) ;
  disp(p) ;
end

---