数値解析と制御のための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

  関連記事

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

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

3次スプライン補間 :splin
3次スプライン補間 splin について!

計測したデータを補間して表示させたい場合には、3次スプライン補間 splinを用 …

行列の計算
行列の和の計算|scilab入門

行列の和の計算についてです。   2つの行列の和について考えます。例として3行3 …

複素行列の足算
複素行列の作り方|scilab入門

複素行列の作り方についてです。 複素数とは  a + bi 、 ここで a,b …

no image
逆行列 inv|scilab入門

逆行列は、方程式を解くときに使います。 逆行列 inv を使用して解く場合につい …

scilab 5.5.2 バージョンアップ
scilabが5.5.2 にバージョン アップになりました!

scilab は常にバージョンアップされています。 しかも無料です。 ちょっとシ …

線形システムの定義
多項式を用いた線形システムの定義 syslin

多項式による伝達関数や線形すステムの定義は、syslinを用いて定義できます。 …

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

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

行列の計算
行列の差の計算

行列の引き算についてです。 2つの行列の和について考えします。例として3行3列の …

scilab 階乗
階乗 prod|scilab入門

階乗の計算には prod を使用します。 例えば 1から10までの階乗を計算する …