problem with numerical integration and differentiation with scipy

Hi,

I was trying to learn Scipy, using it for integration and differentiation, but at the very initial step I encontered the following problems. First, the following integration has a complex value as its numerical result when is solved analytically, but the numerical solution only gives its real part without any note for caution:

reset()
forget()

import scipy as sp
from scipy import integrate

var('y')
f=lambda x: 10^10*sin(x)
g=lambda x,y: f(x+y^2)
I=integrate.dblquad(g(x,y), 0, pi, lambda y:0,lambda y:1)[0]; show(I)
show( integral(integral(g(x,y),x,0,1),y,0,pi) )
show( integral(integral(g(x,y),x,0,1),y,0,pi).n() )


is this natural or I am doing something wrong that may cause my even Real calculations erratic? Next, about numerical differentiation, it seems that the only function that works for callable functions is scipy.derivative() if I'm right!? However, I couldn't work with it: 1st) when I am not going to specify the point at which the differentiation is to be taken, e.g. when the differentiation is under an integral so that it is the integral that should assign the numerical values to its integrand's variable, not me. As a simple example I tried this code:

import scipy as sp
from scipy import integrate, derivative
I2=integrate.quad( sp.derivative(f(y),y, dx=0.00001, n=1, order=7) , 0, pi,)[0]; show(I2)
show( integral(diff(f(y),y),y,0,1).n() )


also it gives the warning that "Warning: The ocurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated."

2nd) when I want to find the derivative of a multivariable function like g(x,y) in the above example, and something like sp.derivative(g(x,y),(x,0.5), dx=0.01, n=1, order=3) gives error as was easily expected.

Looking forward to hearing from you if that integration result is natural and how to resolve the above cited problems with numerical differentiation.

Best Regards

edit retag close merge delete

I feel like maybe you should be asking this on a Scipy list... I don't see any Sage in here, though hopefully someone here will know Scipy well enough to answer nonetheless!

( 2012-08-28 05:35:50 -0600 )edit

Yes it is not actually a Sage question, but I am using Scipy from within Sage --with the hope that a code written in symbolic language, being way too slow to run, can be used in numerical calculations modified to some very small extents, as well. Thanks anyway for you cared, I also ask stack overflow as you suggested

( 2012-08-28 12:36:35 -0600 )edit

Sort by » oldest newest most voted

This may not answer your question, but this Stackoverflow question could be relevant.

more

thanks, well this says "scipy.integrate.quad requires the integrated function return floats (aka real numbers) for the algorithm it uses" and my function f(y)=sin(y) is a real function, however, yes I guess sin(y) might be translated to its complex form [exp(iy)-exp(-iy)] / 2 before integration and this is the cause for such a problem! Reasonable! Thanks again .

( 2012-08-28 12:22:54 -0600 )edit