[SciPy-user] odeint -- segmentation fault

Giovanni Samaey Giovanni.Samaey at cs.kuleuven.ac.be
Tue Aug 31 13:27:41 CDT 2004


Hi,

I have *much* more information now.
A simple Python rhs function suffices to get the segmentation fault.
A sample code that produces the segfault is attached:
The function integrate calls odeint with the "equation being defined through
dydt.  This dydt initialized another integration, around each point x, 
and calls
a function step, which executes odeint with "equation" rhs.

The crucial point to get the segfault is that rhs has a different number 
of odes than dydt.
The sample code added here gives print statements that should help 
understand execution,
and a comment is added that would indicate the segfault.

I am working on linux with SciPy version '0.3.0_266.4239', Numeric 23.3 
and python 2.3.4


Robert Kern wrote:

> Giovanni Samaey wrote:
>
> [snip]
>
>> Since I do not have access to odeint's code, I cannot tell what is 
>> going on, nor how to solve this...
>>
>> How could this problem be solved effectively?
>
>
> First, try an example where rhs is not a SWIGed function. That might 
> be the problem itself. I tried a simple example with a Python rhs (on 
> a Powerbook with CVS SciPy) and I didn't get a segfault.
>
> If you want to see odeint's source code to help you debug, just 
> download it: http://www.scipy.org/download/
>
> When you have more information, also please provide your platform and 
> SciPy version so we know where to look.
>
> Thanks.
>


-- 
Giovanni Samaey		 	http://www.cs.kuleuven.ac.be/~giovanni/ 
Katholieke Universiteit Leuven 	      email: giovanni at cs.kuleuven.ac.be 
Departement Computerwetenschappen                  phone: +32-16-327081
Celestijnenlaan 200A, B-3001 Heverlee, Belgium       fax: +32-16-327996
Office: A04.36


-------------- next part --------------
import scipy
from scipy.integrate import odeint,simps

def step(y0,t,x):
    print "In step!"
    print "y0: ",scipy.shape(y0)
    y0=scipy.ravel(y0)
    y=odeint(rhs,y0,scipy.array([t,t+1e-5]),ml=1,mu=1,args=(x,))
    print "y: ", scipy.shape(y)
    ret=scipy.reshape(y[-1,:],(len(y0),))
    print "ret: ", scipy.shape(ret)
    return scipy.reshape(y[-1,:],(len(y0),))
  
def rhs(y,t,x):
    print "In rhs!"
    ny=scipy.size(y)
    ydot=scipy.zeros(ny,scipy.Float)
    ydot[1:-1]=y[2:]-2*y[1:-1]+y[:-2]
    return ydot

def dydt(y,t,x):
    print "In dydt!"
    y=scipy.reshape(y,(1,len(y)))
    print "y: ", scipy.shape(y)
    ynew=scipy.zeros(scipy.shape(y))
    for i in scipy.arange(len(x)):
        print i
        xx=scipy.arange(x[i]-5e-3,x[i]+5e-3,1e-7)
        yy=scipy.ones(scipy.shape(xx),scipy.Float)
        yynew=step(yy,t,xx)
        print "yynew: ", scipy.shape(yynew)
        ynew[0,i]=yynew[500]
    print "ynew: ",scipy.shape(y)
    print "dydt shape: ", scipy.shape(scipy.ravel((ynew-y)/1e-5))
    # the segmentation fault occurs when this return statement is executed!
    return scipy.ravel((ynew-y)/1e-5)

def integrate(y,t,x):
    print "In integrate!"
    print "Initial condition shape: ", scipy.shape(scipy.ravel(y))
    return odeint(dydt,scipy.ravel(y),t,ml=1,mu=1,args=(x,))

if __name__=="__main__":

    x=scipy.arange(0,1.05,0.1)
    y0=1-4*(x-0.5)**2
    result=integrate(y0,t=scipy.arange(0,1e-2,5e-3),x=x)


    
    
    

    

        
    

        
        
    
        




More information about the SciPy-user mailing list