[SciPy-User] odespy AttributeError

Johann Cohen-Tanugi johann.cohentanugi@gmail....
Thu Feb 21 10:02:26 CST 2013


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@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@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list