数値解析と制御のためのScilab入門

数値解析と制御のためのScilab入門では、行列、微分方程式、制御などの基礎的な計算を紹介しています。また、計算結果をScilabを用いたグラフィック表示による可視化について紹介しています。

*

LU分解法で方程式を解く

   

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

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

LU分解法

LU分解法

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

SPONSORED LINK

3行3列の行列を解いてみることにします、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の値は下のようになります。

–>x =
1.
1.
– 1.

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

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

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


 - 数値計算, 方程式の解法, 行列計算

        

SPONSORED LINK

SPONSORED LINK

  関連記事

関数のオンライン定義deff
関数のオンライン定義 deff を使ってみる!

関数のオンライン定義を用いると、数式だけを呼び出して計算ができます。 関数 fu …

no image
行列の積|scilab入門

今回は行列の積について実際に計算してみました。 簡単な計算でも実際計算してみない …

scilab 絶対値
絶対値の計算 abs|scilab入門

絶対値の計算:abs の計算の紹介です。  5 – 9 の計算は   …

行列
行列演算の基礎計算 |scilab入門

行列演算はScilabの基本となっています。行列を用いることにより数式レベルで扱 …

行列
行列 3行3列の表示 |scilab入門

行列 Bが3行3列の場合  行列 が であるとします。これをscilabで表す場 …

gyoretu2-13 単位行列
ゼロ行列や定数行列を作る |scilab入門

ゼロ行列や定数行列を作ることにします。  ゼロ行列:zeros  定数行列:on …

scilab パイi
円周率 π を表す|scilab入門

円周率 π を表す場合には、パーセントの後にpiとし  %pi のように表記しま …

連立方程式
連立一次方程式を解いてみる|scilab入門

逆行列 inv と (¥)を用いて、簡単なの連立一次方程式を解いてみました。 次 …

多項式
多項式が有利式の場合の計算

フリーのシュミレーションソフトscilabを用いて色々な計算ができます。 常にバ …

行列
行列の対角要素を取り出す diag|scilab入門

  行列の対角要素を取り出す場合には、diagを用います。 行列の対角要素や角線 …