# Animating the trajectory of solutions to ordinary differential equations

Hello Blender Community,

I am new to the software and I am looking for some advice on the following problem.

Suppose that we are given the position and orientation of a rigid object in Euclidean space $$(x,y,z)$$ which are represented by the coordinates $$(\boldsymbol{r},\Lambda)$$ respectively. Further, let us assume that the evolution of the position and orientation is governed by an ordinary differential equation of the form $$(\dot{\boldsymbol{r}},\dot{\Lambda}) = f(\boldsymbol{r},\Lambda)$$ where dot's refer to time derivatives $$\frac{d}{dt}$$.

We then can solve this ODE numerically using, say, the Runge-Kutta algorithm. This then gives us a time indexed sequence of vectors and matrices, namely,
$$\boldsymbol{r}_i = \begin{bmatrix}x_i\\y_i \\ z_i \end{bmatrix}$$ for the object's position and

$$\Lambda_i = \begin{bmatrix}\Lambda_{11}^i& \Lambda_{12}^i & \Lambda_{13}^i\\\Lambda_{21}^i& \Lambda_{22}^i & \Lambda_{23}^i\\ \Lambda_{31}^i& \Lambda_{32}^i & \Lambda_{33}^i \end{bmatrix}$$ for it's orientation. If one were to take enough points $$i$$, then this sequence approximates the trajectory of the object (its position and orientation in space) by applying these sets of transformations to a chosen $$(x,y,z)$$ coordinate system.

I would like to animate the trajectory but I need some help/examples as to how to implement. I think a possible solution is as follows:

Step 1: Either solve the ODE in Python and import an array of solution points as a csv file OR set up the problem using C# (I think this is what Blender uses though I am at loss as to how to interface with it), implement Runga-Kutta and get an time indexed array of state-points.

Step 2: Interpolate the missing points to obtain a smooth trajectory.

Step 3: Fix a coordinate system to the animated object

Step 4: Somehow apply coordinate transformations $$(\boldsymbol{r},\Lambda)$$ to the object

Step 5: run the simulation from which we obtain an animation the trajectory of a rigid body in space.

I know how to solve ODE's in Python or MATLAB (Step 1) but I am at a loss as to how to complete steps 2-5. If anyone has any expertise with such a problem or can suggest any reading material/resources that would be greatly appreciated. Thank you.

Let me give an example. Let us consider the motion of a free rigid body. Euler tells us that the evolution of the body's motion is governed by $$\dot{\Omega} = \mathbb{I}\Omega\times\Omega$$ $$\dot{\Lambda} = \widehat{\Omega}\Lambda$$ where $$\widehat{\Omega} = \Omega\times = \epsilon_{ijk}\Omega_k$$, the antisymmetic tensor. To solve this equation we use the following routine in python:

from scipy.integrate import odeint
import numpy as np
# define the hat map
# takes in 3 vector and returns anti symmetric matrix
def hat(x):
hat = np.array([[0,-x[2],x[1],[x[2],0,-x[0]],[-x[1],x[0],0]])
return hat

# define vector field
# x[0] - x[2] Omega1-3
# x[3] - x[11] Lambda
# I moment of inertia

def f(x,t,I):
dfdt = np.zeros(12)
Om1 = x[0]
Om2 = x[1]
Om3 = x[2]
La11 = x[3]
La12 = x[4]
La13 = x[5]
La21 = x[6]
La22 = x[7]
La23 = x[8]
La31 = x[9]
La32 = x[10]
La33 = x[11]
Om = np.array([Om1,Om2,Om3])
La = np.array([[La11,La12,La13],[La21,La22,La23],[La31,La32,La33]])
Pi = np.matmult(I,Om)
X = np.cross(Pi,Om)
# Dynamics
dfdt[0:4] = X
# Reconstruction equations
hatOm = hat(Om)
Y = np.matmul(hatOm,La)
Ry = np.reshape(Y,(9))
dfdt[4:12] = Ry
return dfdt

# simulate rigid body
# time interval
tspan = np.linspace(0,1,100)
# initial conditions
IC = np.zeros(12)
IC[0] = 1
# moment of inertia
I = np.array([[1,0,0],[0,1,0],[0,0,1]])
x = odeint(f,IC,tspan,args=(I))
Omega1 = x[0:,0]
Omega2 = x[0:,1]
Omega3 = x[0:,2]
La11 = x[0:,3]
La12 = x[0:,4]
La13 = x[0:,5]
La21 = x[0:,6]
La22 = x[0:,7]
La23 = x[0:,8]
La31 = x[0:,9]
La32 = x[0:,10]
La33 = x[0:,11]


My question is how to attach these simulated values to an object to be animated. I have tried the following code, yet the animation it produces does not seem realistic. I am not sure where the issue is.

# lets look at the motion in velocity space
import bpy
ob = bpy.data.objects["Sphere"]
start_position = (1,0,0)
frame_num = 0
t0 = len(tspan)
for i in range(t0-1):
bpy.context.scene.frame_set(frame_num)
ob.location = (Omega1[i],Omega2[i],Omega3[i])
ob.keyframe_insert(data_path='location',index = -1)
frame_num +=20

• Good question, the first part of your code has errors. Could you fix them please? Or even better, could your replace it with a simpler equation? – Leander Oct 28 '19 at 16:07
• What I don't understand is, that you have timesteps of 0.01 but jump 20 frames when settings keyframes. – Leander Oct 28 '19 at 16:08
• As for the motion path, simple create a new object and populate it with your samples, then connect the vertices. Bmesh. – Leander Oct 28 '19 at 16:12