Skip to content
Browse files

Consistent naming of options, convert slidingpendulum to new API

  • Loading branch information...
1 parent 28b281a commit e826a411b6bb3eeb352f8c4666e55e068f1fad71 @bmcage committed Jan 24, 2012
View
4 docs/src/examples/ddaspk/planarpendulum_ddaspk.py
@@ -75,8 +75,8 @@ def res(self, t, x, xdot, tmp):
ig = dae('ddaspk', problem.res)
#first compute the correct initial condition from the values of z0
-ig.set_options(algebraic_var=array([1,1,1,1,algvar]),
- compute_initcond='yode0',
+ig.set_options(algebraic_vars_idx=[4],
+ compute_initcond='y0',
first_step=1e-9,
atol=1e-6, rtol=1e-6)
#initialize solver and obtain correct values of initial condition
View
8 docs/src/examples/doublependulum.py
@@ -334,6 +334,7 @@ def main():
first_step=1e-18,
atol=ATOL*jfac,
rtol=RTOL*jfac,
+ max_steps=1500,
jacfn=jac,
algebraic_vars_idx=problem.algvar_idx,
exclude_algvar_from_error=problem.exclalg_err,
@@ -392,11 +393,12 @@ def main():
problem.set_jac(jac)
ig.set_options(jacfn=problem.ddaspk_jac)
#first compute the correct initial condition from the values of z0
- ig.set_options(algebraic_var=problem.algvar,
- compute_initcond='yode0',
+ ig.set_options(
+ algebraic_vars_idx=problem.algvar_idx,
+ compute_initcond='yp0',
first_step=1e-18,
exclude_algvar_from_error=problem.exclalg_err,
- atol=ATOL, rtol=RTOL, nsteps=1500)
+ atol=ATOL, rtol=RTOL, max_steps=1500)
tinit = ig.init_step(0., problem.z0, problem.zprime0, ddaspkz[0], ddaspkzprime[0])
print('ddaspk started from z0 = ', problem.z0)
View
152 docs/src/examples/slidingpendulum.py
@@ -39,7 +39,7 @@
except:
pass
#imports
-from numpy import (arange, zeros, array, sin, cos, asarray, sqrt, pi)
+from numpy import (arange, zeros, array, sin, cos, asarray, sqrt, pi, empty)
from scikits.odes import dae
import pylab
import os
@@ -50,7 +50,7 @@ class Slidingpendulum():
#default values
deftend = 120.
- deftstep = 1e-2
+ deftstep = 1e-3
defg = 9.8
defm = 1.0
defM = 1.0
@@ -95,7 +95,7 @@ def __init__(self, data=None, type='index2'):
self.g = data.defg
self.omega = data.defomega
- self.stop_t = arange(.0, self.tend, self.tstep)[1:]
+ self.stop_t = arange(.0, self.tend, self.tstep)
lambdaval = 0.0
@@ -105,25 +105,24 @@ def __init__(self, data=None, type='index2'):
0., 0., lambdaval, lambdaval])
self.zprime0 = array([0., 0., 0., 0., -self.g, 0., 0., 0.], float)
self.res = self.resindex2
- self.algvar = array([1,1,1,1,1,1,-1,-1])
+ self.algvaridx = [6,7]
self.exclalg_err = True
else:
self.neq = 7
self.z0 = array([self.x0, self.y0, self.theta0, 0.,
0., 0., lambdaval])
self.zprime0 = array([0., 0., 0., 0., -self.g, 0., 0.], float)
self.res = self.resindex1
- self.algvar = array([1,1,1,1,1,1,-1])
- self.exclalg_err = False
+ self.algvaridx = [6]
+ self.exclalg_err = True
- def resindex2(self, tres, yy, yp):
+ def resindex2(self, tres, yy, yp, res):
m = self.m
M = self.M
l = self.l
omega = self.omega
g = self.g
k = self.k
- tmp = zeros(self.neq)
## x = yy[0]
## y = yy[1]
## theta = yy[2]
@@ -133,36 +132,36 @@ def resindex2(self, tres, yy, yp):
## it is needed to scale the constraint by blowing up
## the lagrange multiplier with 1e10
- tmp[0]= ((M+m)*yp[3] + m*l*cos(yy[2])*yp[5] \
+ res[0]= ((M+m)*yp[3] + m*l*cos(yy[2])*yp[5] \
- m*l*sin(yy[2])*yy[5]**2 \
+ yy[6]*1e10*(-2*yy[0]+omega/3.*sin(omega*yy[0]))\
+k*yy[3])
- tmp[1]= ((M+m)*yp[4] + m*l*sin(yy[2])*yp[5] \
+ res[1]= ((M+m)*yp[4] + m*l*sin(yy[2])*yp[5] \
+ m*l*cos(yy[2])*yy[5]**2 + (M+m)*g \
+ yy[6]*1e10\
+k*yy[4] )
- tmp[2]= m*l*cos(yy[2])*yp[3] + m*l*sin(yy[2])*yp[4] + m*l**2 * yp[5] \
+ res[2]= m*l*cos(yy[2])*yp[3] + m*l*sin(yy[2])*yp[4] + m*l**2 * yp[5] \
+ m*g*sin(yy[2])
#stabalizing factor
- tmp[3]= yp[0] - yy[3] + yy[7]*yy[0]
- tmp[4]= yp[1] - yy[4] + yy[7]*yy[1]
- tmp[5]= yp[2] - yy[5]
+ res[3]= yp[0] - yy[3] + yy[7]*yy[0]
+ res[4]= yp[1] - yy[4] + yy[7]*yy[1]
+ res[5]= yp[2] - yy[5]
#use both constraint and it's derivative
- tmp[6]= (yy[1]-yy[0]**2-1/3.*cos(omega*yy[0]))/10000.
- tmp[7]= (-2*yy[0] + omega/3.*sin(omega*yy[0]))*yy[3] + yy[4]
+ res[6]= (yy[1]-yy[0]**2-1/3.*cos(omega*yy[0]))/10000.
+ res[7]= (-2*yy[0] + omega/3.*sin(omega*yy[0]))*yy[3] + yy[4]
##print tres, 'yy', yy
##print tres, 'yp', yp
- ##print tres, 'res', tmp
- return tmp
+ ##print tres, 'res', res
+ return 0
- def resindex1(self, tres, yy, yp):
+ def resindex1(self, tres, yy, yp, res):
m = self.m
M = self.M
l = self.l
omega = self.omega
g = self.g
k = self.k
- tmp = zeros(self.neq)
+ res = zeros(self.neq)
## x = yy[0]
## y = yy[1]
## theta = yy[2]
@@ -172,33 +171,33 @@ def resindex1(self, tres, yy, yp):
## it is needed to scale the constraint by blowing up
## the lagrange multiplier with 1e10
- tmp[0]= (M+m)*yp[3] + m*l*cos(yy[2])*yp[5] \
+ res[0]= (M+m)*yp[3] + m*l*cos(yy[2])*yp[5] \
- m*l*sin(yy[2])*yy[5]**2 \
+ yy[6]*1e10*(-2*yy[0]+omega/3.*sin(omega*yy[0]))\
+k*yy[3]
- tmp[1]= (M+m)*yp[4] + m*l*sin(yy[2])*yp[5] \
+ res[1]= (M+m)*yp[4] + m*l*sin(yy[2])*yp[5] \
+ m*l*cos(yy[2])*yy[5]**2 + (M+m)*g \
+ yy[6]*1e10\
+k*yy[4]
- tmp[2]= m*l*cos(yy[2])*yp[3] + m*l*sin(yy[2])*yp[4] + m*l**2 * yp[5] \
+ res[2]= m*l*cos(yy[2])*yp[3] + m*l*sin(yy[2])*yp[4] + m*l**2 * yp[5] \
+ m*g*sin(yy[2])
- tmp[3]= yp[0] - yy[3]
- tmp[4]= yp[1] - yy[4]
- tmp[5]= yp[2] - yy[5]
+ res[3]= yp[0] - yy[3]
+ res[4]= yp[1] - yy[4]
+ res[5]= yp[2] - yy[5]
#following equation is wrong, we need to derive another time and
#eliminate the dot{xdot} and dot{ydot} so that lagrange multiplier appears
- tmp[6]= (-2*yy[0] + omega/3.*sin(omega*yy[0]))*yy[3] + yy[4]
+ res[6]= (-2*yy[0] + omega/3.*sin(omega*yy[0]))*yy[3] + yy[4]
##sys.exit()
##print tres, 'yy', yy
##print tres, 'yp', yp
- ##print tres, 'res', tmp
- return tmp
+ ##print tres, 'res', res
+ return 0
def main():
"""
The main program: instantiate a problem, then use odes package to solve it
"""
- uinput = input("Solve as\n 1 = index 2 problem\n 2 = index 1 problem\n"
+ uinput = input("Solve as\n 1 = index 2 problem \n 2 = index 1 problem (Not working!)\n"
" \n 4 = info\n\n"
"Answer (1,2 or 4) : ")
if uinput == '1':
@@ -209,80 +208,45 @@ def main():
print(__doc__)
return
- input1 = input("Solve with\n 1 = ddaspk\n 2 = ida\n\n"
+ input1 = input("Solve with\n 1 = ida\n 2 = ddaspk\n\n"
"Answer (1 or 2) : ").strip()
if input1 not in ["1", "2"]:
print("Invalid solution method given")
return
-
- z = [0]*(1+len(problem.stop_t)); zprime = [0]*(1+len(problem.stop_t))
- ig = dae(problem.res, problem.jac)
- #first compute the correct initial condition from the values of z0
- if input1 == "2":
- ig.set_integrator('odesIDA',algebraic_var=problem.algvar,
- compute_initcond='yode0',
- first_step=1e-9,
- atol=1e-5,rtol=1e-5)
- if input1 == "1":
- #first compute the correct initial condition from the values of z0
- ig.set_integrator('ddaspk',algebraic_var=problem.algvar,
- compute_initcond='yode0',
- first_step=1e-9,
- atol=1e-2,rtol=1e-2)
- ig.set_initial_value(problem.z0, problem.zprime0, t=0.0)
+ if input1 == '1':
+ ig = dae('ida', problem.res, atol=1e-5,rtol=1e-4)
+ elif input1 == '2':
+ ig = dae('ddaspk', problem.res, atol=1e-5,rtol=1e-4)
+ ig.set_options(jacfn=problem.jac,
+ algebraic_vars_idx=problem.algvaridx,
+ compute_initcond='yp0',
+ exclude_algvar_from_error=problem.exclalg_err,
+ max_steps = 15000,
+ first_step=1e-9)
- i=0
- z[i], zprime[i] = ig.solve(1e-9);
- assert ig.successful(), (problem,)
+ #Solve it
+ result= ig.solve(problem.stop_t, problem.z0, problem.zprime0)
+ time = result[1]
+ z = result[2]
+ zprime = result[3]
+
+ #some user output
print('started from z0 = ', problem.z0)
print('initial condition calculated, [z,zprime] = [', z[0], zprime[0], ']')
- if input1 == "2":
- ig.set_integrator('odesIDA',algebraic_var=problem.algvar,
- first_step=1e-9,
- atol=1e-5,
- rtol=1e-5,
- exclude_algvar_from_error=problem.exclalg_err,
- nsteps = 1500)
- if input1 == "1":
- ig.set_integrator('ddaspk',algebraic_var=problem.algvar,
- first_step=1e-9,
- atol=1e-2,
- rtol=1e-2,
- exclude_algvar_from_error=problem.exclalg_err,
- nsteps = 1500)
- ig.set_initial_value(z[0], zprime[0], t=0.0)
-
- i=1
+ print('last sol at time', time[-1])
+ print(' has solution ', z[-1], zprime[-1])
+ res = empty(problem.neq, float)
+ problem.res(problem.stop_t[-1], z[-1], zprime[-1], res)
+ print(' has residual: ', res)
- error = False
- for time in problem.stop_t:
- #print 'at time', time
- z[i], zprime[i] = ig.solve(time)
- #print 'sol at ', time, z[i]
- i += 1
- if not ig.successful():
- error = True
- print('Error in solver, breaking solution at time %g' % time)
- break
+ xt = z[:,0]
+ yt = z[:,1]
+ thetat = z[:,2]
+ nr = len(xt)
-
- print('last sol', z[i-1], zprime[i-1])
- print('has residual: ', problem.res(problem.stop_t[i-2], z[i-1],
- zprime[i-1]))
-
- nr = i
- xt = [z[i][0] for i in range(nr)]
- yt = [z[i][1] for i in range(nr)]
- thetat = [z[i][2] for i in range(nr)]
- time = zeros(nr,float)
- time[0] = 0.0
- if error:
- time[1:] = problem.stop_t[:nr-1]
- else:
- time[1:] = problem.stop_t[:nr]
-
+ pylab.ion()
pylab.figure(1)
pylab.subplot(111)
pylab.scatter(xt, yt)
@@ -367,7 +331,7 @@ def create_animation(sizex, sizey, ext):
secs = 0
frame = 0
print('Generating output ...\n')
- for solnr in range(0,nr,5):
+ for solnr in range(0,nr,50):
drawonesol(solnr, sizex, sizey, frame)
frame += 1
if solnr // 500 != secs :
View
81 scikits/odes/ddaspkint.py
@@ -69,13 +69,13 @@
in packed format, jac_packed[i-j+lband, j] = jac[i,j].
- tcrit : None or float. If given, tcrit is a critical time point
beyond which no integration occurs
-- nsteps : int
+- max_steps : int
Maximum number of (internally defined) steps allowed during one
call to the solver.
- first_step : float
Set initial stepsize. DAE solver can suffer on the first step, set this
to circumvent this.
-- max_step : float
+- max_step_size : float
Limits for the step sizes used by the integrator. This overrides the
internal value
- order : int
@@ -84,11 +84,11 @@
- enforce_nonnegativity: bool
Enforce the nonnegativity of Y during integration
Note: best is to run code first with no condition
-- compute_initcond: None or 'yprime0' or 'yode0'
+- compute_initcond: None or 'yp0' or 'y0'
DDASPK may be able to compute the initial conditions if you do not know them
precisely.
- If yprime0, then y0 will be calculated
- If yode0, then the differential variables will be used to solve for the
+ If y0, then y0 will be calculated
+ If yp0, then the differential variables will be used to solve for the
algebraic variables and the derivative of the differential variables.
Which are the algebraic variables must be indicated with algebraic_var method
- exclude_algvar_from_error: bool
@@ -105,9 +105,14 @@
-2: y0[i] < 0
0: y0[i] not constrained
Alternatively, pass only one integer that then applies to all unknowns
-- algebraic_var: integer array of length the number of unknowns, indicating the
- algebraic variables in y. Use -1 to indicate an algebraic variable, and +1 for
- a differential variable.
+- algebraic_vars_idx: an array or None (= default)
+ Description:
+ If given problem is of type DAE, some items of the residual
+ vector returned by the 'resfn' have to be treated as
+ algebraic variables. These are denoted by the position
+ (index) in the residual vector.
+ All these indexes have to be specified in the
+ 'algebraic_vars_idx' array.
"""
__doc__ += integrator_info
@@ -179,24 +184,24 @@ class ddaspk(DaeBase):
supports_run_relax = 0
supports_step = 1
- def __init__(self,resfn, **options):
+ def __init__(self, resfn, **options):
default_values = {
- 'rtol':1e-6,
- 'atol':1e-12,
- 'lband':None,
- 'uband':None,
- 'tcrit':None,
+ 'rtol': 1e-6,
+ 'atol': 1e-12,
+ 'lband': None,
+ 'uband': None,
+ 'tcrit': None,
'order' : 5,
- 'nsteps' : 500,
- 'max_step' : 0.0, # corresponds to infinite
+ 'max_steps' : 500,
+ 'max_step_size' : 0.0, # corresponds to infinite
'first_step' : 0.0, # determined by solver
- 'enforce_nonnegativity':False,
- 'nonneg_type':None,
- 'compute_initcond':None,
- 'constraint_init':False,
- 'constraint_type':None,
- 'algebraic_var':None,
- 'exclude_algvar_from_error':False,
+ 'enforce_nonnegativity': False,
+ 'nonneg_type': None,
+ 'compute_initcond': None,
+ 'constraint_init': False,
+ 'constraint_type': None,
+ 'algebraic_vars_idx': None,
+ 'exclude_algvar_from_error': False,
'rfn': None,
'jacfn': None,
}
@@ -228,8 +233,8 @@ def _init_data(self):
if self.options['order'] > 5 or self.options['order'] < 1:
raise ValueError('order should be >=1, <=5')
self.order = self.options['order']
- self.nsteps = self.options['nsteps']
- self.max_step = self.options['max_step']
+ self.nsteps = self.options['max_steps']
+ self.max_step = self.options['max_step_size']
self.first_step = self.options['first_step']
self.nonneg =0
if self.options['enforce_nonnegativity'] and self.options['constraint_init']:
@@ -244,18 +249,17 @@ def _init_data(self):
else: self.constraint_type = self.options['constraint_type']
if self.options['compute_initcond'] is None:
self.compute_initcond = 0
- elif re.match(self.options['compute_initcond'], r'yprime0', re.I):
+ elif re.match(self.options['compute_initcond'], r'y0', re.I):
self.compute_initcond = 2
- elif re.match(self.options['compute_initcond'], r'yode0', re.I):
+ elif re.match(self.options['compute_initcond'], r'yp0', re.I):
self.compute_initcond = 1
else:
raise ValueError('Unknown init cond calculation method %s' %(
self.options['compute_initcond']))
- if self.compute_initcond == 1 and self.options['algebraic_var'] is None:
- raise ValueError('Give integer array indicating which are the '\
- 'algebraic variables, +1 for diffential var, '\
- '-1 for algebraic var')
- self.algebraic_var = self.options['algebraic_var']
+ if self.compute_initcond == 1 and not self.options['algebraic_vars_idx']:
+ raise ValueError('Give array indicating which are the '\
+ 'algebraic variables with the algebraic_vars_idx option')
+ self.algvaridx = self.options['algebraic_vars_idx']
self.excl_algvar_err = self.options['exclude_algvar_from_error']
self.success = 1
@@ -289,6 +293,10 @@ def init_step(self, t0, y0, yp0, y_ic0_retn = None, yp_ic0_retn = None):
def _reset(self, n, has_jac):
# Calculate parameters for Fortran subroutine ddaspk.
self.neq = n
+ self.algebraic_var = [1]*self.neq
+ if self.algvaridx is not None:
+ for ind in self.algvaridx:
+ self.algebraic_var[ind] = -1
self.info = zeros((20,), int32) # default is all info=0
self.info[17] = 2 # extra output on init cond computation
if (isscalar(self.atol) != isscalar(self.rtol)) or (
@@ -397,14 +405,21 @@ def solve(self, tspan, y0, yp0, hook_fn = None):
t_retn[0] = tinit
y_retn[0,:] = y0[:]
yp_retn[0, :] = yp0[:]
+ indbreak = None
for ind, time in enumerate(tspan[1:]):
if not self.success:
+ indbreak = ind + 1
break
- result = self.__run([2, 3], y_retn[ind], yp_retn[0], t_retn[ind], time)
+ result = self.__run([2, 3], y_retn[ind], yp_retn[ind], t_retn[ind], time)
t_retn[ind+1] = result[2]
y_retn[ind+1][:] = result[0][:]
yp_retn[ind+1][:] = result[1][:]
self.t = t_retn[-1]
+ if indbreak is not None:
+ self.t = t_retn[indbreak-1]
+ return self.success, t_retn[:indbreak], y_retn[:indbreak],\
+ yp_retn[:indbreak], t_retn[indbreak], y_retn[indbreak],\
+ yp_retn[indbreak]
return self.success, t_retn, y_retn, yp_retn, None, None, None
def __run(self, states, *args):
View
31 scikits/odes/lsodiint.py
@@ -27,11 +27,12 @@
- rband=None|int
- method='adams'|'bdf'
- with_jacobian=0|1
-- nsteps = int
+- max_steps = int
+- max_step_size = float
- (first|min|max)_step = float
- tcrit=None|float
- order = int # <=12 for adams, <=5 for bdf
-- compute_initcond = None|'yode0'
+- compute_initcond = None|'yp0'
Details:
@@ -67,11 +68,12 @@ def adda(t, y, ml, mu, p, nrowp):
- jacfn is not a supported option for lsodi.
-- compute_initcond: None or 'yode0'
+- compute_initcond: None or 'yp0'
LSODI may be able to compute the initial conditions if you do not know them
precisely and the problem is not algebraic at t=0.
- If yode0, then the differential variables (y of the ode system at time 0)
- will be used to solve for the derivatives of the differential variables.
+ If yp0, then the differential variables (y of the ode system at time 0)
+ will be used to solve for the derivatives of the differential variables,
+ so yp0 will be calculated.
Source: http://www.netlib.org/ode/lsodi.f
"""
@@ -125,8 +127,8 @@ def __init__(self, resfn, **options):
'uband': None,
'tcrit': None,
'order': 0,
- 'nsteps': 500,
- 'max_step': 0.0, # corresponds to infinite
+ 'max_steps': 500,
+ 'max_step_size': 0.0, #corresponds to infinite
'min_step': 0.0,
'first_step': 0.0, # determined by solver
'method': "adams",
@@ -163,16 +165,16 @@ def _init_data(self):
self.tcrit = self.options['tcrit']
self.order = self.options['order']
- self.nsteps = self.options['nsteps']
- self.max_step = self.options['max_step']
+ self.nsteps = self.options['max_steps']
+ self.max_step = self.options['max_step_size']
self.min_step = self.options['min_step']
self.first_step = self.options['first_step']
if re.match(self.options['method'], r'adams', re.I): self.meth = 1
elif re.match(self.options['method'], r'bdf', re.I): self.meth = 2
else: raise ValueError('Unknown integration method %s'%(self.options['method']))
if self.options['compute_initcond'] is None:
self.compute_initcond = 0
- elif re.match(self.options['compute_initcond'], r'yode0', re.I):
+ elif re.match(self.options['compute_initcond'], r'yp0', re.I):
self.compute_initcond = 1
else:
raise ValueError('Unknown init cond calculation method %s' %(
@@ -312,14 +314,21 @@ def solve(self, tspan, y0, yp0, hook_fn = None):
else:
itask = self.call_args[3]
self.call_args[3] = 4
+ intbreak = None
for ind, time in enumerate(tspan[1:]):
- result = self.__run(y_retn[ind], yp_retn[0], t_retn[ind], time)
+ result = self.__run(y_retn[ind], yp_retn[ind], t_retn[ind], time)
if not self.success:
+ intbreak = ind+1
break
t_retn[ind+1] = result[2]
y_retn[ind+1][:] = result[0][:]
yp_retn[ind+1][:] = result[1][:]
self.t = t_retn[-1]
+ if indbreak is not None:
+ self.t = t_retn[indbreak-1]
+ return self.success, t_retn[:indbreak], y_retn[:indbreak],\
+ yp_retn[:indbreak], t_retn[indbreak], y_retn[indbreak],\
+ yp_retn[indbreak]
return self.success, t_retn, y_retn, yp_retn, None, None, None
def __run(self, *args):
View
9 scikits/odes/sundials/ida.pyx
@@ -225,9 +225,12 @@ cdef class IDA:
'algebraic_vars_idx':
Values: numpy vector or None (= default)
Description:
- If given problem is of type DAE, some items of the residual vector returned by the 'resfn' have to
- be treated as algebraic variables. These are denoted by the position (index) in the residual vector.
- All these indexes have to be specified in the 'algebraic_vars_idx' array.
+ If given problem is of type DAE, some items of the residual
+ vector returned by the 'resfn' have to be treated as
+ algebraic variables. These are denoted by the position
+ (index) in the residual vector.
+ All these indexes have to be specified in the
+ 'algebraic_vars_idx' array.
'exclude_algvar_from_error':
Values: False (= default), True
Description:

0 comments on commit e826a41

Please sign in to comment.
Something went wrong with that request. Please try again.