方程式の数値解析編

作成:2000-10-02

更新:2000-10-03

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


// 連立方程式

次の連立方程式を解く。

 
 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の逆行列)
で解ける。

(以下SCILABのプロンプトは>>で示す。)
>> M=[[1 , 2 , -1]
[0 , 3 , 4]
[1 , 0 , -1]]

M = ! 1. 2. - 1. ! ! 0. 3. 4. ! ! 1. 0. - 1. !

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

A = ! 2. ! ! 18. ! ! - 2. !

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

>> inv(M)*A

ans = ! 1. ! ! 2. ! ! 3. ! >> M^-1 * A ans = ! 1. ! ! 2. ! ! 3. ! //いずれにせよ、x=1,y=2,z=3が求まる。


// 高次方程式の解法

>> x=poly(0,"x");

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

y = 2 3 1 + x + x + x

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

ans = ! - 5.441E-17 + i ! ! - 5.441E-17 - i ! ! - 1. !

>> clean(roots(y))

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


一般の(非線型)方程式の解 // 一般の(非線型)方程式の解

// x-cos(x)=0 を解いてみる。
// まずグラフを書いてみる。
>> x=[0:0.01:1];
>> y=x-cos(x);
>> plot(x,y)

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

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

//fsolveで収束計算させてみる。
//さっき求めた0.7を初期値にして解いてみる。xfに解が入り、termcodeに収束コードが入る。
>> deff('y=ff(x)' , 'y=x-cos(x)' ) >>fsolve(0.7,ff) ans = .7390851 //ちなみに、正しい書式(?)は下記の通りであり、infoは収束判定、vが最終的なxを使ったff(x)の値、xが最終値となっている。 >>[x,v,info]=fsolve(0,ff) info = 1. v = 0. x = .7390851 //答えは 0.7390851


return