2016年1月1日金曜日

[Python] 方程式を解くプログラム

方程式(今回の例では f(x) = x^2 - a )を解くプログラム例を示します.

特別なモジュールを利用せずに,Python単体の言語機能のみで,2分法(Bisection method)アルゴリズムを使ったプログラム例:

# -*- coding: utf-8 -*-

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

def f(x):       # Function definition
    return x*x - a

# Main execution unit
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                # Computation of new intermediate value
    if f(xmid) > 0:
        xp = xmid                       # Update xp if intermediate value is positive
    else:
        xn = xmid                       # Update xn unless the intermediate value is positive

    print("{:.15f} {:.15f}".format(xn, xp))

実行例:
$ python bisec.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 モジュール)を用いて平方根を求めるプログラム例.

# -*- coding: utf-8 -*-

import math     # Import module

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

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

実行例:
$ python sqrt.py
Input a value to find a positive square root: 2
sqrt( 2.0 ) =  1.4142135623730951
$ python sqrt.py
Input a value to find a positive square root: 3
sqrt( 3.0 ) =  1.7320508075688772

Pythonモジュール(simply モジュール)を用いて方程式の解を求めるプログラム例:

# -*- coding: utf-8 -*-

from sympy import *     # Import module

var("x")    # Use variable x
equation = Eq(x**2 - 2, 0# Setting equation
answer = solve(equation)    # Solve equation

print(answer)               # Output result

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

上記の例だと答えが

と表示されています.試しに,

    x**3 + 2 * x**2 - 5 * x - 6

という方程式を解くプログラム例を以下に示します.
# -*- coding: utf-8 -*-

from sympy import *     # Import module

var("x")    # Use variable x
equation = Eq(x**3 + 2 * x**2 - 5 * x - 6, 0# Setting equation
answer = solve(equation)    # Solve equation

print(answer)               # Output result

実行例:
$ python solve2.py
[-3, -1, 2]
ちゃんと,-3, -1, 2 という答えが得られています.

0 件のコメント :

コメントを投稿