方程式の数値解析編

作成:2000-03-29

更新:2000-04-20

連立方程式高次方程式非線型方程式の数値解析


// 連立方程式

//==================================================
/*
次の連立方程式を解く。

 
 x + 2y -  z =  2
3y + 4z = 18
x - z = -2

これを

(1 2 -1 )(x) ( 2)
(0 3 4 )(y) = (18)
(1 0 -1 )(z) (-2)

の MX=A  という形にしておけば
X=inv(M)*A (inv(M)はMの逆行列)
で解ける。

(以下MaTXのプロンプトは>>で示す。)
*/

>> M=[[1 , 2 , -1]
[0 , 3 , 4]
[1 , 0 , -1]]

=== [M] : ( 3, 3) ===
( 1) ( 2) ( 3)
( 1) 1.00000000E+000 2.00000000E+000-1.00000000E+000
( 2) 0.00000000E+000 3.00000000E+000 4.00000000E+000
( 3) 1.00000000E+000 0.00000000E+000-1.00000000E+000

>> A=[[ 2]
[18]
[-2]]

=== [A] : ( 3, 1) ===
( 1)
( 1) 2.00000000E+000
( 2) 1.80000000E+001
( 3)-2.00000000E+000

// Mの逆行列の記述法には何種類かあるようだ。

>> M~*A

=== [ans] : ( 3, 1) ===
( 1)
( 1) 1.00000000E+000
( 2) 2.00000000E+000
( 3) 3.00000000E+000

>> inv(M)*A

=== [ans] : ( 3, 1) ===
( 1)
( 1) 1.00000000E+000
( 2) 2.00000000E+000
( 3) 3.00000000E+000

>> M^-1 * A

=== [ans] : ( 3, 1) ===
( 1)
( 1) 1.00000000E+000
( 2) 2.00000000E+000
( 3) 3.00000000E+000

//いずれにせよ、x=1,y=2,z=3が求まる。



// 高次方程式の解法

>> x=Polynomial("x");

>> y=x^3+x^2+x^1+1

y = x^3 + x^2 + x + 1

// roots(y)でy=0の解を求める。 round2z()は誤差を丸める関数
>> round2z(roots(y))

=== [ans] : ( 3, 1) ===
[ ( 1)-Real ( 1)-Imag ]
( 1) 0.00000000E+000 1.00000000E+000
( 2) 0.00000000E+000-1.00000000E+000
( 3)-1.00000000E+000 0.00000000E+000

//答えは (i , -i , -1 )


一般の(非線型)方程式の解 //==================================================
// fsolvedemo.mm 一般の(非線型)方程式の解
// 注意:MaTXのバージョンが5.08未満だとfsolveを使うことは出来ないので
// http://www.matx.org/ にあるパッチを当てること。
// x-cos(x)=0 を解いてみる。
// まずグラフを書いてみる。
>> x=[0:0.01:1];
>> y=x-cos(x);
>> gplot(x,y); gplot_cmd("set grid") ; gplot_cmd("replot")

image/eqation.png

//なんとなくx=0.7付近でy=0になるので、このに答えがあるような。
// minimum でabs(y)の中から一番小さい(0に近い)物をあぶりだす
>> minimum(abs(y))
ans = {0.00153144, 75}

//というわけで、x=0.75あたりかな? (y=0.00153144)

//fsolveで収束計算させてみる。
//fsolveの引数はMatrixなので、Matrixな関数を定義しなければならない。
//

>> Func Matrix ff(x)
Matrix x;
{
Matrix y;
y(1)=x(1)-cos(x(1));
return y;
}
// //さっき求めた0.7を初期値にして解いてみる。xfに解が入り、termcodeに収束コードが入る。
>> x0=[0.7];
>> {xf, termcode, path} = fsolve(ff,x0)
ans = {MATRIX, 1, MATRIX}

// termcode = 1 なので収束したらしい。答えはxfに入っている。

>> print xf
=== [xf] : ( 1, 1) === ( 1)
( 1) 7.39082039E-001
//答えは 0.73908203