Как пройти через несколько разложений LU, чтобы решить Ax = b с разными векторами b?

Question

В настоящее время я занимаюсь практикой проблем с MatLab, и я застрял в этой проблеме, которая включает в себя мост фермы Он просит меня прибавить вес в грузовике, припаркованном в узле 9, с шагом 0,1, начиная с веса w9 = 13000 и увеличивая до некоторого большого веса, пока мост не рухнет. Каждый луч может выдержать 20000 ньютонов силы (положительного или отрицательного). Следовательно, мост рухнет, когда какая-либо одна сила будет больше, чем 20000 по величине, или эквивалентно, когда максимум абсолютных значений всех сил больше или равен 20000.

Я использую цикл for для перебора весов W9 = 13000.0, 13000.1, 13000.2,. , , и так далее, чтобы решить Система Ax = b многократно для сил внутри моста, пока сила, по крайней мере, на одной из балок не превысит точку разрыва. Я использую разложение LU для эффективного решения системы.

A = [-0.5 1 0 0 0 0 0 0 0 0.5 0 0 0 0 0; 
-sqrt(3)/2 0 0 0 0 0 0 0 0 -sqrt(3)/2 0 0 0 0 0;
0 -1 1 0 0 0 0 0 0 0 -0.5 0.5 0 0 0; 
0 0 0 0 0 0 0 0 0 0 -sqrt(3)/2 -sqrt(3)/2 0 0 0;
0 0 -1 1 0 0 0 0 0 0 0 0 -0.5 0.5 0; 
0 0 0 0 0 0 0 0 0 0 0 0 -sqrt(3)/2 -sqrt(3)/2 0; 
0 0 0 -1 0.5 0 0 0 0 0 0 0 0 0 -0.5; 
0 0 0 0 -sqrt(3)/2 0 0 0 0 0 0 0 0 0 -sqrt(3)/2; 
0 0 0 0 -0.5 -1 0 0 0 0 0 0 0 0 -0.5; 
0 0 0 0 0 1 -1 0 0 0 0 0 0 -0.5 0.5; 
0 0 0 0 0 0 0 0 0 0 0 0 0 sqrt(3)/2 sqrt(3)/2; 
0 0 0 0 0 0 1 -1 0 0 0 -0.5 0.5 0 0; 
0 0 0 0 0 0 0 0 0 0 0 sqrt(3)/2 sqrt(3)/2 0 0; 
0 0 0 0 0 0 0 1 -1 -0.5 0.5 0 0 0 0; 
0 0 0 0 0 0 0 0 0 sqrt(3)/2 sqrt(3)/2 0 0 0 0];

w7 = 800;
w8 = 900;
w9 = 13000;
W = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; w7; 0; w8; 0; w9];

for kk = w9:0.1:20000
    for jj = 1:15
        [L,U,P] = lu(A);
         Y = L\(P*W);
         X = U\Y;
         if abs(X(jj,1)) > 20000
             break
         end
    end
end
save('A.dat','w9','-ascii')
save('B.dat','X','-ascii')

Когда я запустил его, w9 вообще не изменился, и у меня получился тот же вектор X. Я ожидал, что он будет повторяться некоторое время, прежде чем цикл действительно прервется, но в итоге он этого не сделает. Может кто-нибудь помочь мне решить эту проблему?

Answer 1

Вам необходимо добавить дополнительный «вес» в w9 или в любом другом месте, где это уместно, и переопределить W. Иначе в вашем цикле ничего не изменится ...

Я немного оптимизировал ваш код, так как вам не нужен внутренний цикл.

Наконец, вам нужно перейти к более высокому w9. Ваше состояние не будет выполнено до w9 = 22228.

w7 = 800;
w8 = 900;
w9_initial = 13000;

for w9 = w9_initial:0.1:25000  
    % Define W with the increased w9   
    W = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; w7; 0; w8; 0; w9];
    % Calculate. These calculations aren't dependent on jj, so no point having a loop
    [L,U,P] = lu(A);
    Y = L\(P*W);
    X = U\Y;
    if any( abs( X ) > 20000 )
         break
    end
end