daily pastebin goal
5%
SHARE
TWEET

Untitled

a guest Apr 26th, 2018 56 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import sys
  2. import numpy as np
  3. from sympy import diff, symbols, sin
  4.  
  5. ideal = (1 + (5 ** 0.5)) / 2
  6.  
  7.  
  8. def bisection(f, a, b, eps, i=0):
  9.     c = (a + b) / 2
  10.     dist = np.abs(a - b)
  11.     i += 1
  12.     if dist < 2 * eps:
  13.         return c, i
  14.     elif f(c) < f(c + eps):
  15.         return bisection(f, a, c, eps, i)
  16.     else:
  17.         return bisection(f, c, b, eps, i)
  18.  
  19.  
  20. def helper(b, a):
  21.     global ideal
  22.     return (a - b) / ideal
  23.  
  24.  
  25. def golden(f, x0, x1, eps, i=0):
  26.     alpha = x1 - helper(x0, x1)
  27.     betha = x0 + helper(x0, x1)
  28.     i += 1
  29.     if np.abs(x1 - x0) < eps:
  30.         return (x0 + x1) / 2, i
  31.     elif f(alpha) >= f(betha):
  32.         return golden(f, alpha, x1, eps, i)
  33.     else:
  34.         return golden(f, x0, betha, eps, i)
  35.  
  36.  
  37. def newton(eps=10 ** (-5)):
  38.     x, y = symbols('x y')
  39.     func = x ** 2 + 10 * (y - sin(x)) ** 2
  40.     f_x = diff(func, x)
  41.     f_y = diff(func, y)
  42.     f_xx = diff(func, x, x)
  43.     f_xy = diff(func, x, y)
  44.     f_yx = diff(func, y, x)
  45.     f_yy = diff(func, y, y)
  46.     val = [1, 1]
  47.  
  48.     def grad(a, b):
  49.         s = {x: a, y: b}
  50.         return np.array([f_x.evalf(subs=s), f_y.evalf(subs=s)],
  51.                         dtype=np.float64)
  52.  
  53.     def gesse(a, b):
  54.         s = {x: a, y: b}
  55.         return np.array([[f_xx.evalf(subs=s), f_xy.evalf(subs=s)],
  56.                          [f_yx.evalf(subs=s), f_yy.evalf(subs=s)]],
  57.                         dtype=np.float64)
  58.  
  59.     def next_iter():
  60.         p = np.linalg.solve(gesse(1, 1), grad(1, 1))
  61.         next_el = val - (np.linalg.solve(gesse(1, 1), grad(1, 1)))
  62.         norm = np.amax(next_el - p)
  63.         while norm < eps:
  64.             grad(next_el[0], next_el[1])
  65.             gesse(next_el[0], next_el[1])
  66.         return p
  67.     return next_iter()
  68.  
  69.  
  70. np.set_printoptions(precision=14)
  71. sys.setrecursionlimit(11000)
  72. # print(golden(math.sin, 0, 1, 0.01))
  73. # print(bisection(math.sin, 0, 1, 0.01))
  74. print(newton())
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top