In [1]:
import sympy as sp
from sympy.plotting import plot
%matplotlib inline
sp.init_printing()
<img src="http://civil.engr.siu.edu/cheval/engr351/Images/ENGR351.jpg" width="500px" height="300px" >
Задано функцію $f(x)$, потрібно знайти корінь цієї функції, тобто хоча б одне значення параметру $x=x_0$, при якому $f(x_0)=0$. Якщо такого значення не існує повернути $null$.
Розглянемо три різні методи розвязку даної задачі:
Кожен з цих методів має свої недоліки і переваги, тому немає однозначно найкращого методу для розвязання цїєї задачі.
Для початку введемо декілька загальнопринятих позначень: $\epsilon$ та $x$ як символи бібліотеки SymPy
In [54]:
EPS = sp.Rational("1e-3")
x = sp.Symbol("x")
Визначимо функцію $fun$, для якої ми збираємося шукати корінь
In [53]:
fun = x * x * x - 2 * x
plot(fun, (x, -2, 2))
Out[53]:
Та її похідну $der$, що необхідна для коректної роботи деяких методів
In [55]:
der = sp.diff(fun, x)
plot(der, (x, -2, 2))
Out[55]:
Метод полягає у зменшені відрузку що розглядається вдвічі на кожній ітерації. Необхідна умова для застосування цього метода $f(a) \cdot f(b) <= 0$
Покладемо $l = a, r = b$, тоді виконується інваріант $f(l) \cdot f(r) <=0$. Покажемо що він зберігається на кожній ітерації.
На кожній ітерації циклу вибирається точка $m = \large\frac{l + r}{2}$, і перевіряється умова $f(a) \cdot f(m) <= 0$. Якщо вона виконується, тоді корінь знаходиться на проміжку $[a; m]$, інакше корінь треба шукати на проміжку $[m; b]$.
Рекурсивно виконуємо функцію пошуку для одного з вище вказаних проміжків.
In [51]:
def dih(a, b, f=fun, eps=EPS):
print("[{}; {}]".format(a, b))
if f.subs(x, a) * f.subs(x, b) > 0:
return None
if a > b:
a, b = b, a
if (b - a).evalf() <= EPS / sp.Integer(2):
return a
m = a + (b - a) / sp.Integer(2)
if f.subs(x, a) * f.subs(x, m) <= 0:
return dih(a, m, f, eps)
else:
return dih(m, b, f, eps)
In [52]:
res = dih(a=-5, b=sp.Rational('-0.1'))
"Result {}".format(sp.N(res))
Out[52]:
In [61]:
def newton(x0, f=fun, d=der, eps=EPS):
x1 = x0 - f.subs(x, x0) / d.subs(x, x0)
print(x1)
while sp.Abs(x1 - x0).evalf() > EPS / sp.Integer(2):
x0, x1 = x1, x1 - f.subs(x, x1) / d.subs(x, x1)
print(x1)
return x1
In [62]:
res = newton(x0=sp.Rational("0.7"))
"Result {}".format(sp.N(res, 10))
Out[62]:
In [63]:
alpha = sp.Symbol("alpha")
h = x - fun * alpha
h
Out[63]:
In [69]:
def simple(x0, alpha, f=fun, eps=EPS):
h = x - alpha * f
x1 = h.subs(x, x0)
print("[{}; {}]".format(x0, x1))
while abs(x1 - x0) > EPS / sp.Integer(2):
x0, x1 = x1, h.subs(x, x1)
print("[{}; {}]".format(x0, x1))
return x1
In [72]:
res = simple(x0=-3, alpha=1/10)
"Result {}".format(sp.N(res, 10))
Out[72]: