Python 数学

Pythonで方程式を解く方法（SciPy、ニュートン法、二分法による計算）

本記事の目的と対象

本記事の目的は、Pythonで方程式を解くプログラムを作成することである。

SciPyモジュールの関数を使って解く方法と、方程式の数値解法であるニュートン法、および二分法を素直にコーディングする方法の3種類を提示する。

本記事では例として、以下の方程式の解を求める。

この方程式は因数分解などの単純な方法では解くことが出来ない。 \begin{eqnarray*} x^2-5x+2=e^x \end{eqnarray*}
対象者は以下の条件を満たす読者である。
  • Pythonを書ける環境が整っている人
  • 高校レベルの数学（微分）を理解している人
  • Pythonで方程式を解く方法を知りたい人
また、以下の記事の内容を前提とする。

なお、本記事の内容は下記書籍の内容を参考にしているため、合わせて参照してほしい。
みんなのPython 第4版
posted with カエレバ
柴田 淳 SBクリエイティブ 2016-12-22
Amazonで購入
楽天市場で購入
Yahooショッピングで購入
7netで購入
モンテカルロ法によるリアル・オプション分析―事業計画の戦略的評価(CD-ROM付き)
posted with カエレバ
大野 薫 きんざい 2012-12-25
Amazonで購入
楽天市場で購入
Yahooショッピングで購入
7netで購入

目次

Python/SciPy.Otptimizeを用いる方法

SciPyは科学技術計算に用いられるPythonライブラリの一つである。

SciPyには最適化計算（関数の最大値や等式の数値解を求めること）に用いられるモジュールscipy.optimizeが含まれており、これを利用して方程式を解くことが出来る。

まず、必要なモジュールをインポートする。

方程式に含まれる指数関数も、同時に呼び出している*1。
from scipy import optimize, exp

次に解きたい方程式を定義する。

等式\( x^2-5x+2=e^x\)の左辺を\( f(x)\)、右辺を\( g(x)\)とおき、その差\( h(x)=f(x)-g(x)\)が\(0 \)に等しいという方程式を解く。
def f(x):
    return x**2-5*x+2
def g(x):
    return exp(x)
def h(x):
    return f(x)-g(x)

定義した関数を引数として、SciPyのoptimize.fsolve関数を使って方程式を解く。 
>>> optimize.fsolve(h,0) #0近辺から解を探索
array([ 0.16890436])

こうして、方程式 \begin{eqnarray*} x^2-5x+2&=&e^x\\ \Leftrightarrow f(x)&=&g(x)\\ \Leftrightarrow h(x)&=&0 \end{eqnarray*} の解は\( x=0.16890436\)であることがわかった。

ニュートン法

次にニュートン法で方程式を解く。

