diff --git a/_code/notes/diffsim/derivations.py b/_code/notes/diffsim/derivations.py index 42aeb549cfcab035d7f5de2c7de64bf4966aeb93..499a86beb30fc691c779167da66dfb6c1ffac092 100644 --- a/_code/notes/diffsim/derivations.py +++ b/_code/notes/diffsim/derivations.py @@ -1,20 +1,23 @@ -from sympy import * -from sympy.assumptions.assume import global_assumptions +import sympy as sp +#from sympy.assumptions.assume import global_assumptions def print_expr(name, expr): print(name + ": ") - print(latex(expr)) + print(sp.latex(expr)) print("") -t = Symbol('t', real=True) -dt = Symbol('\Delta t', real=True) -x = Symbol('x', real=True) -xn = Function('x_n', real=True) -f = Function('f', real=True) +t = sp.Symbol('t', real=True) +dt = sp.Symbol('\Delta t', real=True) +x = sp.Symbol('x', real=True) +xn = sp.Function('x_n', real=True) +f = sp.Function('f', real=True) + +x_1, y_1, x_2, y_2 = sp.symbols("x_1 y_1 x_2 y_2", real=True) +G, m_2 = sp.symbols("G m_2", real=True) def euler(): expr = xn(x) + dt * f(t, xn(x)) - return simplify(Derivative(expr, x)) + return sp.simplify(sp.diff(expr, x)) def rk4(): k1 = dt * f(t, xn(x)) @@ -22,7 +25,24 @@ def rk4(): k3 = dt * f(t + dt/2, xn(x) + k2/2) k4 = dt * f(t + dt, xn(x) + k3) expr = xn(x) + k1/6 + k2/3 + k3/3 + k4/6 - return simplify(Derivative(expr, x)) + return sp.simplify(sp.diff(expr, x)) + +def nbody(): + distance = sp.sqrt((x_2 - x_1)**2 + (y_2 - y_1)**2) + f12_x = G * (x_2 - x_1) / distance**3 + f12_y = G * (y_2 - y_1) / distance**3 + df12_x_x1 = sp.simplify(sp.diff(f12_x, x_1)) + df12_x_y1 = sp.simplify(sp.diff(f12_x, y_1)) + df12_y_x1 = sp.simplify(sp.diff(f12_y, x_1)) + df12_y_y1 = sp.simplify(sp.diff(f12_y, y_1)) + print_expr("f12_x", f12_x) + #print_expr("f12_x", f12_y) + print_expr("df12_x_x1", df12_x_x1) + print_expr("df12_x_y1", df12_x_y1) + print_expr("df12_y_x1", df12_y_x1) + print_expr("df12_y_y1", df12_y_y1) + -print_expr("euler", euler()) -print_expr("rk4", rk4()) +#print_expr("euler", euler()) +#print_expr("rk4", rk4()) +nbody()