Initial coefficients in Thomas's algorithm-CodePudding

x = 0.0:0.1:1.0;
N = length(x);
h = x(2)-x(1);
a = 1/hˆ2*ones(size(x));
b = -2/hˆ2*ones(size(x));
c = 1/hˆ2*ones(size(x));
r = x.*cos(x);
a(1) = -1; b(1) = 1; r(1) = 0;
c(N) = 0; b(N) = 1; r(N) = 0;
% Forward sweep
for j = 2:N
    b(j) = b(j)-c(j)*a(j-1)/b(j-1);
    r(j) = r(j)-c(j)*r(j-1)/b(j-1);
end
% Final equation
y(N) = r(N)/b(N);
for j = (N-1):-1:1
    y(j) = r(j)/b(j)-a(j)*y(j 1)/b(j);
end

This is the code of Thomas's algorithm for solving y''=xcosx with initial points y'(0)=0 and y(1)=0. I don't understand why a(1)=-1 and b(1)=1. Would someone explain it to me? I get confused.

CodePudding user response：

You are solving the tri-diagonal system

c(k)*y(k-1)   b(k)*y(k)   a(k)*y(k 1) = r(k),  k=2...N-1

The boundary conditions are implemented as

y(2)-y(1) = h*y_func'(0) = 0  ==> a(1)=1, b(1)=-1, r(1)=0 (and implicitly c(1)=0)

and

y(N)=y_func(1)=0 ==> c(N)=0, b(N)=1, r(N)=0 (and implicitly a(N)=0).

The algorithm computes an LU decomposition where both triangular matrices are bi-diagonal. The inversion of the L factor is applied in-place in the forward sweep, the back-substitution according to the r factor then in the reverse sweep.