1

I am working through The NURBS Book by Piegl and Tiller. For the global interpolation algorithm, they require you to provide two utility routines for solving a system of linear equations:

LUDecomposition(A, q, sbw) to decompose the q x q coefficient matrix with semibandwidth sbw into lower and upper triangular components; for simplicity we assume A is an q x q square array, but a utility shoule be used which only stores the nonzero band.

ForwardBackward(A, q, sbw, rhs, sol) to perform the forward/backward substitution (see [Press88]); rhs[] is the right hand side of the system (the coordinates of the Q_k), and sol[] is the solution vector (coordinates of the P_i).

Checking the reference Press88, I found that it is Numerical Recipes in C. I should be able to rework the algorithm in that book to get the ForwardBackward function, but as far as the LUDecomposition goes, where can I find one that works for the special case of a matrix with diagonal bands?

nonremovable
  • 798
  • 1
  • 6
  • 19

1 Answers1

0

Matlab code for LU decomp of Tridiagonal matrix

 function [u1,d1,l1] = decomt(u,d,l)
    %length of diagonal
    n=length(d);
    u1=u;
    d1 = d;
    l1=l;

    %perform LU decomp

    d1(1) = d(1);
    for i =2:n
        l1(i-1) = l(i-1)/d1(i-1); %update lower triangle
        d1(i)= d(i) - (l(i-1)/d1(i-1))*u(i-1); % update diagonal
    end
    end

In order to solve, front sub and back sub.

function [x] = solvet(u,d,l,b) 

n=length(d);
x = (1:n);
y =(1:n);

% Solve tridiag LUx=b

% Step 1 Solve Ly=b for y

y(1) = b(1)

for i=2:n
    y(i) = b(i) - l(i-1)*y(i-1);
end

%Step 2 : Solve Ux=y for x
x(n) = y(n)/d(n);
for i=(n-1):-1:1
    x(i) = (y(i)-u(i)*x(i+1))/d(i);
end
end