Прямая замена в упакованном массиве

Question

Итак, я пытаюсь создать скрипт, который возьмет массив, преобразует его в упакованную форму, затем запустит факторизацию Cholesky и, наконец, воспользуемся форвардной и обратной подстановкой для ее решения (Ax=b).

Проблема im im создает алгоритм прямой подстановки (столбца) для упакованных массивов. Я правильно выполнил обратную замену, но я не могу правильно преобразовать алгоритм.

Это текущий назад алгоритм замещения у меня есть (который идеально подходит)

x = b; 
p = n*(n+1)/2; 

for j = n:-1:1 
    x(j) = x(j)/u(p); 
    p = p - 1; 
    for i = j-1:-1:1 
    x(i) = x(i) - u(p) * x(j); 
    p = p - 1; 
    end 
end 

Где n является размер исходного A матрицы (которая является квадратным и СПД). Отображение из A матрицы к u выглядит следующим образом:

A(i,j) = u(i+j*(j-1)/2) 

Моя текущая итерация алгоритма прямого замещения заключается в следующем:

x = b; 
p = 1; 

for j = 1:n 
    x(j) = x(j)/u(p); 
    p = p + 1; 
    for i = j+1:n 
     x(i) = x(i) - u(p) * x(j); 
     p = p + 1; 
    end 
end 

Я не могу показаться, чтобы выяснить, что я делаю неправильно , Что бы я ни пытался, я просто давал мне NaN или Inf в качестве ответов для x. Может ли кто-нибудь умнее меня понять, как должен выглядеть алгоритм?

Answer 1

Насколько я понимаю, вы пытаетесь решить две проблемы, которые приводят к решению проблемы Ax = b.

Где A = LL' то есть A равна нижней треугольной матрицы (L), умноженная на его Транспонирование (L')

Во-первых, решая Ly = b с прямой замены.

И наконец, путем решения вопроса L'x = y с обратной заменой.

Где для спины substituion вашего u(p) определяется как таковые (при п = 3):

 | u(1) u(2) u(4) | 
L' = | 0 u(3) u(5) | 
    | 0 0 u(6) | 

вопрос я думаю, что вы бежите в с передней substition приходит образует то, что L будет определяться как :

| u(1) 0 0 | 
L = | u(2) u(3) 0 | 
    | u(4) u(5) u(6) | 

Ваше значение p нужно будет иметь следующую прогрессию, чтобы выполнить то, что вы хотите с передней substition:

{1, 2, 4, 3, 5, 6}

Я изменил ваш текущий код для прямого подзапроса, чтобы внедрить этот шаблон.

x = b; 
p = 1; 

for j = 1:n 
    d = j-1; 
    p = j + (d^2+d)/2; 
    x(j) = x(j)/u(p); 
    d = d + 1; 
    for i = j+1:n 
     p = j + (d^2+d)/2; 
     x(i) = x(i) - u(p) * x(j); 
     d = d + 1; 
    end 
end 

Я смог проверить это, и он дает ожидаемый результат. Надеюсь, это вам поможет.