bvp_per.m

LU decomposition routines in matlab

%
%
% compute an approximate solution of the BVP
%
%     - alpha  y''(x) + beta y'(x) + gamma y(x) = f(x)  x in [a,b],   
%
%                         y'(a) = y'(b),
%
%                          y(a) = y(b),
%
% on [a,b] = [0,1]
%
% Matthias Heinkenschloss
% Feb 6, 2001
%

% set problem data
alpha = 1;
beta  = 0;
gamma = 1;



f1 = inline( ' ones(size(x)) ', 'x' );
f2 = inline( ' x.^2 ', 'x' );          % f(x) = x^2  (more interesting example)

% set mesh points
n   = 100;
dx  = 1/(n+1);
x   = [0:dx:1]';


% diagonal elements
i = (1:n)';
j = (1:n)';
s = [3*alpha/2+dx*dx; (2*alpha+dx*dx)*ones(n-2,1); 3*alpha/2+dx*dx];

% superdiagonal elements
i = [i; (1:n-1)'];
j = [j; (2:n)'];
s = [s; -(alpha-0.5*dx*beta)*ones(n-1,1)];

% subdiagonal elements
i = [i; (2:n)'];
j = [j; (1:n-1)'];
s = [s; -(alpha+0.5*dx*beta)*ones(n-1,1)];

i = [i; 1; n];
j = [j; n; 1];
s = [s; -0.5*(alpha+0.5*dx*beta); -0.5*(alpha-0.5*dx*beta)];

% system matrix
A = sparse( i, j, s, n, n);

figure(1)
subplot(2,2,1)
spy(A)
title('A')

[L,U] = lu(A);
subplot(2,2,3)
spy(L); title('L')
subplot(2,2,4)
spy(U); title('U')




disp(' Solve BVP with right hand side f(x) = 1')
disp(' Hit return to continue ....'); pause
% right hand side
b = dx*dx*f1(x(2:n+1));

y = U\(L\b);

figure(2)
plot((0:dx:1),[0.5*(y(1)+y(n)); y; 0.5*(y(1)+y(n))], '--')
xlabel('x')
%print -deps bvp-ex2.eps


disp(' Solve BVP with right hand side f(x) = x^2')
disp(' Hit return to continue ....'); pause
% right hand side
b = dx*dx*f2(x(2:n+1));

y = U\(L\b);

figure(2)
plot((0:dx:1),[0.5*(y(1)+y(n)); y; 0.5*(y(1)+y(n))], '--')
xlabel('x')
%print -deps bvp-ex2.eps