supporting recitation i did on wed 4/4

intro-num-methods.py

# newton raphson table generator

import math
def f(x):
    return x**2 + 2*x - 8

def df(x):
    return 2*x + 2

def newton_raphson(x0, f, df, n):
    x = x0
    for i in range(n):
        # print k, xk, f(xk), df(xk), f(xk)/df(xk)
        print(f"k = {i}, x = {x}, f(x) = {f(x)}, df(x) = {df(x)}, f(x)/df(x) = {f(x)/df(x)}")
        x = x - f(x)/df(x)


# newton_raphson(0, f, df, 10)

def foo(x1, x2):
  return x2 * math.cos(x1) + x2**2 * x1


# do symmetric difference approximation of f'(x_0) = (f(x_0 + h) - f(x_0 - h))/(2h)
def df_dx1(x1, x2, h):
  return (foo(x1+h, x2) - foo(x1-h, x2))/(2*h)

def df_dx2(x1, x2, h):
  return (foo(x1, x2+h) - foo(x1, x2-h))/(2*h)

def jacobian(x1, x2, h):
  return [[df_dx1(x1, x2, h), df_dx2(x1, x2, h)]]

print(jacobian(0, -1, 0.01))

print('fofofofo')

def poly(x):
  # 2x^4 - 3x^3 + 2x^2 + 4
  return 2*x**4 - 3*x**3 + 2*x**2 + 4

# newtons algorithm for scalar for minding x that minimizes f(x)
def newtons_method(x0, f, df, d2f, n):
  x = x0
  for i in range(n):
    if df(x) == 0:
      print("div by zero")
      break
    if d2f(x) == 0:
      print("div by zero")
      break
    print(f"k = {i}, x = {x}, f(x) = {f(x)}, df(x) = {df(x)}, df2(x) = {d2f(x)}, f(x)/df(x) = {f(x)/df(x)}")
    x = x - 1 * df(x)/d2f(x)

newtons_method(3, poly, lambda x: 8*x**3 - 9*x**2 + 4*x, lambda x: 24*x**2 - 18*x + 4, 10)

# hessian matrix using symmetric difference approximation
def d2f_dx1dx1(x1, x2, h):
  return (df_dx1(x1+h, x2, h) - df_dx1(x1-h, x2, h))/(2*h)

def d2f_dx1dx2(x1, x2, h):
  return (df_dx2(x1+h, x2, h) - df_dx2(x1-h, x2, h))/(2*h)

def d2f_dx2dx1(x1, x2, h):
  return (df_dx1(x1, x2+h, h) - df_dx1(x1, x2-h, h))/(2*h)

def d2f_dx2dx2(x1, x2, h):
  return (df_dx2(x1, x2+h, h) - df_dx2(x1, x2-h, h))/(2*h)

def hessian(x1, x2, h):
  return [[d2f_dx1dx1(x1, x2, h), d2f_dx1dx2(x1, x2, h)], [d2f_dx2dx1(x1, x2, h), d2f_dx2dx2(x1, x2, h)]]

print(hessian(math.pi/2, -1, 0.1))

# jacbian calc

# forward difference
def f1(x1, x2):
  # sinx1^7 + x2^3
  return math.cos(x1**7) + x2**2*x1

def f2(x1, x2):
  # 4x1*sqrt(x2)
  return 4*math.sqrt(x1)*x2**2

# forward difference approximation f'(x_0) = (f(x_0 + h) - f(x_0))/h
def df1_dx1(x1, x2, h):
  return (f1(x1+h, x2) - f1(x1, x2))/h

def df1_dx2(x1, x2, h):
  return (f1(x1, x2+h) - f1(x1, x2))/h

def df2_dx1(x1, x2, h):
  return (f2(x1+h, x2) - f2(x1, x2))/h

def df2_dx2(x1, x2, h):
  return (f2(x1, x2+h) - f2(x1, x2))/h

def jacobian_forward_diff(x1, x2, h):
  return [[df1_dx1(x1, x2, h), df1_dx2(x1, x2, h)], [df2_dx1(x1, x2, h), df2_dx2(x1, x2, h)]]

print(jacobian_forward_diff(1, 2, 0.001), " forward diff approx") 

# symmetric difference approximation f'(x_0) = (f(x_0 + h) - f(x_0 - h))/(2h)
def df1_dx1_sym(x1, x2, h):
  return (f1(x1+h, x2) - f1(x1-h, x2))/(2*h)

def df1_dx2_sym(x1, x2, h):
  return (f1(x1, x2+h) - f1(x1, x2-h))/(2*h)

def df2_dx1_sym(x1, x2, h):
  return (f2(x1+h, x2) - f2(x1-h, x2))/(2*h)

def df2_dx2_sym(x1, x2, h):
  return (f2(x1, x2+h) - f2(x1, x2-h))/(2*h)

def jacobian_sym_diff(x1, x2, h):
  return [[df1_dx1_sym(x1, x2, h), df1_dx2_sym(x1, x2, h)], [df2_dx1_sym(x1, x2, h), df2_dx2_sym(x1, x2, h)]]

print(jacobian_sym_diff(1, 2, 0.001), " symmetric diff approx")

# backward difference approximation f'(x_0) = (f(x_0) - f(x_0 - h))/h
def df1_dx1_back(x1, x2, h):
  return (f1(x1, x2) - f1(x1-h, x2))/h

def df1_dx2_back(x1, x2, h):
  return (f1(x1, x2) - f1(x1, x2-h))/h

def df2_dx1_back(x1, x2, h):
  return (f2(x1, x2) - f2(x1-h, x2))/h

def df2_dx2_back(x1, x2, h):
  return (f2(x1, x2) - f2(x1, x2-h))/h

def jacobian_back_diff(x1, x2, h):
  return [[df1_dx1_back(x1, x2, h), df1_dx2_back(x1, x2, h)], [df2_dx1_back(x1, x2, h), df2_dx2_back(x1, x2, h)]]

print(jacobian_back_diff(1, 2, 0.001), " backward diff approx")

import symengine

vars = symengine.symbols('x y')
f = symengine.sympify(['cos(x**7) + y**2', '4*sqrt(x)*y**2'])
# f = symengine.sympify(['sin(x)**7 + y**3', '4*x*sqrt(y)'])
J = symengine.zeros(len(f),len(vars))
for i, fi in enumerate(f):
    for j, s in enumerate(vars):
        J[i,j] = symengine.diff(fi, s)
print(J)
print(J.subs({vars[0]: 1, vars[1]: 2}))

def grad_f(x):
  # x^2+9
  return x**2 - 9

def grad_f1(x):
  # f prime
  return 2*x

# do gradient descent with step size
def gradient_descent(x0, f, df, n, step_size):
  x = x0
  for i in range(n):
    print(f"k = {i}, x = {x}, f(x) = {f(x)}, df(x) = {df(x)}")
    x = x - step_size*df(x)
  
gradient_descent(2, grad_f, grad_f1, 10, .25)