[SciPy-User] odespy AttributeError

[Mara Grahl]

Hi,

thank you for the suggestion, I'll consider that,

best
Mara


-------- Original-Nachricht --------
> Datum: Thu, 21 Feb 2013 17:02:26 +0100
> Von: Johann Cohen-Tanugi <johann.cohentanugi at gmail.com>
> An: SciPy Users List <scipy-user at scipy.org>
> CC: MaraMaus at nurfuerspam.de
> Betreff: Re: [SciPy-User] odespy AttributeError

> hello, best is probably to contact the main developer directly, as this 
> from github does not look like a community effort : 
> https://github.com/hplgit
> 
> good luck,
> Johann
> 
> On 02/21/2013 04:39 AM, MaraMaus at nurfuerspam.de wrote:
> > Hi,
> >
> > since there doesn't seem to be a mailing list for the package odespy, I
> place my question here - I hope this is ok?
> >
> > Whereas the examples from the odespy manual work fine, I encounter a
> strange error when trying to use odespy for some complicated system of
> ordinary differential equations:
> >
> > Traceback (most recent call last):
> > File "ode1.py", line 131, in <module>
> >      u, s = solver.solve(time_points)
> > File "C:\Python27\lib\site-packages\odespy\solvers.py", line 1036, in
> solve
> >      self.u[n+1] = self.advance()   # new value
> > File "C:\Python27\lib\site-packages\odespy\RungeKutta.py", line 194, in
> advance
> >      rms_norm = np.sqrt(np.sum(rms*rms)/self.neq)
> > AttributeError: sqrt
> >
> > Does anyone know what might be the reason for this error message?
> >
> > I would be glad for any help, thank you in advance,
> >
> > Mara
> >
> >
> >
> > For completeness the full code:
> > A short description: most of the code (between the hashs #setup_beg and
> #setup_end) sets up the system of 9 ordinary 1st order differential
> equations for the variables v1, v2, v3, vp1, vp2, vp3, vpp1, vpp2, vpp3. The list
> DLHS contains explicit expressions for
> > [diff(v1,s), diff(v2,s), ... , diff(vpp3,s)]
> >
> >
> > #setup_beg
> > nx=3
> > T=45.1
> > mu=253.9
> > sL=500.0
> > sR=6.0
> > xL=0.0
> > xR=100.0**2
> > dx=(xR-xL)/(nx-1)
> >
> > x={}
> > for i in range(1,nx+1):
> >   x[i]=xL+dx*(i-1)
> >
> > print(x[nx])
> >
> > # vp is first derivative of v, vp1 is first derivative of v at grid
> point 1 etc.
> > # note that vp[0]=vp1, ..., vp[nx-1]=vpnx
> >
> > from sympy import solve, Symbol
> >
> > pts=range(1,nx+1)
> > vp=[Symbol('vp'+str(i)) for i in pts]
> > vpp=[Symbol('vpp'+str(i)) for i in pts]
> > vppp=[Symbol('vppp'+str(i)) for i in pts]
> > vpppp=[Symbol('vpppp'+str(i)) for i in pts]
> > #beachte: vp[0]=vp1,...,vp[nx-1]=vpnx etc.
> >
> > eq1={}
> > for i in range(1,nx):
> >   eq1[i]=vp[i-1]+vpp[i-1]*(x[i+1]-x[i])/2 +
> vppp[i-1]*(x[i+1]-x[i])**2/8+vpppp[i-1]*(x[i+1]-x[i])**3/48 \
> >      -vp[i]+vpp[i]*(x[i+1]-x[i])/2 -
> vppp[i]*(x[i+1]-x[i])**2/8+vpppp[i]*(x[i+1]-x[i])**3/48
> >
> > eq2={}
> > for i in range(1,nx):
> >   eq2[i]=vpp[i-1]+vppp[i-1]*(x[i+1]-x[i])/2 +
> vpppp[i-1]*(x[i+1]-x[i])**2/8 \
> >      -vpp[i]+vppp[i]*(x[i+1]-x[i])/2 - vpppp[i]*(x[i+1]-x[i])**2/8
> >
> > eq3={}
> > eq3[1]=vppp[0]+vpppp[0]*(x[2]-x[1])/2-vppp[1]-vpppp[1]*(x[1]-x[2])/2
> >
> eq3[2]=vppp[nx-2]+vpppp[nx-2]*(x[nx]-x[nx-1])/2-vppp[nx-1]-vpppp[nx-1]*(x[nx-1]-x[nx])/2
> >
> >
> > eqs=[]
> > for i in range(1,nx):
> >   eqs.append(eq1[i])
> > for i in range(1,nx):
> >   eqs.append(eq2[i])
> > eqs.append(eq3[1])
> > eqs.append(eq3[2])
> >
> > vars=[]
> > for i in range(1,nx+1):
> >   vars.append(vppp[i-1])
> > for i in range(1,nx+1):
> >   vars.append(vpppp[i-1])
> >
> > sol=solve(eqs,vars)
> >
> > #beachte: vppp[0]=vppp1,...,vppp[nx-1]=vpppnx etc.
> > # Check z.B. fuer nx=5: (gleiche Loesung mit MM)
> > # vpppp5=
> > #print(sol[vpppp[4]])
> > # vpppp1=
> > #print(sol[vpppp[0]])
> >
> > from sympy import Function, diff
> > from sympy.functions import coth, tanh, sqrt
> >
> > s=Symbol('s')
> > vh=Function('vh')
> > xh=Symbol('xh')
> >
> > Epi= sqrt( s**2 + 2*diff(vh(xh),xh) )
> > Esi= sqrt( s**2 + 2*diff(vh(xh),xh) +4*xh*diff(vh(xh),xh,xh) )
> > Eq= sqrt(s**2 +3.2**2*xh)
> >
> > import math
> > pi=math.pi
> >
> > gl1= s**4/12/pi**2*( 3/Epi*coth(Epi/2/T) +1/Esi*coth(Esi/2/T)
> -12/Eq*(tanh( (Eq-mu)/2/T ) + tanh( (Eq+mu)/2/T ) ) )
> > gl2= diff(gl1,xh)
> > gl3= diff(gl2,xh)
> >
> > DLHS=[]
> >
> > for i in range(1,nx+1):
> >   DLHS.append( gl1.subs({'Derivative(vh(xh), xh,
> xh)':vpp[i-1],'Derivative(vh(xh), xh)':vp[i-1],'xh':x[i]}).subs(sol) )
> > for i in range(1,nx+1):
> >   DLHS.append( gl2.subs({'Derivative(vh(xh), xh, xh,
> xh)':vppp[i-1],'Derivative(vh(xh), xh, xh)':vpp[i-1],'Derivative(vh(xh),
> xh)':vp[i-1],'xh':x[i]}).subs(sol) )
> > for i in range(1,nx+1):
> >   DLHS.append( gl3.subs({'Derivative(vh(xh), xh, xh, xh,
> xh)':vpppp[i-1],'Derivative(vh(xh), xh, xh, xh)':vppp[i-1],'Derivative(vh(xh), xh,
> xh)':vpp[i-1], \
> >                         'Derivative(vh(xh),
> xh)':vp[i-1],'xh':x[i]}).subs(sol) )
> >
> >
> #DLHS[0]=D[v(x_1),s],...,DLHS[nx-1]=D[v(x_nx),s],DLHS[nx]=D[vp(x_1),s],..., DLHS[nx+nx+nx-1]=D[vpp(x_nx),s]
> >
> >
> > #setup_end
> >
> >
> > v=[Symbol('v'+str(i)) for i in pts]
> >
> > #initial conditions
> > isc=[]
> >
> > for i in range(1,nx+1):
> >   isc.append( 5/2*x[i]**2)
> > for i in range(1,nx+1):
> >   isc.append( 5*x[i])
> > for i in range(1,nx+1):
> >   isc.append( 5)
> >
> >
> >
> >
> > # to do: generalize to arbitrary nx
> > def f(u,s):
> >      v1, v2, v3, vp1, vp2, vp3, vpp1, vpp2, vpp3 = u
> >      return DLHS
> >
> > import odespy
> > solver = odespy.RungeKutta.CashKarp(f)
> > solver.set_initial_condition(isc)
> > from numpy import linspace
> > T = 490 # end of simulation
> > N = 30 # no of time steps
> > time_points = linspace(sL, T, N+1)
> > u, s = solver.solve(time_points)
> >
> > #from matplotlib.pyplot import *
> > #first=u[:,0]
> > #plot(s, first)
> > #show()
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User at scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user