python22.sci2u.dk Active sci2u
Loading...
Login
Prefer latest booklet Prefer latest booklet Login
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.
Edit mode

Click on the booklet title or a chapter title to edit (max. 30 characters)


How To's

How to get python How to plot points and lines How to use loops How to draw general shapes How to draw curves How to use functions How to manipulate Polygons How to create animations How to do mathematical analysis How to fit How to go from SymPy to NumPy How to solve a single differential equation How to do coupled and 2nd order ODEs How to combine vectors and ODEs
14.1 ODEs: Vectors
14.2 Second order ODE with vectors
How to do ODE events Add Chapter..

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

Figure 1

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

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



Sci2u Assignment: 820
Delete "How To's"?
Once deleted this booklet is gone forever!
Block is found in multiple booklets!

Choose which booklet to go to:

© 2019-2022 Uniblender ApS