function [x,t,u] = MOL_BEuler(d0,d1,f,G,T,h,m)
% Use the method of lines to solve
% u_t = u_xx + G(x,t), 0 < x < 1, 0 < t < T
% with boundary conditions
% u(0,t) = d0(t), u(1,t) = d1(t)
% and initial condition
% u(x,0) = f(x).
% Use the backward Euler's method to solve the system of
% ODEs. For the discretization, use a spatial stepsize of
% delta=1/m and a time step of h.
%
x = linspace(0,1,m+1)'; delta = 1/m; delta_sqr = delta^2;
t = (0:h:T)'; N = length(t);
% Initialize u.
u = zeros(m+1,N);
u(:,1) = f(x);
u(1,:) = d0(t); u(m+1,:) = d1(t);
% Create tridiagonal coefficient matrix.
a = -(h/delta_sqr)*ones(m-1,1); c = a;
b = (1+2*h/delta_sqr)*ones(m-1,1);
a(1) = 0; c(m-1) = 0; option = 0;
% Solve for u using the backward Euler's method.
for n=2:N
g = G(x(2:m),t(n));
g(1) = g(1) + (1/delta_sqr)*u(1,n);
g(m-1) = g(m-1) + (1/delta_sqr)*u(m+1,n);
f = u(2:m,n-1) + h*g;
switch option
case 0
[v,alpha,beta,message] = tridiag(a,b,c,f,m-1,option);
option = 1;
case 1
v = tridiag(alpha,beta,c,f,m-1,option);
end
u(2:m,n) = v;
end
u = u';
end % MOL_BEuler
function [x, alpha, beta, message] = tridiag(a,b,c,f,n,option)
% function [x, alpha, beta, message] = tridiag(a,b,c,f,n,option)
%
% Solve a tridiagonal linear system M*x=f
%
% INPUT:
% The order of the linear system is given as n.
% The subdiagonal, diagonal, and superdiagonal of M are given
% by the arrays a,b,c, respectively. More precisely,
% M(i,i-1) = a(i), i=2,...,n
% M(i,i) = b(i), i=1,...,n
% M(i,i+1) = c(i), i=1,...,n-1
% option=0 means that the original matrix M is given as
% specified above.
% option=1 means that the LU factorization of M is already
% known and is stored in a,b,c. This will have been
% accomplished by a previous call to this routine. In that
% case, the vectors alpha and beta should have been
% substituted for a and b in the calling sequence.
% All input values are unchanged on exit from the routine.
%
% OUTPUT:
% Upon exit, the LU factorization of M is already known and
% is stored in alpha,beta,c. The solution x is given as well.
% message=0 means the program was completed satisfactorily.
% message=1 means that a zero pivot element was encountered
% and thesolution process was abandoned. This case happens
% only when option=0.
if option == 0
alpha = a; beta = b;
alpha(1) = 0;
% Compute LU factorization of matrix M.
for j=2:n
if beta(j-1) == 0
message = 1; return
end
alpha(j) = alpha(j)/beta(j-1);
beta(j) = beta(j) - alpha(j)*c(j-1);
end
if beta(n) == 0
message = 1; return
end
end
% Compute solution x to M*x = f using LU factorization of M.
% Do forward substitution to solve lower triangular system.
if option == 1
alpha = a; beta = b;
end
x = f; message = 0;
for j=2:n
x(j) = x(j) - alpha(j)*x(j-1);
end
% Do backward substitution to solve upper triangular system.
x(n) = x(n)/beta(n);
for j=n-1:-1:1
x(j) = (x(j) - c(j)*x(j+1))/beta(j);
end
end % tridiag