Plotting ode solution in phase plane

Question

I am trying to plot the solution to my ode using solver integrate.odeint however I am not obtaining a solution when I attempt to plot it. I can't see where the formulation of my code is wrong.

Please find below: import numpy as np import matplotlib.pyplot as pl from numpy import sin, cos import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate import matplotlib.animation as animation from math import*

g = 9.81 
l = 1.6
l_big = 2.0
l_small = 1.6
m = 0.01
alpha = 0.4
k = 100
def sh2(r1,t):

    theta1,omega1 = r1 
    sh2_theta1 = omega1/(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))**2

    sh2_omega1 = -g*(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))*sin(theta1)

    #return sh2_theta1, sh2_omega1

    return np.array([sh2_theta1, sh2_omega1],float)

init_state = np.radians([69.0,0])
dt = 1/40
time = np.arange(0,10.0,dt)
timexo = np.arange(0,10.0,dt)

state2 = integrate.odeint(sh2,init_state,time)
#print(state2)
print(len(state2),len(timexo))

plt.plot(timexo[0:2500],state2[0:2500])
plt.show()

#phase plot attempt 
# initial values
x_0 = 0 # intial angular position
v_0 = 1.0 # initial angular momentum
t_0 = 0 # initial time

# initial y-vector from initial position and momentum
y0 = np.array([x_0,v_0]) 


# max value of time and points in time to integrate to
t_max = 10
N_spacing_in_t = 1000

# create vector of time points you want to evaluate
t = np.linspace(t_0,t_max,N_spacing_in_t)

# create vector of positions for those times
y_result = np.zeros((len(t), 2))

# loop through all demanded time points
for it, t_ in enumerate(t):

    # get result of ODE integration up to the demanded time
    y = integrate.odeint(sh2,init_state,t_)

    # write result to result vector
    y_result[it,:] = y

# get angle and angular momentum
angle = y_result[:,0]
angular_momentum = y_result[:,1]

# plot result
pl.plot(angle, angular_momentum,'-',lw=1)
pl.xlabel('angle $x$')
pl.ylabel('angular momentum $v$')

pl.gcf().savefig('pendulum_single_run.png',dpi=300)

pl.show()

The output of this code :

Plot 1 : correct plot for the ode solution over time
Plot 2 : empty plot for the phase plane which is whaat is causing problems.

Any hints are welcomed. On a less important note my first plot gives two lines I was only expecting the blue line though.

f.wue f.wue · Accepted Answer · 2019-02-23T10:36:32

The plot is empty beacause the integrator returns zero in your for-loop. Why use a for-loop in the first place? If you integrate over time, like you did in the first part of the code, everything works fine. Note: You don't need to import matplotlib.pyplot twice.

Regarding your two-lines in plot1: the object state2[0:2500] is of dimensions 2x2500. Therefore two lines appear. If you only want on of the lines, use np.transpose and then state2[0] or state2[1] to access the lines.

The code below will solve your problem. I added a plt.figure() command to generate both plots without waiting for the first clot to be closed.

import matplotlib.pyplot as plt
import numpy as np
from numpy import sin
import scipy.integrate as integrate
from math import *

g = 9.81
l = 1.6
l_big = 2.0
l_small = 1.6
m = 0.01
alpha = 0.4
k = 100
def sh2(r1,t):

    theta1,omega1 = r1 
    sh2_theta1 = omega1/(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))**2

    sh2_omega1 = -g*(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))*sin(theta1)

    #return sh2_theta1, sh2_omega1

    return np.array([sh2_theta1, sh2_omega1],float)

init_state = np.radians([69.0,0])
dt = 1/40
time = np.arange(0,10.0,dt)
timexo = np.arange(0,10.0,dt)

state2 = integrate.odeint(sh2,init_state,time)

print(len(state2),len(timexo))
state2_plot = np.transpose(state2[0:2500])
plt.plot(timexo[0:2500],state2_plot[1])

#phase plot attempt 
# initial values
x_0 = 0.0 # intial angular position
v_0 = 1.0 # initial angular momentum
t_0 = 0 # initial time

# initial y-vector from initial position and momentum
y0 = np.array([x_0,v_0]) 

# max value of time and points in time to integrate to
t_max = 10
N_spacing_in_t = 1000

# create vector of time points you want to evaluate
t = np.linspace(t_0,t_max,N_spacing_in_t)

# create vector of positions for those times

y_result = integrate.odeint(sh2, init_state, t)

# get angle and angular momentum
angle = y_result[:,0]
angular_momentum = y_result[:,1]

# plot result
fig = plt.figure()
plt.plot(angle, angular_momentum,'-',lw=1)
plt.xlabel('angle $x$')
plt.ylabel('angular momentum $v$')

plt.gcf().savefig('pendulum_single_run.png',dpi=300)
plt.show()

Output:

Plotting ode solution in phase plane

1 Answers