ついにWATLABブログから書籍「いきなりプログラミングPython」が発売しました！
力学系の問題を式展開する場合は文字式のまま書き進める場合がほとんどです。連立方程式を数値的にコンピュータで解くアルゴリズムは多々ありますが、文字式の検算には使えません。ここではPython/Sympyを使って文字式をそのまま扱って連立方程式を解く方法を紹介します。
こんにちは。wat(@watlablog)です。ここではSympyを使って連立方程式を文字式のまま解く方法を紹介します！
コンピュータで連立方程式を文字式のまま解きたい場面
式展開の検算
これまで当WATLABブログでは、連立方程式を解く方法として以下の記事を紹介しました。
- Pythonで線形代数！連立1次方程式を解く
- Numpyだけで書いたガウスの消去法で連立1次方程式を解いてみた
- Pythonで連立方程式をヤコビ法（反復法）で解く方法
- Pythonで連立方程式をガウス・ザイデル法（反復法）で解く方法
- 連立方程式をSOR法で解くPythonコードと緩和係数のパラスタ
逆行列法やガウスの消去法は直接法、ヤコビ法/ガウス・ザイデル法/SOR法は反復法ですが、これらは全て数値解法として大別されます。
未知変数以外の全ての係数が既知（数で表現できる）である場合は上記の手法で連立方程式の答えを求めることができます。
しかしながら、係数行列の一部または全部が文字で表現された連立方程式を解かなければいけない場合、通常は人間が手計算を行います。
筆者のような計算が苦手なタイプの人はこの手計算工程で凡ミスをすることが多いため、常に間違っていないかどうかの検算を求めてしまいます。
このようなモチベーションで、本記事では文字式のまま連立方程式を解く方法を学び、コンピュータに先生になってもらうことを目指します。
Sympyとは？
今回はPythonの中でSympyというライブラリを使って数式を計算します。Sympyとは、Pythonで数式を記号計算できるライブラリです。
工学分野ではMathematicaというソフトが数式処理ソフトとして有名ですが、数式をそのまま処理できるソフトとして類似点は多いと思います。
当WATLABブログは過去にSympy関係で以下の記事を書いてきました。他にもありますので、ご興味のある方はサイト右上にある検索窓から探してみてください。
Pythonで連立方程式を文字式のまま解くコード
動作環境
この記事のプログラムは以下のPython環境で動作確認をしています。
| Python | Python 3.9.6 |
|---|---|
| PyCharm (IDE) | PyCharm CE 2020.1 |
| SymPy | 1.8 |
2元連立1次方程式を解くコード例
まずは式（1）の2元連立1次方程式を解いてみます。Sympyで式を表現する時は右辺を0にするように書きます。
以下がPythonコード例です。solveを使うだけなのでとても簡単です。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
from sympy import* if __name__ == '__main__': # 変数を定義 x, y, a, b, c, d, e, f = symbols('x, y, a, b, c, d, e, f') # 連立方程式を定義(右辺を0にした式として記述する） eq1 = a * x + b * y + e eq2 = c * x + d * y + f # 連立方程式を文字式のまま解く sol = solve([eq1, eq2], [x, y]) # 結果を表示 print(sol) |
sympy.latexを使えば、printで表示する時にLaTeXコードを出力可能です。論文を書く時はもちろん、ブログに式を貼る用途にもばっちりです。
|
1 |
print(latex(sol)) |
以下が結果です。文字式として解を出すことができました。
もちろん未知数x, y以外の値を数値として入力することで、直接数字としての解を求めることも可能です。
2元連立2次方程式を解くコード例
続いて式（2）の2元連立2次方程式を解いてみましょう。
以下がコード例です。記号と式を書き換えただけで、先ほどのコードと同じです。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
from sympy import* if __name__ == '__main__': # 変数を定義 x, y, a, b, c, d = symbols('x, y, a, b, c, d') # 連立方程式を定義(右辺を0にした式として記述する） eq1 = x + y eq2 = a * x ** 2 + b * y ** 2 + c * x * y + d # 連立方程式を文字式のまま解く sol = solve([eq1, eq2], [x, y]) # 結果を表示 print(sol) |
こちらが求められた解です。2次方程式なので2個ずつ解が得られました。
…ただこの程度だとしても式のボリュームがあるので本当か疑わしい。具体的に値を入れて答えを確認してみましょう。
式（2）の形は維持しつつ、式（3）の例題を作りました。
この式は代入法と因数分解により式（4）を得るため、
以下の解セットを得ます。これをSympyで確認しましょう。
こちらが全コードです。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
from sympy import* if __name__ == '__main__': # 変数を定義 x, y = symbols('x, y') # 連立方程式を定義(右辺を0にした式として記述する） eq1 = x + y - 1 eq2 = 2 * x ** 2 + y ** 2 + x * y - 16 # 連立方程式を文字式のまま解く sol = solve([eq1, eq2], [x, y]) # 結果を表示 print(sol) |
そして以下の結果を得ました。事前に出した解と同じ結果を得ましたね。
まとめ
本日は以上です。Short Tips的な記事となりましたが、Sympyを使って連立方程式を文字式のまま解くことができました。
例えば流体力学の問題では連続の式とベルヌーイの定理から得られた保存式を連立させることがありますが、そのような問題をプログラミングで解くことも可能となります。
それにしてもSympyは本当に便利ですね。
Sympyなら連立方程式も数式処理可能でした！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！
WATさま はじめまして．伊藤英樹です．
Python を勉強し始めたところですが，
代数方程式の係数が a,b,c…と文字
のままの場合の解法が分からず，悩んで
いました．
一気に理解が進みました．
ありがとうございました．
ご訪問ありがとうございます。
記事が参考になって本当に良かったです！
ご感想の投稿ありがとうございました。Sympyは非常に便利ですので、どんどんご活用いただければと思います。