Integral result discrepancy between SAGE and MAPLE

[jorges]

Hi,
Some time ago I translated the derivation of a problem from MAPLE to SAGE. It took some time but I got it all working, until last week I found some strange results for a specific value of one of the parameters. Basically MAPLE and SAGE agree on all but one of the values I tried for parameter 'q'. The problem arises for 'q'=1 (See the minimal example below). I tend to trust MAPLE's result, given the meaning of the integral, a normalized form of the power dissipated at a switch when it is off, i.e V^2/R
Any clues of what might be going on here?

This is with SAGE 6.1.1, on Arch linux 64 bits


# Note that d, q, mon, moff, Vdd and omega are independent INPUT parameters, and that the rest are dependent intermediate results
var('d p q phi mon moff Vdd Coff1 Coff2 omega')
params = {'phi': -0.342466921443868, 'd': 1, 'p': 3.69012461123172, 'q': 1.0, 'Vdd': 1, 'Coff2': 621.040280393254 - 1743.07649791142*I, 'moff': 1000, 'Coff1': 621.040280393254 + 1743.07649791141*I, 'omega': 1, 'mon': 0.01}
aoff =(1/2)*(-1+sqrt(1-4*q^2*moff^2))/moff
boff = -(1/2)*(1+sqrt(1-4*q^2*moff^2))/moff
TVcoff(t)=Coff1*exp(aoff*omega*t)+Coff2*exp(boff*omega*t)+Vdd+Vdd*(p*q^2*moff^2*(q-1)*(q+1)*cos(omega*t+phi)+p*sin(omega*t+phi)*moff*q^2)/(1+(q^4-2*q^2+1)*moff^2);
ps=integral(((TVcoff(t))^2/moff),t, d*pi/omega, 2*pi/omega)
Nps = ps.subs(**params).n()
# The value I am interested in is the real part of Nps
print Nps

-1083.52083510533 + 4.70489703060889*I

MAPLE returns

0.168006953631138e-1-0.615350612025753e-19*I



jorge

[John H Palmieri]

I'm guessing that the issue is that your integrand simplifies when q=1 -- at least one term becomes zero -- but you do the integral before doing this simplification. Maybe the symbolic integration is not valid when q=1. If you plug in the parameters before integrating, you get something very different from Sage. Adding these two lines

qs=integral(((TVcoff(t))^2/moff).subs(**params),t, d*pi/omega, 2*pi/omega) # substitute before integrating
Nqs = qs.subs(**params).n()

gives me something close to Maple's answer:

sage: Nps
1057.74513808638 + 3.42074504083530*I
sage: Nqs
0.0164729319512844 + 5.58793544769287e-9*I

   John

[jorges]

On Wednesday, 16 April 2014 01:24:25 UTC-3, John H Palmieri wrote:

I'm guessing that the issue is that your integrand simplifies when q=1 -- at least one term becomes zero -- but you do the integral before doing this simplification. Maybe the symbolic integration is not valid when q=1. If you plug in the parameters before integrating, you get something very different from Sage. Adding these two lines

qs=integral(((TVcoff(t))^2/moff).subs(**params),t, d*pi/omega, 2*pi/omega) # substitute before integrating
Nqs = qs.subs(**params).n()

gives me something close to Maple's answer:

sage: Nps
1057.74513808638 + 3.42074504083530*I
sage: Nqs
0.0164729319512844 + 5.58793544769287e-9*I

Indeed, that looks good. I still don't understand what you mean by "the symbolic integration is not valid when q=1". I would think substituting before or after should not make a difference. Also, I remember that I attempted to simplify before integrating before, mainly to speed up calculations, but had some bad experience with real numbers and maxima. I can't recall the precise details right now, but after that I "learned" to defer substitution to the very end. That is until now. I love SAGE, but these things, i.e. having to know special cases and how to handle them, make using it much more difficult, especially for a non-power user as me.
Thanks for pointing out how to avoid this issue.

Jorge

[John H Palmieri]

On Tuesday, April 15, 2014 9:58:51 PM UTC-7, jorges wrote:

On Wednesday, 16 April 2014 01:24:25 UTC-3, John H Palmieri wrote:

I'm guessing that the issue is that your integrand simplifies when q=1 -- at least one term becomes zero -- but you do the integral before doing this simplification. Maybe the symbolic integration is not valid when q=1. If you plug in the parameters before integrating, you get something very different from Sage. Adding these two lines

qs=integral(((TVcoff(t))^2/moff).subs(**params),t, d*pi/omega, 2*pi/omega) # substitute before integrating
Nqs = qs.subs(**params).n()

gives me something close to Maple's answer:

sage: Nps
1057.74513808638 + 3.42074504083530*I
sage: Nqs
0.0164729319512844 + 5.58793544769287e-9*I

Indeed, that looks good. I still don't understand what you mean by "the symbolic integration is not valid when q=1". I would think substituting before or after should not make a difference.

It's possible, for example, that when doing the integration, it divides by q-1 to "simplify" some algebra. Maybe it does something else which turns out not to be valid when q=1. This is all just a guess, anyway.

  John

[John H Palmieri]

By the way, you can also do

    numerical_integral(((TVcoff(t))^2/moff).subs(**params).real(), pi, 2*pi)
    numerical_integral(((TVcoff(t))^2/moff).subs(**params).imag(), pi, 2*pi)

This function seems to require a real function, so I've put in .real() and .imag() to break it into pieces.

--
John

 

[Robert Dodier]

On 2014-04-16, jorges <jscand...@gmail.com> wrote:

> Indeed, that looks good. I still don't understand what you mean by "the
> symbolic integration is not valid when q=1". I would think substituting
> before or after should not make a difference.

Well, if they're not the same, that's a bug in Maxima. If you have time,
it would be terrific if you could make a bug report.
See: http://sf.net/p/maxima/bugs

> Also, I remember that I
> attempted to simplify before integrating before, mainly to speed up
> calculations, but had some bad experience with real numbers and maxima. I
> can't recall the precise details right now, but after that I "learned" to
> defer substitution to the very end.

A different way to handle that is to change the floats to rational numbers,
via ratsimp. Then Maxima will be happy with expressions which contain
those numbers.

Hope this helps,

Robert Dodier

[jorges]

Thanks John and Robert for your suggestions. I'll definitely fill a bug against maxima