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]
%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')
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]
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
r0=[1,0] v0=[0,2] y0=np.concatenate((r0,v0))
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))
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)
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]
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
Choose which booklet to go to: