Transcript 連立方程式

連立一次方程式の解法（１）
連立二元一次方程式
𝑎11 𝑥1 + 𝑎12 𝑥2 = 𝑏1
𝑎21 𝑥1 + 𝑎22 𝑥2 = 𝑏2
クラメール（Cramer）の公式
𝑥1 = ∆1 ∆
𝑎11
∆= 𝑎
21
𝑎12
𝑎22
𝑏1
∆1 =
𝑏2
𝑥2 = ∆2 ∆
𝑎12
𝑎22
𝑎11
∆2 =
𝑎21
元の数が大きくなると行列式の計算は困難
𝑏1
𝑏2
連立一次方程式の解法（２）
連立n元一次方程式
𝑎11 𝑥1 + 𝑎12 𝑥2 + ⋯ + 𝑎1𝑛 𝑥𝑛 = 𝑏1
𝑎21 𝑥1 + 𝑎22 𝑥2 + ⋯ + 𝑎2𝑛 𝑥𝑛 = 𝑏2
⋮
𝑎𝑛1 𝑥1 + 𝑎𝑛2 𝑥2 + ⋯ + 𝑎𝑛𝑛 𝑥𝑛 = 𝑏𝑛
第一式をa11で割って、ak1倍して第k式から引く
第二式をa’22で割って、a’k2倍して第k式から引く
以降同様の処理を繰り返す
𝑥1 + 𝑎′12 𝑥2 + ⋯ + 𝑎′1𝑛 𝑥𝑛 = 𝑏′1
𝑥2 + ⋯ + 𝑎′2𝑛 𝑥𝑛 = 𝑏′2
⋮
𝑥𝑛 = 𝑏′𝑛
連立一次方程式の解法（３）
𝑥1 + 𝑎′12 𝑥2 + ⋯ + 𝑎′1𝑛 𝑥𝑛 = 𝑏′1
𝑥2 + ⋯ + 𝑎′2𝑛 𝑥𝑛 = 𝑏′2
⋮
𝑥𝑛 = 𝑏′𝑛
第n式をan-1,n倍して上の式（第n-1式）から引く
第n式をan-2,n倍し、第n-1式をan-2,n-1倍して
上の式（第n-2式）から引く
以降同様の処理を繰り返す
𝑥1 = 𝑏′′1
𝑥1 = 𝑏′′2
⋮
𝑥𝑛 = 𝑏′𝑛
ガウス（Gauss）の消去法
連立一次方程式の解法（４）
ガウスの消去法で連立方程式を解く外部関数： gauss.f
コマンドプロンプトで次のコマンドを入力
[s2………]$ cp ~takizawa/gauss.f .
半角スペースとピリオドを忘れずに!!
[s2………]$ ls
でファイルリストを表示、「gauss.f」がコピーされていることを確認
gauss.f の呼び出し利用するメインプログラム（例えば prog10.f）を自分で作成し、
[s2………]$ gfortran prog10.f gauss.f
としてプログラムをコンパイル
連立一次方程式の解法（５）
ガウスの消去法で連立方程式を解く外部関数： gauss.f
関数の引数
function gsnelm(mat, vec, dim, mxdim)
Integer::dim, mxdim, gsnelm
double precision::mat(mxdim,mxdim),vec(mxdim)
gsnelm → 解が求まったら「0」、解がない時は「0」以外の値を返す関数
dim
→ 連立方程式の数（元の数）
mxdim
→ メインプログラムで宣言された配列の数（mxdim ≧ dim）
mat
→ 連立方程式の係数（aij）行列を表す二次元配列
vec
→ 連立方程式の右辺（bk）を表す一次元配列
連立一次方程式の解法（６）
ガウスの消去法で連立方程式を解く外部関数： gauss.f
メインプログラムの例
連立三元方程式を解く場合
integer::dim,mxdm,gsnelm,rc
doubel precision::mat(5,5),vec(5)
dim=3
mxdm=5
mat(1,1)=1.0d+0
mat(1,2)=3.0d+0
・・・
mat(3,3)=2.0d+0
vec(1)=12.0d0
・・・
rc=gsnelm(mat,vec,dim,mxdm)
・・・
if (rc==0) then
print *,’x(1)=‘,vec(1)
print *,’x(2)=‘,vec(2)
・・・
endif
・・・
必ず一致させる！