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

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


  ガウスの消去法 連立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 !ピボットが0なら
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 
  │└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 

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