How to combine vectors and ODEs

14.1 ODEs: Vectors

Sometimes your coupled first order differential equations involve vectors. As an example, consider an object at position $\mathbf{r}$ moving with some time and position dependent velocity, $\mathbf{v}(t,\mathbf{r})$ :

$\frac{d}{dt}\mathbf{r} = \mathbf{v}(t,\mathbf{r})$

It could be a robot at $\mathbf{r}$ that moves with speed

$v = \sqrt{1-\frac{y^2}{y_0^2}}v_0$

( $y$ is the $y$ -coordinate of $\mathbf{r}$ , $y_0$ is a length scale, and $v_0$ is reference speed) in the direction of a moving human, who at time $t$ is at position $\mathbf{r}_h$ :

$\mathbf{r}_h = \left[\begin{array}{c}v_0\cdot t\\0\end{array}\right].$

The human thus moves with constant speed along the $x$ -axis.

The situation thus looks like this:

Figure 1

The vector $\mathbf{d}$ connects the robot with the human:

$\mathbf{d} = \mathbf{r}_h - \mathbf{r},$

and the vector $\mathbf{\hat{d}}$ is a vector of unit length in the direction of $\mathbf{d}$ :

$\mathbf{\hat{d}} = \frac{\mathbf{d}}{\left|\mathbf{d}\right|}.$

Now, wrapped together we have:

$\frac{d}{dt}\mathbf{r} = \sqrt{1-\frac{y^2}{y_0^2}}v_0\>\mathbf{\hat{d}}.$

Solving this equation in Python using solve_ivp turns out to be quite simple if the formulation in terms of vectors is maintained. It may be done like shown here:

import numpy as np
from scipy.integrate import solve_ivp

x0 = 0
y0 = 1
rinit = [x0,y0 - 1e-6] # must be smaller than y0 to get non-zero velocity

tinit = 0
tfinal = 20
trange = [tinit,tfinal]

ts = np.linspace(tinit,tfinal,5)

v0 = 0.1

def drdt(t, r):
    r = np.array(r)
    rp = np.array([v0 * t, 0])
    d = rp - r
    dhat = d/np.linalg.norm(d)
    y = r[1]                        # here y is extracted
    v = np.sqrt(1-y**2/y0**2) * v0
    drdt = v * dhat
    return drdt
    
mysol = solve_ivp(drdt, trange, rinit, max_step=1e-2,
                  t_eval=ts)
ts = mysol.t
xs = mysol.y[0]
ys = mysol.y[1]

Now the solution may be plotted:

%matplotlib inline
import matplotlib.pyplot as plt
plt.rc('font', size=16)

fig, ax = plt.subplots()
ax.set_aspect('equal', adjustable='box')
ax.grid()
ax.plot(xs,ys,linestyle='dashed',marker='o')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xticks(list(np.arange(0,y0*2.2,y0*0.5)))
ax.set_yticks(list(np.arange(0,y0*1.2,y0*0.5)))

ax.text(xs[-1]+0.1,ys[-1],'({:.3f},{:.3f})'.format(xs[-1],ys[-1]),ha='left')

Figure 2

And colored lines connecting the robot and the human may be plotted and their lengths be evaluated and printed:

lens = []
for t,x,y in zip(ts,xs,ys):
    lxs = [x,v0*t]
    lys = [y,0]
    ax.plot(lxs,lys)
    lens.append(np.linalg.norm([np.diff(lxs),np.diff(lys)]))
print(lens)

[0.999999, 0.9993456911625648, 0.9989227148074187, 0.9987202276021728, 0.9986372994529613]

Figure 3

It is seen that with the particular choice of speed, the distance between the robot and the human remains constant.

14.2 Second order ODE with vectors

Consider two objects in space, a sun at the origin of a coordinate system, and a planet at position $\mathbf{r}$ . The sun excerts a gravitational force on the planet according to:

$\mathbf{F} = -\frac{Gm_{sun}m}{\left|\mathbf{r}\right|^3}\mathbf{r}$

where $m_{sun}$ and $m$ are the masses of the sun and the planet, respectively. $G$ is the gravitional force constant.

Assuming that the sun remains at the origin, Newtons 2nd law for this system is:

$m\frac{d^2}{dt^2}\mathbf{r}=\mathbf{F}$

which can be cast in the form of two 1st order differential equations involving the position and velocity vectors, $\mathbf{r}$ and $\mathbf{v}$ :

$\begin{array}{rcl} \displaystyle\frac{d}{dt}\mathbf{r}&=&\mathbf{v}\\ &&\\ \displaystyle\frac{d}{dt}\mathbf{v}&=&-\displaystyle\frac{Gm_{sun}}{\left|\mathbf{r}\right|^3}\mathbf{r} \end{array}$

These first order differential equations may be solved using solve_ivp while keeping strictly to the vector formulation. This will be shown below.

We start by the usual import statements and some constants:

%%capture

import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import norm


import matplotlib.pyplot as plt
from matplotlib import animation, rc
rc('animation', html='jshtml')

G=1
m_sun=10

Then follows the declaration of the initial condition. Since the two first order vectorial differential equations involve a total of four dependent variables, the initial condition is a list of four numbers. Here it is decided that they come in the order $r_x$ , $r_y$ , $v_x$ , and $v_y$ :

r0=[1,0]
v0=[0,2]
y0=np.concatenate((r0,v0))

Having decided on the order in which the dependent variables are stored in the four membered list, the list must be unpacked like that when it arrives as y in the function dydt(t,y). So writing this function, it should start by acknowledging that its second argument, y, is a list whose first two elements are the $\mathrm{r}$ -vector, and whose two last elements are the $\mathrm{v}$ -vector. This is done here and is followed by the mathematical operations on the vectors according to the two differential equations:

def dydt(t,y):
    r = y[:2]
    v = y[2:]
    drdt = v
    dvdt = -G*m_sun/np.power(norm(r),3)*r
    return np.concatenate((drdt,dvdt))

The function ends by assembling with np.concatenate the two vectorial differentials, drdt and dvdt, into a full four-membered list which is returned to the calling solve_ivp.

With the chosen constants, the orbital period of the planet around the sun is $1$ time unit. Here we declare that 150 images are wanted per revolution and that one revolution is desired:

t_end = 1
steps=150*t_end
times=np.linspace(0,t_end,steps)

Finally solve_ivp is called and the resulting $r_x(t)$ and $r_y(t)$ coordinates are extracted as the first two lists in the .y attribute on object returned from solve_ivp:

solution = solve_ivp(dydt, [0, t_end], y0, max_step=1e-3, t_eval=times)
xs=solution.y[0]
ys=solution.y[1]

Had one desired the $v_x(t)$ and $v_y(t)$ solutions, they could have been extracted as solution.y[2] and solution.y[3].

Finally follows an animation:

def update(i):
    line.set_data(xs[0:i+1], ys[0:i+1])
    planet_dot.set_data(xs[i],ys[i])
    return line, planet_dot

fig, ax = plt.subplots(figsize=(4, 4))
line, = ax.plot([], [])
planet_dot, = ax.plot([],[],'ro')
ax.plot(0,0,'bo')
ax.set_xlim([-1.2, 1.2])
ax.set_ylim([-1.2, 1.2])
ax.set_xticks([])
ax.set_yticks([])

anim = animation.FuncAnimation(fig,
                               update,
                               frames=steps,
                               interval=10,
                               blit=True,
                               repeat_delay=0)
anim

and the whole thing eventually looks like this:

Figure 4