Search code examples
matlablinear-algebranumericalspline

Methods for Linear system solution with matlab

I have a linear system Ax = b , which is created by natural splines and looks like this:

enter image description here

where enter image description here

The code in matlab which is supposed to solve the system is the following: 

 clear; 
 clc;

x = [...] ;
a = [...];
x0 = ...;
n = length(x) - 1 ; 

for i = 0 : (n-1) 
    h(i+1) = x(i+2) - x(i+1) ; 
end


b= zeros( n+1 , 1 ) ;
for i =2: n 
    b(i,1) = 3 *(a(i+1)-a(i))/h(i) - 3/h(i-1)*(a(i) - a(i-1) ) ;
end


%linear system solution. 

l(1) =0 ; m(1) = 0 ; z(1) = 0 ;

for i =1:(n-1) 
    l(i+1) = 2*( x(i+2) - x(i) ) - h(i)* m(i) ;
    m(i+1) = h(i+1)/l(i+1);
    z(i+1) = ( b(i+1) - h(i)*z(i) ) / l ( i+1) ;
end

l(n+1) =1;
z(n+1) = 0 ;
c(n+1) = 0 ;

for j = ( n-1) : (-1) : 0 
    c(j+1) = z(j+1) - m(j+1)*c(j+2) ;
end

but I can't understand which method is being used for solving the linear system. If I had to guess I would say that the LU method is used, adjusted for tridiagonal matrices, but I still can't find the connection with the code...

Any help would be appreciated!!!

Solution

  • The coefficients look a little odd (particularly that 2 in the l equation), but it looks like a specialized Thomas Algorithm where:

    1. The second-to-last loop performs a forward elimination of the subdiagonal to bring the matrix into upper triangular form.
    2. The last loop performs the back substitution for the solution.

    The code doesn't seem to match one-to-one with the general algorithm since the solution is using the vectors that compose the diagonals instead of the diagonals themselves with no apparent preallocation of memory. So I can't say if this method is "better" than the general one off the bat.