アルゴリズム入門ということで、、、
筆算（手計算）の参考のため、トレース出力しています。
!ガウスの消去法 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
|
新しく発言する
EXIT
インデックスへ
返事を書く 
ノートメニュー