Link to booklet is copied to clipboard!

Saved to booklet!

Removed from booklet!

Started capturing to booklet.

Stopped capturing to booklet.

Preferring latest booklet.

No longer prefers latest booklet.

Entered edit mode.

Exited edit mode.

Sometimes your coupled first order differential equations involve vectors.
As an example, consider an object at position
moving with some time and position dependent velocity,
:

It could be a robot at that moves with speed

( is the -coordinate of , is a
length scale, and is reference speed)
in the direction of a moving human, who at time is at
position :

The human thus moves with constant speed along the -axis.
The situation thus looks like this:

The vector connects the robot with the human:

and the vector is a vector
of unit length in the direction of :

Now, wrapped together we have:

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')

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]

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

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

where and are the masses of the sun and the planet,
respectively. is the gravitional force constant.
Assuming that the sun remains at the origin, Newtons 2nd law for this
system is:

which can be cast in the form of two 1st order differential equations
involving the position and velocity vectors, and :

These first order differential equations may be solved using
We start by the usual import statements and some constants:

`solve_ivp`

while keeping strictly to the vector
formulation. This will be shown below.
%%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 , ,
, and :

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 -vector, and whose two last
elements are the -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
With the chosen constants, the orbital period of the planet around the
sun is time unit. Here we declare that 150 images are wanted per
revolution and that one revolution is desired:

`np.concatenate`

the two vectorial differentials, `drdt`

and
`dvdt`

, into a full four-membered list which is returned
to the calling `solve_ivp`

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

Finally

`solve_ivp`

is called and the resulting and
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 and solutions, they could have
been extracted as
Finally follows an animation:

`solution.y[2]`

and
`solution.y[3]`

.
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:

Once deleted this booklet is gone forever!

Choose which booklet to go to: