[SciPy-User] odespy AttributeError

[MaraMaus at nurfuerspam.de]

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()