ニュートン法とは、漸化式 \begin{eqnarray*} x_{n+1}=x_n-\frac{h(x_n)}{h'(x_n)} \end{eqnarray*} で表される数列が方程式の解に収束することを利用した方程式の数値解法である（詳細は[1]参照）。

まず、解きたい方程式に含まれる指数関数を呼び出す。
from math import exp

次に解きたい方程式を定義する。

等式\( x^2-5x+2=e^x\)の左辺を\( f(x)\)、右辺を\( g(x)\)とおき、その差\( h(x)=f(x)-g(x)\)が\( 0\)に等しいという方程式を解く。

ニュートン法は導関数\( h'(x)\)も必要になるため、これを（今回は頭の中で）計算し、定義しておく。
def f(x):
    return x**2-5*x+2
def g(x):
    return exp(x)
def h(x):
    return f(x)-g(x)
def h_prime(x):
    return 2*x-5-exp(x)

最後に、ニュートン法のプログラムを作成する。

ここでは解の探索の開始値であるiniと、誤差errを与えて\( h(x)=0\)を解くような関数として定義する。 
def Newton(ini,err):
    x_n=ini
    x_n_scc=0
    count=0
    while True:
        count=count+1
        x_n_scc=x_n-h(x_n)/h_prime(x_n)
        if h(x_n_scc)<err:
            break
        x_n=x_n_scc
    print(\
        "数値解は",x_n_scc,\
        "\nその時の関数の値は",h(x_n_scc),\
        "\n計算回数は",count,"です")

結果は次の通り。 
>>> Newton(0,1/100000)
数値解は 0.16890400780768014 
その時の関数の値は 2.04672386860949e-06 
計算回数は 2 です

ニュートン法による計算は収束が速く（つまり少ない回数で解にたどりつくことができ）、\( 2\)回の計算で数値解\( x=0.16890400\)を得た。

 scipy.optimizeによる計算とは若干異なっているが、これはどのくらいの精度になったら計算を停止するかという設定が異なっているためであると考えられる。

二分法（バイナリーサーチ法）

最後に、二分法（バイナリーサーチ法）で方程式を解く。

二分法とは、
"解を含む区間の中間点を求める操作を繰り返すことによって方程式を解く"
方法である[2]。

この方法は、まず最初に解が含まれるであろう区間を指定し、区間の中間点の前（正の方向）と後ろ（負の方向）のどちらに解が含まれるかを確かめる。


仮に中間点の前に解が含まれていれば、次に当該中間点を下限としたより小さな区間を指定する。

この区間に対して、同じように中間点を求め、その前後いずれに解が含まれるかを確かめる。

この作業を繰り返すことで、目的となる数値解を求めるのである（詳細は[3]参照）。

二分法による解の計算をするPythonのプログラムは、以下のようになる。

まず、解きたい方程式に含まれる指数関数を呼び出す。
from math import exp

次に解きたい方程式を定義する。

等式\( x^2-5x+2=e^x\)の左辺を\( f(x)\)、右辺を\( g(x)\)とおき、その差\( h(x)=f(x)-g(x)\)が\( 0\)に等しいという方程式を解くことになる。
def f(x):
    return x**2-5*x+2
def g(x):
    return exp(x)
def h(x):
    return f(x)-g(x)

最後に、二分法のプログラムを作成する。

ここでは解の探索区間の上限upperと下限lower、そして誤差errを与え、\( h(x)=0\)を解くような関数として定義する。 
def BinarySearch(upper,lower,err):
    up=upper
    low=lower
    count=0
    while True:
        count=count+1
        mid=(up+low)/2
        y=h(mid)
        if abs(y)<err: #解が発見された
            break
        elif h(low)*y<0: #解は下限と中間点の間にある
            up=mid
        else:                     #解は上限と中間点の間にある
            low=mid
    print(\
        "数値解は",mid,\
        "\nその時の関数の値は",y,\
        "\n計算回数は",count,"です")

結果は次の通り。 
>>> BinarySearch(1,-1,1/10000000)
数値解は 0.16890335083007812 
その時の関数の値は 5.88754549157855e-06 
計算回数は 19 です

二分法はニュートン法に比べ収束が遅く、\( 19\)回の計算で数値解\( x=0.16890335\)を得た。

前述の2つの方法とは、計算の精度と方法が異なるため、若干違った数値解が得られている。

まとめ

本記事では、Pythonで方程式を解く方法として、scipy.optimizeモジュールを使う方法、ニュートン法、そして二分法を示し、コードの例を示した。

これら方法のよれば、通常の方法では解くことのできない複雑な方程式であっても、数値計算によって解くことが出来る。

参考文献

[3]大野 薫, 2012, モンテカルロ法によるリアル・オプション分析―事業計画の戦略的評価, きんざい
モンテカルロ法によるリアル・オプション分析―事業計画の戦略的評価(CD-ROM付き)
posted with カエレバ
大野 薫 きんざい 2012-12-25
Amazonで購入
楽天市場で購入
Yahooショッピングで購入
7netで購入

また、下記を参考にした。

みんなのPython 第4版
posted with カエレバ
柴田 淳 SBクリエイティブ 2016-12-22
Amazonで購入
楽天市場で購入
Yahooショッピングで購入
7netで購入

*1: ここでmath.expを使うとエラーとなる。expはscipyの関数を使う必要がある。https://docs.scipy.org/doc/numpy/reference/generated/numpy.exp.html

0 件のコメント :

コメントを投稿