2014年1月1日水曜日

[Python] 方程式を解く

Pythonで方程式を解くプログラム例.例として解く関数 f(x) は
    f(x) = x^2 - 2
とします.

まずは,2分法によって方程式を解くプログラム例(bisec.py).

# Global variables
a = 2           # f(x) = x * x - a
LIMIT = 1e-20   # Exit conditions

# Subcontract function
# f() function
def f(x):
    return x * x - a
# End of f() function

# Main execution unit
# Initial setting
xp = float(input("Input xp: "))
xn = float(input("Input xn: "))

# Iterative processing
while (xp - xn) * (xp - xn) > LIMIT:    # Repeat until the exit condition is satisfied.
    xmid = (xp + xn) / 2                # Calculating new intermediate values
    if f(xmid) > 0:                     # If the intermediate value is positive,
        xp = xmid                       # Update xp.
    else:                               # If the intermediate value is not positive,
        xn = xmid                       # Update xn
    print("{:.15f} {:.15f}".format(xn, xp))

実行例:
$ python bysec.py
Input xp: 1.5
Input xn: 1.3
1.400000000000000 1.500000000000000
1.400000000000000 1.450000000000000
1.400000000000000 1.425000000000000
1.412500000000000 1.425000000000000
1.412500000000000 1.418750000000000
1.412500000000000 1.415625000000000
1.414062500000000 1.415625000000000
1.414062500000000 1.414843750000000
1.414062500000000 1.414453125000000
1.414062500000000 1.414257812500000
1.414160156250000 1.414257812500000
1.414208984375000 1.414257812500000
1.414208984375000 1.414233398437500
1.414208984375000 1.414221191406250
1.414208984375000 1.414215087890625
1.414212036132813 1.414215087890625
1.414213562011719 1.414215087890625
1.414213562011719 1.414214324951172
1.414213562011719 1.414213943481445
1.414213562011719 1.414213752746582
1.414213562011719 1.414213657379150
1.414213562011719 1.414213609695435
1.414213562011719 1.414213585853576
1.414213562011719 1.414213573932648
1.414213562011719 1.414213567972183
1.414213562011719 1.414213564991951
1.414213562011719 1.414213563501835
1.414213562011719 1.414213562756777
1.414213562011719 1.414213562384248
1.414213562197983 1.414213562384248
1.414213562291116 1.414213562384248
$

上記のように,2分法で解いても良いのですが,ライブラリを活用する方法があります.
Pythonで正の平方根を求めるには,mathモジュールをインポートする必要があります.mathモジュールを用いると,xの正の平方根 √x a は以下のようになります.

# Import module
import math

# Main execution unit
# Input
x = float(input("Input the value you want to find a positive square root: "))

# Output
print ("sqrt(", x, ")=", math.sqrt(x))

実行例:
$ python sqrt.py
Input the value you want to find a positive square root: 2
sqrt( 2.0 )= 1.4142135623730951

Pythonでは,単に平方根を求めるだけではなく,方程式の解を求めるモジュールも用意されています.

# Import module
from sympy import *

# Main execution unit
var("x")                                        # Use variable x
equation = Eq(x**3 + 2 * x**2 - 5 * x -6, 0)    # Set equation
answer = solve(equation)                        # Solve the equation
print (answer)                                  # Output

実行例:
$ python solve.py
[-sqrt(2), sqrt(2)]
$

例として,別の方程式
    f(x) = x**3 + 2 * x**2 - 5 * x - 6
を解くプログラム例を以下に示します.

# Import module
from sympy import *

# Main execution unit
var("x")                                        # Use variable x
equation = Eq(x**3 + 2 * x**2 - 5 * x - 6, 0)   # Set equation
answer = solve(equation)                        # Solve the equation
print (answer)                                  # Output

実行例:
$ python solve2.py
[-3, -1, 2]

0 件のコメント :

コメントを投稿