Computational Science Asked by Dipole on February 3, 2021
The code below basically illustrates my problem. It is a test code for a pendulum. I solve it using a method suggested on https://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-scipy-integrate-ode). Now I want to plot points at a range of evenly spaced times to show the fact the dynamic nature of the pendulum swings. For example I want to get the solution and plot the corresponding solution state space points at the times
[0,0.2,0.4,0.6,0.8,1.0]
Is there any automatic way Python can do this (using interpolation of the internal solution points which are obtained during the ODE integration)?
Here is the main code:
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
import warnings
def f(t,y):
l = 1
m = 1
d = 1
g = 9.8
return [y[1], -np.sin(y[0])*g/l-y[1]*d/m]
y0, t0 = [np.pi/2, 0], 0
t1 = 500
backend = 'vode'
solver = ode(f).set_integrator(backend, nsteps=1)
solver.set_initial_value(y0, t0)
# suppress Fortran-printed warning
solver._integrator.iwork[2] = -1
y, t = [], []
warnings.filterwarnings("ignore", category=UserWarning)
while solver.t < t1:
solver.integrate(t1, step=True)
y.append(solver.y)
t.append(solver.t)
warnings.resetwarnings()
y = np.array(y)
t = np.array(t)
plt.plot(y[:,0], y[:,1], 'b-')
plt.plot(y[0::5,0], y[0::5,1], 'b.')
plt.show()
The example given in the docs of scipy.integrate.ode hints at using successive integration between your points of interest:
>>> r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
>>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
>>> t1 = 10
>>> dt = 1
>>> while r.successful() and r.t < t1:
>>> r.integrate(r.t+dt)
>>> print("%g %g" % (r.t, r.y))
The scipy.integrate.odeint function takes an input argument 't': A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.
Answered by GertVdE on February 3, 2021
See also: Python: Grid with step control ODE solver
GertVdE's answer covers the user-facing implementation well.
Essentially, you want two different things to happen:
Here's how that works in practice:
You can see all of these details in the SUNDIALS Usage Notes and in the various SUNDIALS Users Guides, which cover the theory behind the integrators as well as their use. Although these sources of documentation cover the successors to VODE (CVODE and CVODES), the theory should be the same, and the interfaces should be similar, except that CVODE and CVODES have additional functionality.
Answered by Geoff Oxberry on February 3, 2021
Another solution is to solve with solve_ivp
and use the dense_output
option which allows interpolating between solution steps:
import numpy as np
from scipy.integrate import ode, solve_ivp
import matplotlib.pyplot as plt
import warnings
def f(t,y):
l = 1
m = 1
d = 1
g = 9.8
return [y[1], -np.sin(y[0])*g/l-y[1]*d/m]
y0, t0 = [np.pi/2, 0], 0
t1 = 100
t = np.linspace(t0,t1,2000)
sol = solve_ivp(f, (t0, t1), y0, dense_output=True)
y_real = sol.y
y_interp = sol.sol(t)
plt.clf()
plt.plot(y_real[0,:],y_real[1,:],'or')
plt.plot(y_interp[0,:],y_interp[1,:],'.-b')
plt.show()
Answered by Arthur Bauville on February 3, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP