deep down

in_my_heart.py

import numpy as np
from sympy import *
from sympy.parsing.sympy_parser import parse_expr

def newton(f, x0, eps):
    x = Symbol('x')
    sym_f = parse_expr(f)
    f = lambdify(x, sym_f)
    df = lambdify(x, diff(sym_f, x))
    while 1:
        x1 = x0 - f(x0) / df(x0)
        if abs(x1 - x0) < eps:
            return x1
        x0 = x1

def newton_system(f1, f2, x0, y0, eps):
    x, y = symbols('x y')
    f1_sym, f2_sym = parse_expr(f1), parse_expr(f2)
    f1, f2 = lambdify((x, y), f1_sym), lambdify((x, y), f2_sym)
    df1dx = lambdify((x, y), diff(f1_sym, x))
    df2dx = lambdify((x, y), diff(f2_sym, x))
    df1dy = lambdify((x, y), diff(f1_sym, y))
    df2dy = lambdify((x, y), diff(f2_sym, y))
    while 1:
        g, h = np.linalg.solve(
            np.array([
                [df1dx(x0, y0), df1dy(x0, y0)],
                [df2dx(x0, y0), df2dy(x0, y0)]
            ]),
            np.array([
                -f1(x0, y0),
                -f2(x0, y0)
            ]))
        x1, y1 = x0 + g, y0 + h
        if abs(x1 - x0) + abs(y1 - y0) < eps:
            return x1, y1
        x0, y0 = x1, y1

    
    
def test():
    f = 'tan(3*x) - x**2 + .4'
    # x0 justification: https://www.wolframalpha.com/input/?i=plot+y+%3D+tan(3x)+-+x%5E2+%2B+.4 
    x0 = 0
    eps = 10**(-4)
    print('Root of equation:\nx = %f\n'%(newton(f, x0, eps)))

    f1 = 'sin(y + 2) - x - 1.5'
    f2 = 'y + cos(x - 2) - .5'
    # x0, y0 justif.: https://www.wolframalpha.com/input/?i=plot+%7Bsin(y%2B2)+-+x+%3D1.5,+y+%2B+cos(x-2)+%3D+.5%7D
    x0, y0 = 0, 0
    print('Solution of system:\nx = %f\ny = %f'%(newton_system(f1, f2, x0, y0, eps)))


if __name__ == '__main__':
    test()