Solving coupled ODES which varies with time in Python

Question

While trying to solve 6 couple equation using odeint it say "p = odeint(odes,r0,θ0,ϕ0,x0,z0,t) NameError: name 'r0' is not defined"

Code: from scipy.integrate import odeint import numpy as np import matplotlib.pyplot as plt from math import sin, cos, pi

def odes(p, t):
# constants
Ω=9.74e-3
B_θ=-8.6e-6*sin(θ)
B_r=25893.2e-9*cos(θ)
β=-9.36e-10
# assigning each ODE to a vector element
r = p[1]
θ = p[3]
ϕ = p[5]
x = p[2]
y = p[4]
z = p[6]

# defining the ODEs
drdt = x
dxdt = r*(y**2 + (z+Ω)**2 * sin(θ)**2 - β*z*sin(θ)*B_θ)
dθdt = y
dydt = (-2*x*y-r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r
dϕdt = z
dzdt = (-2*x*(z+Ω)*sin(θ)-2*r*y*(z+Ω)*cos(θ)+β*(x*B_θ-r*y*B_r))/(r*sin(θ))

return [drdt, dxdt, dθdt, dydt, dϕdt, dzdt]

# initial conditions
r0 = 0.71e+8
θ0 = 0.5*pi
ϕ0 = 0
x0 = 0
y0 = 0
z0 = 0

# time window
t = np.linspace(0,15,1000)
p = odeint(odes,r0,θ0,ϕ0,x0,z0,t)

r=p[:,1]
θ=p[:,3]
ϕ=p[:,5]
x=p[:,2]
y=p[:,4]
z=p[:,6]

# plot the results
plt.semilogy(t,r)
plt.semilogy(t,x)
plt.semilogy(t,θ)
plt.semilogy(t,y)
plt.semilogy(t,ϕ)
plt.semilogy(t,z)

plt.show()

Answer 1

There are multiple errors in your code that prevent it from running:

• The indentation is wrong. This could be a copy-paste error, but in python indentation is a syntax element.

• θ is used in "constants" before it is defined as "vector element" (or state component)

• python is zero-based, the array indices go from 0 to 5, not 1 to 6

• odeint only takes one initial state variable. In the vector case this has to be a tuple or array

Next the components have wildly different scales. This is of no import for a fixed-step solver, but can lead to distortions in a step-size controller, as some components get shadowed by the largest or fastest changing components. Scale the components individually, so that the solver sees only uniformly scaled components. With these changes I get to a plot, but that does not look realistic, the radius increases rapidly while the simulations are for a stable radius with small oscillations.

scales = np.array([1e7, 0.1, 1, 1e-20, 10, 0.1])

def odes(p, t):
    # assigning each ODE to a vector element
    r,x,θ,y,ϕ,z = p*scales
    # constants
    Ω=9.74e-3
    B_θ=-8.6e-6*sin(θ)
    B_r=25893.2e-9*cos(θ)
    β=-9.36e-10

    # defining the ODEs
    drdt = x
    dxdt = r*(y**2 + (z+Ω)**2 * sin(θ)**2 - β*z*sin(θ)*B_θ)
    dθdt = y
    dydt = (-2*x*y-r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r
    dϕdt = z
    dzdt = (-2*x*(z+Ω)*sin(θ)-2*r*y*(z+Ω)*cos(θ)+β*(x*B_θ-r*y*B_r))/(r*sin(θ))

    return np.array([drdt, dxdt, dθdt, dydt, dϕdt, dzdt])/scales

# initial conditions
r0 = 7.1e+07 
θ0 = 0.5*pi
ϕ0 = 0
x0 = 0
y0 = 0
z0 = 0

# time window
t = np.arange(0,100*60,10)
p0 = np.array([r0,x0,θ0,y0,ϕ0,z0])/scales
p = odeint(odes,p0,t, atol=1e-8, rtol=1e-10)

r,x,θ,y,ϕ,z = p.T*scales[:,None]

# plot the results
fig,ax=plt.subplots(2,3,figsize=(6,4))
for a,u in zip(ax.flatten(),[r,x,θ,y,ϕ,z]):
    a.plot(t,u); a.grid()

plt.tight_layout(); plt.show()