mdvsh · April 4, 2023 23:36 · Apr 4, 2023
diff --git a/intro-num-methods.py b/intro-num-methods.py
@@ -0,0 +1,167 @@
+# 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)
No results found