新しく発言する EXIT インデックスへ

ガウスの消去法連立1次方程式の解Ax=b

　　ガウスの消去法　連立1次方程式の解　Ax=b 山中和義 2007/11/25 10:21:25 
　　├!ガウス・ジョルダン（Gauss-Jordan）法N元... 山中和義 2007/11/25 10:22:50 
　　│└FullBASICなら逆行列による計算がお手軽です... 山中和義 2007/11/25 10:28:23 
　　├ヤコビ反復法（Jacobi） 山中和義 2008/09/29 10:28:16 
　　│└ガウス・ザイデル反復法（Gauss-Seidel） 山中和義 2008/09/29 10:32:12 
　　└!共役勾配法（conjugategradientmethod）N元... 山中和義 2008/10/06 20:58:41 

　　ガウスの消去法　連立1次方程式の解　Ax=b 山中和義 2007/11/25 10:21:25   ツリーへ
ガウスの消去法　連立1次方程式の解　Ax=b 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/11/25 10:21:25
アルゴリズム入門ということで、、、
筆算（手計算）の参考のため、トレース出力しています。



!ガウスの消去法　N元連立1次方程式の解　Ax=b

OPTION ARITHMETIC rational !有理数（分数）

!連立方程式を記述する
LET N=3 !変数（元）の数

DATA 2,3, 1, 4 !2*x1 +3*x2 +x3 = 4
DATA 4,1,-3,-2 !4*x1 +x2 -3*x3 = -2
DATA -1,2, 2, 2 ! -x1 +2*x2 +2*x3 = 2

DIM a(N,N+1) !係数のみ配列で記録する
MAT READ a


!前進消去により、i行のXiを消していく（掃き出し）。
!
!a11*x1+a12*x2+a13*x3=b1　−　(1)
!a21*x1+a22*x2+a23*x3=b2　−　(2)
!a31*x1+a32*x2+a33*x3=b3　−　(3)
!
!(2)-(1)*a21/a11、(3)-(1)*a31/a11より
!a11*x1+a12*x2+a13*x3=b1　−　(1)　←　ピボット
! a22'*x2+a23'*x3=b2'　−　(2)'
! a32'*x2+a33'*x3=b3'　−　(3)'
!
!(3)'-(2)'*a32'/a22'より
!a11*x1+a12*x2+a13*x3=b1　−　(1)
! a22'*x2+a23'*x3=b2'　−　(2)'　←　ピボット
! a33''*x3=b3''　−　(3)''
!

FOR k=1 TO N-1 !ピボットの選択　a(k,k)
LET mx=0
LET s=k
FOR i=k TO N !絶対値が最大となる行を探す
IF ABS(a(i,k))>mx THEN
LET mx=ABS(a(i,k))
LET s=i
END IF
NEXT i
IF mx=0 THEN !ピボットが０なら
PRINT "解けない。";k
STOP
END IF
IF s>k THEN
FOR j=1 TO N+1 !k行とs行を入れ替える
LET t=a(k,j)
LET a(k,j)=a(s,j)
LET a(s,j)=t
NEXT j
PRINT "(";STR$(k);")式と(";STR$(s);")式を入れ替える" !確認
MAT PRINT a !確認
END IF


FOR i=k+1 TO N !ビボット行以降の掃き出し
LET d=a(i,k)/a(k,k) !掃き出しのための係数
PRINT "(";STR$(i);")式 - (";STR$(k);")式 × ";d;"より、x";STR$(k);"を消去する。" !確認
FOR j=k TO N+1 !i行のXiを消去
LET a(i,j)=a(i,j)-a(k,j)*d !i行-k行×Aik/Akk
NEXT j
MAT PRINT a !確認
NEXT i
NEXT k


