iOS Pythonista, Python3 で非線形方程式を解く

はじめに

Pythonista は，iOS 用 Python (2.7) のプログラミング環境であり，numpy および matplotlib が同梱されており，大変便利なものですが，scipy を使うことができません．
よって，非線形方程式を解くには，基本に立ち戻ってコーディングしてやる必要があります．

これに対し，OS X や Windows 用 Python（私の場合は Python 3.4）では，scipy を使うことができるため，非線形方程式を scipy の機能を用いて簡単なコーディングにより解くことができます．

ここでは，Pythonista (scipy 無し) と Python + scipy による非線形方程式の解法の事例を紹介します．

なお，コーディングが簡単というだけであり，当然のことながら，scipyを用いた場合でも，初期値によっては必ずしも望むべき領域にある解が得られるとは限らないことには，注意が必要です．

Pythonista (iOS Python 2.7)での二分法

以下のプログラムは，越流堰において，流量$Q$が与えられた時の越流水深$x$を二分法により求めるものです．

解くべき方程式は，以下のとおり．$B$, $W$は定数，$Q$は既知の値，$C$は$x$の関数であり，$x$が求めるべき越流水深です．

$$
Q-C\cdot B\cdot x^{1.5}=0 \qquad (C=1.28+1.42\cdot x/W)
$$

from math import *

def func(x,Q):
    B=19.0
    W=2.453
    C=1.28+1.42*x/W
    return Q-C*B*x**(1.5)

def cal(x01,y02,Q):
    imax=100
    eps=1.0e-6
    x1=x01
    x2=x02
    for i in range(0,imax):
        xi=0.5*(x1+x2)
        f1=func(x1,Q)
        f2=func(x2,Q)
        fi=func(xi,Q)
        if abs(x2-x1) < eps: break
        if f1*fi < 0.0: x2=xi
        if fi*f2 < 0.0: x1=xi
        print '{0:5d}{1:15.7e}'.format(i,xi)
    return xi

#main routine
Q=30.0
x01=0.001
x02=10.000
x=cal(x01,x02,Q)
print '{0:15.7e}'.format(x)

Pythonista (iOS Python 2.7)での Newton-Raphson法

連立非線型方程式
$$
\begin{cases}
f(x,y)=0 \\
g(x,y)=0
\end{cases}
$$

を解くためには，以下の増分形連立方程式を逐次解き，収束点を求めます．

$$
\begin{equation}
\begin{bmatrix}
\cfrac{\partial f(x_i,y_i)}{\partial x} & \cfrac{\partial f(x_i,y_i)}{\partial y} \\
\cfrac{\partial g(x_i,y_i)}{\partial x} & \cfrac{\partial g(x_i,y_i)}{\partial y} \\
\end{bmatrix}
\begin{Bmatrix}
\Delta x \\
\Delta y
\end{Bmatrix}
=\begin{Bmatrix}
-f(x_i,y_i) \\
-g(x_i,y_i)
\end{Bmatrix}
\end{equation}
$$

$$
\begin{cases}
x_{i+1}=x_i+\Delta x \\
y_{i+1}=y_i+\Delta y
\end{cases}
$$

以下の事例は，２つの円の交点を求めるためのプログラムであり，関数 $f$ および $g$ は以下のものです．

$$
\begin{cases}
f(x,y)=(x-a_1)^2+(y-b_1)^2-r_1^2=0 \\
g(x,y)=(x-a_2)^2+(y-b_2)^2-r_2^2=0
\end{cases}
$$

from math import *

def func(x,y,a1,b1,a2,b2,r1,r2):
    f=(x-a1)**2+(y-b1)**2-r1**2
    g=(x-a2)**2+(y-b2)**2-r2**2
    fx=2*(x-a1)
    fy=2*(y-b1)
    gx=2*(x-a2)
    gy=2*(y-b2)
    return f,g,fx,fy,gx,gy

def newton(a1,b1,r1,a2,b2,r2,x0,y0):
    itmax=100
    eps=1.0e-6
    x=x0
    y=y0
    for i in range(0,itmax):
        f,g,fx,fy,gx,gy=func(x,y,a1,b1,a2,b2,r1,r2)
        coef=fx*gy-fy*gx
        dx=(-f*gy+g*fy)/coef
        dy=( f*gx-g*fx)/coef
        x=x+dx
        y=y+dy
        print '{0:5d}{1:15.7e}{2:15.7e}{3:15.7e}{4:15.7e}'.format(i,x,y,dx,dy)
        if abs(dx) < eps and abs(dy) < eps: break
    return x,y

#main routine
a1=0.000; b1= 5.000; r1=10.000
a2=0.000; b2=-5.000; r2= 3.000
x0=-5.0; y0=-0.1
x,y=newton(a1,b1,r1,a2,b2,r2,x0,y0)
print '{0:15.7e}{1:15.7e}'.format(x,y)

Scipyの利用 （Mac：Python 3.4）

Scipyを用いて，上記で扱った２つの問題を解いてみます．

堰を超える流れの水深を求める（Brent法）

# Water depth calculation
from scipy import optimize

def func1(x,Q):
    B=19.0
    W=2.453
    C=1.28+1.42*x/W
    return Q-C*B*x**(1.5)

#Main routine
#Solution of non-linear equation
Q=30.0
res=optimize.brentq(func1,0.001,10,args=(Q))
print('{0:15.7e}'.format(res))

２つの円の交点を求める（２元連立非線形方程式を解く）

# Intersection coordinates of 2 circles
from scipy import optimize

def func2(x,a1,b1,r1,a2,b2,r2):
    return [(x[0]-a1)**2+(x[1]-b1)**2-r1**2,
            (x[0]-a2)**2+(x[1]-b2)**2-r2**2]

#Solution of nonlinear simultaneous equation (2x2)
a1=0.000; b1= 5.000; r1=10.000
a2=0.000; b2=-5.000; r2= 3.000
sol=optimize.root(func2, [-5.0, -0.1], args=(a1,b1,r1,a2,b2,r2), method='broyden1')
print('{0:15.7e}{1:15.7e}'.format(sol.x[0],sol.x[1]))

Scipyによるコーディングは以下を参考にしました．

• "http://org-technology.com/posts/scipy-root-finding-scalar-functions.html"
• "http://org-technology.com/posts/scipy-root-finding-multidimensional-functions.html"

以上