NURBS: Где я могу найти эти две служебные функции линейной алгебры?

Я работаю над Книгой NURBS Пигла и Тиллера. Для алгоритма глобальной интерполяции они требуют, чтобы вы предоставили две утилиты для решения системы линейных уравнений:

LUDecomposition(A, q, sbw) разложить матрицу коэффициентов q x q с полушириной sbw на нижнюю и верхнюю треугольные составляющие; для простоты мы предполагаем, что A представляет собой массив из q x q квадратов, но следует использовать утилиту, которая хранит только ненулевую полосу.

ForwardBackward(A, q, sbw, rhs, sol) для выполнения прямой/обратной замены (см. [Press88]); rhs[] — правая часть системы (координаты Q_k), а sol[] — вектор решения (координаты P_i).

Проверив ссылку Press88, я обнаружил, что это Числовые рецепты на C. Я должен быть в состоянии переработать алгоритм из этой книги, чтобы получить функцию ForwardBackward, но что касается LUDecomposition, где я могу найти алгоритм, который работает для особого случая матрицы с диагональными полосами?


linear-algebra matrix-decomposition
person nonremovable    schedule 16.05.2018    source источник
comment
Просто чтобы уточнить, A предполагается tridiagonal?   -  person StaticBeagle    schedule 29.05.2018
comment
Я считаю, что они диагональны, но не обязательно трехдиагональны. Отсюда и условие полуширины.   -  person nonremovable    schedule 29.05.2018
comment
Да, вы правы, спасибо за разъяснение.   -  person StaticBeagle    schedule 29.05.2018

Ответы (1)


arrow_upward
0
arrow_downward

Код Matlab для разложения LU трехдиагональной матрицы

 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

Чтобы решить, передний сабвуфер и задний сабвуфер.

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
person Community    schedule 01.06.2018