!後退代入により、Xiを順に求めていく。
!
!(3)''より、x3が求められる。
!　x3=b3''/a33''　−　(3)''の変形
!x3を(2)'に代入して、x2を求める。
!　x2=(b2'-a23'*x3)/a22'　−　(2)'の変形（x2以外は右辺へ）
!x3、x2を(1)に代入して、x1を求める。
!　x1=(b1-a12*x2-a13*x3)/a11　−　(1)の変形（x1以外は右辺へ）

FOR i=N TO 1 STEP -1
LET b=a(i,N+1)
FOR j=i+1 TO N !既に求めたXjを代入しながら
LET b=b-a(i,j)*a(j,N+1) !Xi以外は右辺へ（i行の変形）
NEXT j
LET a(i,N+1)=b/a(i,i) !係数で割って、Xiを得る

PRINT "x";STR$(i);"=";a(i,N+1) !結果を表示する
NEXT i


END



!参照　SAMPLEフォルダ内 invmat0.bas
 
　　├!ガウス・ジョルダン（Gauss-Jordan）法N元... 山中和義 2007/11/25 10:22:50   ツリーへ
Re: ガウスの消去法　連立1次方程式の解　Ax=b 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/11/25 10:22:50
!ガウス・ジョルダン（Gauss-Jordan）法　N元連立1次方程式（Linear Equations）の解　Ax=b

OPTION ARITHMETIC rational !有理数（分数）

!連立方程式を記述する
LET N=3 !変数（元）の数

DATA 2,3, 1, 4 !2*x1 +3*x2 +x3 = 4
DATA 4,1,-3,-2 !4*x1 +x2 -3*x3 = -2
DATA -1,2, 2, 2 ! -x1 +2*x2 +2*x3 = 2

DIM A(N,N+1) !係数のみ配列で記録する
MAT READ A


!a11*x1+a12*x2+a13*x3=b1　−　(1)
!a21*x1+a22*x2+a23*x3=b2　−　(2)
!a31*x1+a32*x2+a33*x3=b3　−　(3)
!
!(1)をa11で割る
!x1+a12'*x2+a13'*x3=b1'　−　(1)'　←　ピボット
!(2)-(1)*a21より
! a22'*x2+a23'*x3=b2'　−　(2)'
!(3)-(1)*a31より
! a32'*x2+a33'*x3=b3'　−　(3)'
!
!(2)'をa22'で割る
! x2+a23''*x3=b2''　−　(2)''　←　ピボット
!(1)'-(2)'*a12'より
!x1 +a13''*x3=b1''　−　(1)''
!(3)'-(2)'*a32'より
! a33''*x3=b3''　−　(3)''
!
!(3)''をa33''で割る
! x3=b3'''　−　(3)'''　←　ピボット
!(1)''-(3)''*a13''より
!x1 =b1'''　−　(1)'''
!(2)''-(3)''*a23''より
! x2 =b2'''　−　(2)'''
!
FOR k=1 TO N
LET mx=0
LET s=k
FOR i=k TO N !絶対値が最大となる行を探す
IF ABS(a(i,k))>mx THEN
LET mx=ABS(a(i,k))
LET s=i
END IF
NEXT i
IF mx=0 THEN !ピボットが０なら
PRINT "解けない。";k
STOP
END IF
IF s>k THEN
FOR j=1 TO N+1 !k行とs行を入れ替える
LET t=a(k,j)
LET a(k,j)=a(s,j)
LET a(s,j)=t
NEXT j
PRINT "(";STR$(k);")式と(";STR$(s);")式を入れ替える" !確認
MAT PRINT a !確認
END IF


LET p=a(k,k) !ピボットの選択
PRINT "(";STR$(k);")式 ÷ ";p;"より、x";STR$(k);" の係数を１にする。" !確認
FOR j=k TO N+1 !ピボット行をpで割る（Xkの係数を１にする）
LET a(k,j)=a(k,j)/p
NEXT j
MAT PRINT a !結果を表示する


