Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions CHANGELOG
Original file line number Diff line number Diff line change
@@ -1,4 +1,7 @@
--- CHANGELOG ---
--- Assimulo-3.4.3 ---
* Improved compliance with newer scipy version by instead using corresponding numpy calls when more suitable.

--- Assimulo-3.4.2 ---
* Updated an import statement in one of the examples to resolve compliance issues with scipy 1.10.1.

Expand Down
81 changes: 33 additions & 48 deletions examples/cvode_gyro.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,7 @@
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.

from scipy import *
from numpy import hstack as numpy_hstack
import numpy as np
import nose
from assimulo.problem import Explicit_Problem
from assimulo.solvers import CVode
Expand All @@ -35,74 +34,60 @@ def run_example(with_plots=True):
"""

def curl(v):
return array([[0,v[2],-v[1]],[-v[2],0,v[0]],[v[1],-v[0],0]])
return np.array([[0, v[2], -v[1]], [-v[2], 0, v[0]], [v[1], -v[0], 0]])

#Defines the rhs
def f(t,u):
# Defines the rhs
def f(t, u):
"""
Simulations for the Gyro (Heavy Top) example in Celledoni/Safstrom:
Journal of Physics A, Vol 39, 5463-5478, 2006
"""
I1=1000.
I2=5000.
I3=6000.
u0=[0,0,1.]
pi=u[0:3]
Q=(u[3:12]).reshape((3,3))
Qu0=dot(Q,u0)
f=array([Qu0[1],-Qu0[0],0.])
f=0
omega=array([pi[0]/I1,pi[1]/I2,pi[2]/I3])
pid=dot(curl(omega),pi)+f
Qd=dot(curl(omega),Q)
return numpy_hstack([pid,Qd.reshape((9,))])
I1 = 1000.
I2 = 5000.
I3 = 6000.
u0 = [0, 0, 1.]
pi = u[0:3]
Q = (u[3:12]).reshape((3, 3))
Qu0 = np.dot(Q, u0)
omega = np.array([pi[0]/I1, pi[1]/I2, pi[2]/I3])
pid = np.dot(curl(omega), pi)
Qd = np.dot(curl(omega),Q)
return np.hstack([pid, Qd.reshape((9,))])

def energi(state):
energi=[]
for st in state:
Q=(st[3:12]).reshape((3,3))
pi=st[0:3]
u0=[0,0,1.]
Qu0=dot(Q,u0)
V=Qu0[2] # potential energy
T=0.5*(pi[0]**2/1000.+pi[1]**2/5000.+pi[2]**2/6000.)
energi.append([T])
return energi

#Initial conditions
y0=numpy_hstack([[1000.*10,5000.*10,6000*10],eye(3).reshape((9,))])
# Initial conditions
y0 = np.hstack([[1000.*10, 5000.*10, 6000*10], np.eye(3).reshape((9,))])

#Create an Assimulo explicit problem
exp_mod = Explicit_Problem(f,y0, name="Gyroscope Example")
# Create an Assimulo explicit problem
exp_mod = Explicit_Problem(f, y0, name="Gyroscope Example")

#Create an Assimulo explicit solver (CVode)
exp_sim=CVode(exp_mod)
# Create an Assimulo explicit solver (CVode)
exp_sim = CVode(exp_mod)

#Sets the parameters
exp_sim.discr='BDF'
exp_sim.iter='Newton'
exp_sim.maxord=2 #Sets the maxorder
exp_sim.atol=1.e-10
exp_sim.rtol=1.e-10
# Sets the parameters
exp_sim.discr = 'BDF'
exp_sim.iter = 'Newton'
exp_sim.maxord = 2 # Sets the maxorder
exp_sim.atol = 1.e-10
exp_sim.rtol = 1.e-10

#Simulate
# Simulate
t, y = exp_sim.simulate(0.1)

#Plot
# Plot
if with_plots:
import pylab as P
P.plot(t,y/10000.)
P.plot(t, y/10000.)
P.xlabel('Time')
P.ylabel('States, scaled by $10^4$')
P.title(exp_mod.name)
P.show()

#Basic tests
nose.tools.assert_almost_equal(y[-1][0],692.800241862)
nose.tools.assert_almost_equal(y[-1][8],7.08468221e-1)
nose.tools.assert_almost_equal(y[-1][0], 692.800241862)
nose.tools.assert_almost_equal(y[-1][8], 7.08468221e-1)

return exp_mod, exp_sim


if __name__=='__main__':
mod,sim = run_example()
mod, sim = run_example()