##はじめに
こんにちは、麻菜結です。必要に迫られてしこしこ作ったプログラムが良い感じに作れたので紹介する記事です。作ったのはそこそこ実用的な連立一次方程式を解くプログラムです。バグがあったら優しく教えてください。なお、掃き出し方そのものの説明は行わないのでご了承ください。
##環境
環境はPython3系です。
> python --version
Python 3.9.0
##プログラム
プログラムは以下に示すものです。
import sys
def hakidashi(m, show_matrix):
num_line = len(m)
num_column = len(m[0])
div = lambda a, b: 0 if b==0 else a / b
for p in range(num_line):
m[p] = [div(n, m[p][p]) for n in m[p]]
for i in range(num_line):
x = div(m[i][p], m[p][p])
for j in range(num_column):
if(i != p):
m[i][j] = m[i][j] - x * m[p][j]
if(show_matrix):
for e in m:
print(e)
return list(map(lambda m: m[-1], m))
def dec_to_frac(f, max_denominator, calc_accuracy):
for i in range(1, max_denominator + 1):
k1 = int(f * i)
k2 = k1 + 1
if(abs((k1 / i) - f) < calc_accuracy):
return f"{k1}/{i}"
elif(abs((k2 / i) - f) < calc_accuracy):
return f"{k2}/{i}"
else:
return "not found"
def read_file(file_name, file_codec):
c, m = [], []
with open(file_name, encoding=file_codec) as f:
lists = list(f.readlines())
c = split_line(lists[0])
for l in lists[1:]:
m.append(split_line(l, num=True))
return c, m
def split_line(line, num=False):
line = line.strip()
line += "#"
data_line = []
word = ""
is_firstblank = True
for e in line:
if(" " == e and is_firstblank):
data_line.append(word)
is_firstblank = False
word = ""
elif("," == e):
data_line.append(word)
word = ""
elif("#" == e):
data_line.append(word)
break
else:
is_firstblank = True
word += e
else:
data_line.append(word)
return [
float(x.strip())
if num else
x.strip()
for x in data_line
if x not in ["", " "]
]
if __name__ == "__main__":
file_name = sys.argv[1]
file_encode = "utf-8"
max_denominator = 10000000
calc_accuracy = 0.0000000001
show_matrix = True
c, m = read_file(file_name, file_encode)
if(len(c) + 1 == len(m[0])):
ans = hakidashi(m, show_matrix)
for c, ans in zip(c, ans):
f = dec_to_frac(ans, max_denominator, calc_accuracy)
print(f"{c} = {ans} ({f})")
else:
print("data error", file=sys.stderr)
以下はそれぞれの説明です。
##掃き出し方
hakidashi()が該当関数です。特に工夫はありませんが、ゼロで割っちゃうエラーを回避するためにdivという関数を内部で定義しており、ゼロで割りそうなときは0を返すようにしてあります。リターンの直前にprintしているのは、もし解が無いような方程式を渡されてしまったとき、定数項が正しく1にならないのでそこで見て判断できます。
##分数表現
dec_to_frac()が該当関数です。分母を変えつつ検証していって、許せる精度になったら分数表現にした文字列を返します。もし範囲に見当たらなかったときはnot foundを表示します。ここがforループで計算回数が多いのであんまりに重かったらmax_denominatorやcalc_accuracyの値を調整してください。
こちらのサイトを参考にさせていただきました。
FloatToFraction (JavaScript版)
##ファイルの読み込み
read_file()とsplit_line()が該当関数です。
read_fileではコマンドラインからもらったファイル名からファイルデータを開き、split_lineに渡しています。cは文字ラベル(xとかyとか)を入れる配列、mは数値の配列データです。
split_lineは、カンマセパレートでもスペースセパレートでも読み込めるように、またはコメントを書き込めるように丁寧に分けるやつです。
##使い方
たとえば、以下のようなテキストを記述したファイルを使います。
x y z
4 2 1 26
1 6 3 48
4 1 5 34
上のは以下の連立方程式を表しています。
$$
4x + 2y + z = 26 \
x + 6y + 3z = 48 \
4x + y + 5z = 34
$$
> python solve_PE.py mat1.txt
[1.0, 0.0, 0.0, 2.727272727272727]
[0.0, 1.0, 0.0, 5.818181818181818]
[0.0, 0.0, 1.0, 3.454545454545455]
x = 2.727272727272727 (30/11)
y = 5.818181818181818 (64/11)
z = 3.454545454545455 (38/11)
また、以下のように解が無い場合を考えます。
x1, x2 # 下式が上式の
3, 6, 9 # ただの二倍になっているので
6, 12, 18 # 解はない
なお、#はコメントを表し、実際にファイルに記述することが可能です。
> python solve_PE.py mat2.txt
[1.0, 2.0, 3.0]
[0, 0, 0]
x1 = 3.0 (3/1)
x2 = 0 (0/1)
値は適当に入れられていますが、上の行列が三角っぽくないのでおかしいことがわかります。
##おわりに
お疲れさまでした。上手く動きましたか？百行も満たないのに結構使えるプログラムになったのでPythonすごいというのと、先人の考えたアルゴリズムすごいなという気持ちです（小学生並の感想）。やっぱ自分で組めると楽しいですね、こういうことに幸せを感じて生きていきたい。
最後まで読んでくださってありがとうございます。役に立ったなら幸いです。どうでもいいですけど、私の次は何か書こうかなというのは本当にあてになりませんね。自分で過去の記事を見てびっくりします。手がかじかんできたのでこの辺にします。ありがとうございました。