FOR i=1 TO N !各行のXkを掃き出す
IF i<>k THEN !ピボット行は除く
LET d=a(i,k) !i行のXkの係数
PRINT "(";STR$(i);")式 - (";STR$(k);")式 × ";d;"より、x";STR$(k);" を消去する。" !確認
FOR j=k TO N+1 !Xkの係数を０にする
LET a(i,j)=a(i,j)-a(k,j)*d !i行-ピボット行*係数
NEXT j

MAT PRINT a !結果を表示する
END IF
NEXT i

PRINT
NEXT k

END
 
　　│└FullBASICなら逆行列による計算がお手軽です... 山中和義 2007/11/25 10:28:23   ツリーへ
Re: !ガウス・ジョルダン（Gauss-Jordan）法N元... 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/11/25 10:28:23
FullBASICなら逆行列による計算がお手軽です。
行列計算が可能ですから、記述が簡潔になります。


!逆行列による　N元連立1次方程式の解　Ax=b

OPTION ARITHMETIC rational !有理数（分数）

! 3*x1 -4*x2 +2*x3 = 0
! 2*x1 +5*x2 +3*x3 = 6
!-2*x1 +3*x2 -2*x3 = 1

!連立方程式を記述する
LET N=3 !変数（元）の数

DATA 3,-4, 2 !左辺
DATA 2, 5, 3
DATA -2, 3,-2

DIM A(N,N) !係数のみ配列で記録する
MAT READ A

DATA 0 !右辺
DATA 6
DATA 1
DIM b(N)
MAT READ b

DIM x(N) !解

!行列Aの逆行列をA'とする。Ax=bに左からかけると、A'Ax=A'b　∴x=A'b

DIM Ai(N,N) !逆行列を求める
MAT Ai=INV(A)
MAT PRINT Ai !確認

MAT x=Ai*b !左からかける

MAT PRINT x !結果を表示する

END
 
　　├ヤコビ反復法（Jacobi） 山中和義 2008/09/29 10:28:16   ツリーへ
Re: ガウスの消去法　連立1次方程式の解　Ax=b 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2008/09/29 10:28:16
ヤコビ反復法（Jacobi）



!ヤコビ反復法（Jacobi）によるN元連立1次方程式の解　Ax=b

!OPTION ARITHMETIC RATIONAL !有理数（分数）

LET cEps=1e-14 !誤差


!連立方程式を記述する
LET N=4 !変数（元）の数

DIM A(N,N) !左辺
DATA 10,-1,2,0 !係数のみ配列で記録する
DATA -1,11,-1,3
DATA 2,-1,10,-1
DATA 0,3,-1,8
MAT READ A

DIM b(N) !右辺
DATA 6,25,-11,15
MAT READ b


DIM x(N),x0(N)
MAT x0=ZER

LET cMax=100 !反復回数
FOR k=1 TO cMax

FOR i=1 TO N
LET s=0
FOR j=1 TO i-1 !s=Σ[j=1,N、i≠j]Aij*X0j
LET s=s+A(i,j)*x0(j)
NEXT j
FOR j=i+1 TO N
LET s=s+A(i,j)*x0(j)
NEXT j
LET x(i)=(b(i)-s)/A(i,i) !近似解 Xi=(bi-ΣAij*X0j)/Aii
NEXT i

LET s=0 !収束を確認する
FOR i=1 TO N
LET s=s+ABS(x(i)-x0(i)) !|x-x0|の和
NEXT i
IF s<cEps THEN EXIT FOR

MAT x0=x !次へ

NEXT k
IF k>cMax THEN
PRINT "収束しません。"
STOP
END IF

PRINT k;
MAT PRINT x; !結果を表示する



!検算
MAT x0=A*x
MAT PRINT x0; !b


END
 
　　│└ガウス・ザイデル反復法（Gauss-Seidel） 山中和義 2008/09/29 10:32:12   ツリーへ
Re: ヤコビ反復法（Jacobi） 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2008/09/29 10:32:12
ガウス・ザイデル反復法（Gauss-Seidel）


