Code covered by the BSD License

lbvp.m

Mark Holmes (view profile)

A simple solver for a 2nd order linear BVP

lbvp
```function lbvp

%  uses centered differences to solve the linear BVP
%       y'' + p(x)y' + q(x)y = f(x)   for xL < x < xR
% where
%      a0*y(xL) + b0*y'(xL) = c0   and  a1*y(xR) + b1*y'(xR) = c1

% it is required that a0 and b0 are not both zero (and same for a1 and b1)

% an explanation of how to use this code is given at the end of this file

% set boundary conditions
xL=0; a0=1; b0=0; c0=1;
xR=1; a1=1; b1=0; c1=2;

% specify number of grid points
nx=100;

x=linspace(xL,xR,nx);
h=x(2)-x(1);
% calculate the coefficients of finite difference equation
for i=2:nx-1
a(i)=-2+h*h*q(x(i));
b(i)=1-0.5*h*p(x(i));
c(i)=1+0.5*h*p(x(i));
f(i)=h*h*rhs(x(i));
end;
if b0==0
a(1)=1;  b(1)=0;  c(1)=0;  f(1)=c0/a0;
else
a(1)=-2+h*h*q(x(1))+2*h*a0*(1-0.5*h*p(x(1)))/b0;
b(1)=0;
c(1)=2;
f(1)=h*h*rhs(x(1))+2*h*c0*(1-0.5*h*p(x(1)))/b0;
end;
if b1==0
a(nx)=1;  b(nx)=0;  c(nx)=0;  f(nx)=c1/a1;
else
a(nx)=-2+h*h*q(x(nx)) - 2*h*a1*(1+0.5*h*p(x(nx)))/b1;
b(nx)=2;
c(nx)=0;
f(nx)=h*h*rhs(x(nx)) - 2*h*c1*(1+0.5*h*p(x(nx)))/b1;
end;

% solve the tri-diagonal matrix problem
y=tri(a,b,c,f);

% plot solution
plot(x,y)
hold on
box on
grid on
%axis([0 1.02 -2 2])

xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);

% define functions in BVP

function g=p(x)
g=100*(3*x-1);

function g=q(x)
g=100*x;

function g=rhs(x)
g=0;

% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
v(i-1) = c(i-1)/w;
w = a(i) - b(i)*v(i-1);
y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
y(j) = y(j) - v(j)*y(j+1);
end

%  This code is run by simply entering the command  lbvp  in MATLAB (making
%  sure that this file is in the path for MATLAB).  The user needs to first
%  enter certain information about the problem into this file before running it.

%  Needed Input From User

%  Boundary conditions (BCs):  xL, a0, b0, c0 (line 14)  and xR, a1, b1, c1 (line 15)
%  Number of grid points:  nx (line 18)
%
%  p(x):  this is the coefficient of y' (line 66)
%  q(x):  this is the coefficient of y (line 69)
%  rhs(x):  this is the right hand side of the equation (line 72)

%  Example:   y'' + 100(3x-1)y' + 100xy = 0  for 0 < x < 1  with y(0) = 1 and y(1) = 2
%    xL=0; a0=1; b0=0; c0=1;
%    xR=1; a1=1; b1=0; c1=2;
%    nx = 100
%    p = 100*(3x-1)
%    q = 100*x
%    rhs = 0

% What to do if it doesn't work (aside from checking on input errors)
%  1.  If the problem has a boundary or interior layer you should consider
%      increasing  nx  (don't be shy about making nx very large).  If the layer
%      is sharp then you might want to try  sbvp.m
%  2.  If q(x) is negative in the interval, then the possibility arises that
%      the solution is either non-unique or non-existent

%  Background
%  This is a simple code that solves linear 2nd order BVPs. It does not
%  use adaptive methods, but instead uses 2nd order centered differences on
%  a uniform grid.  The specific algorithm is described in the book Intro to Numerical
%  Methods for Differential Equations (Holmes, 2007) on pages 46-49.  It is very
%  robust, and was used to solve all of the linear BVPs in Intro to Perturbation
%  Methods (Holmes, 1995).  It is also relatively fast, and as an illustration
%  is solves the above example in 0.004 sec on a laptop.  If you replace the
%  100's in the DE with 10000's, and take nx = 10000, it takes only 0.48 sec.

%  It has been tested on MATLAB, version R2009b
%  April 8, 2010

```