LU分解法で方程式を解く

方程式を　LU分解法 (decomposition)　で解くことにします。

元連立一次方程式(system of linear equation)の係数行列 A の階数（rank）が n で、正則(nonsingular)である場合には、係数行列 A を下三角行列 L と上三角行列 U の積に分解できます。

連立一次方程式 分解を行う段階で対角要素、 が0になる場合や非常に小さな絶対値を持つことがあります。このような場合には演算不能や大きな誤差を生じる恐れがあるので、これを防止するためにピボット操作を行います。ピボット操作には、行または列を入れ替えて行う部分ピボットティングと、行と列を入れ替えて行う完全ボットティングがあります。ピボット操作については、他書を参考にして下さい。

３行３列の行列を解いてみることにします、Aは係数行列です。
Ax=b
の形式です。根本的には　x=inv(A)b　でxを求めることになります。
プログラマは以下のようになります。

1. //function [L,U] = ludecomp(A)
2. // LU decomposition
3. // A: matrix to decompose
4. // L: lower triangular matrix
5. // U: upper triangular matrix
6. A=[2 1 1 ; 2 3 1 ; 1 1 3]
7. B=[ 2 4 -1]’
8. [n,n] = size(A);
9. toler = 10^(-9);
10. L = zeros(n,n);
11. U = zeros(n,n);
12. for i=1:n-1
13. if abs(A(i,i))< toler 14. disp('The pivot is too small') 15. else 16. aii= 1 / A(i,i); 17. for j=i+1:n 18. aji = - A(j,i) * aii; 19. A(j,i+1:n) = A(j,i+1:n) + aji * A(i,i+1:n); // 下三角行列 20. end 21. A(i,i+1:n) = aii * A(i,i+1:n); // 上三角行列 22. end 23. end 24. for i=1:n-1 25. L(i,1:i) = A(i,1:i); 26. U(i,i+1:n) = A(i,i+1:n); U(i,i) = 1; 27. end 28. U(n,n) = 1; L(n,1:n) = A(n,1:n); 29. // end of lu decomp 30. y(1) = B(1)/L(1,1); 31. for i=2:n 32. B1 = 0; 33. for k=1:i-1 34. B1=B1+L(i,k)*y(k); 35. nd 36. y(i)= (B(i)-B1)/L(i,i); 37. end 38. x(n) = y(n) /U(n,n); 39. for i=n-1:-1:1 40. y1 = 0; 41. for k=i+1:n 42. y1=y1+U(i,k)*x(k); 43. end 44. x(i)= (y(i)-y1); 45. end 46. x

方程式の解ｘの値は下のようになります。

–>x =
1.
1.
– 1.

SciNotesに１～４６行までを記入して実行すると解が得られます。４６行はｘの値を表示しています。
行列AとBの数値を変えて実行すると、ある程度の大きさの方程式を解くことができます。本プログラムが、実際に実行したリストを載せておりますので会が得られるはずです。

解が得られない場合には、しっかりとディバックしてみて下さい。

方程式と解くための離散化するためのアルゴリズムは、割愛させていただきました。
他書を参考にしてプログラムと照らし合わせてみて下さい。