!ガウス・ザイデル反復法（Gauss-Seidel）によるN元連立1次方程式の解　Ax=b

!OPTION ARITHMETIC RATIONAL !有理数（分数）

LET cEps=1e-14 !誤差


!連立方程式を記述する
LET N=4 !変数（元）の数

DIM A(N,N) !左辺
DATA 10,-1, 2, 0 !係数のみ配列で記録する
DATA -1,11,-1, 3
DATA 2,-1,10,-1
DATA 0, 3,-1, 8
MAT READ A

DIM b(N) !右辺
DATA 6,25,-11,15
MAT READ b


DIM x(N),x0(N)
MAT x0=ZER

LET cMax=100 !反復回数
FOR k=1 TO cMax

FOR i=1 TO N
LET s=0
FOR j=1 TO i-1 !Σ[j=1,i-1]Aij*Xj
LET s=s+A(i,j)*x(j)
NEXT j
FOR j=i+1 TO N !Σ[j=i+1,N]Aij*X0j
LET s=s+A(i,j)*x0(j)
NEXT j
LET x(i)=(b(i)-s)/A(i,i) !近似解 Xi=(bi-Σ[j=1,i-1]Aij*Xj-Σ[j=i+1,N]Aij*X0j)/Aii
NEXT i

LET s=0 !収束を確認する
FOR i=1 TO N
LET s=s+ABS(x(i)-x0(i)) !|x-x0|の和
NEXT i
IF s<cEps THEN EXIT FOR

MAT x0=x !次へ

NEXT k
IF k>cMax THEN
PRINT "収束しません。"
STOP
END IF

PRINT k;
MAT PRINT x; !結果を表示する



!検算
MAT x0=A*x
MAT PRINT x0; !b


END
 
　　└!共役勾配法（conjugategradientmethod）N元... 山中和義 2008/10/06 20:58:41   ツリーへ
Re: ガウスの消去法　連立1次方程式の解　Ax=b 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2008/10/06 20:58:41
!共役勾配法（conjugate gradient method）　N元連立1次方程式の解　Ax=b

!2*x1 +3*x2 +x3 = 4
!4*x1 +x2 -3*x3 = -2
! -x1 +2*x2 +2*x3 = 2

!連立方程式を記述する
LET N=3 !変数（元）の数

DIM A(N,N) !左辺
DATA 2,3, 1 !係数のみ配列で記録する
DATA 4,1,-3
DATA -1,2, 2
MAT READ A

DIM b(N) !右辺
DATA 4,-2,2
MAT READ b


DIM x(N) !解
MAT x=ZER

DIM tA(N,N), c(N),p(N),q(N),r(N),s(N), T1(N)
MAT T1=A*x !c0=b-A*x0
MAT c=b-T1
MAT tA=TRN(A)
MAT r=tA*c
MAT p=r !p0=r0

LET c2=DOT(c,c) !c2=(c0,c0)
DO
LET c1=c2

MAT q=A*p !qk=A*pk
MAT s=tA*q !sk=tA*qk

LET alpha=DOT(p,r)/DOT(q,q) !αk=(pk,rk)/(qk,qk)
MAT T1=(alpha)*p !xk+1=xk + αk*pk
MAT x=x+T1
MAT T1=(alpha)*q !ck+1=ck - αk*qk
MAT c=c-T1
MAT T1=(alpha)*s !rk+1=rk - αk*sk
MAT r=r-T1

LET beta=-DOT(r,s)/DOT(q,q) !β=-(rk+1,sk)/(qk,qk)
MAT T1=(beta)*p !pk+1=rk+1 + βk*pk
MAT p=r+T1

LET c2=DOT(c,c) !c2=(c+1,ck+1)
LOOP UNTIL c1<=c2 !収束するまで

MAT PRINT x; !結果を表示する


END

インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。