Commit a029ac50 authored by Erik Strand's avatar Erik Strand
Browse files

Add derivations of n-body jacobians

Specifically, partial derivatives of the state gradient (so velocities
and accelerations) w.r.t. the state variables (so displacements and
velocities).
parent f8eb49e4
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()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment