| Re: Integration qpint [message #81057] |
Mon, 13 August 2012 05:30 |
Gompie
Messages: 76
Registered: August 2012
|
Member |
|
|
Hi Craig , Helder
Thanks a lot for helping me out with this incredibly useful function.
Gompie
|
|
|
|
| Re: Integration qpint [message #81059 is a reply to message #81057] |
Mon, 13 August 2012 02:01  |
Helder Marchetto
Messages: 520
Registered: November 2011
|
Senior Member |
|
|
On Sunday, August 12, 2012 7:42:57 PM UTC+2, Gompie wrote:
> Hi
>
> Can anyone help. I am using qpint1d to integrate but i get
>
> "Program caused arithmetic error: Floating overflow:"
>
>
>
> The code is
>
> C1 = 1.1910427E-12
>
> c2=1.43883
>
> T=300
>
> c=c1*(c2^(-4))*T^4
>
> llimit=c2*2380/T
>
> ulimit=c2*2940/T
>
> print,"X=",llimit,ulimit
>
> lin= c*qpint1d('x^3/(EXP(X)-1)', llimit, +inf, /expr)
>
> rin=c*qpint1d('x^3/(EXP(X)-1)', ulimit, +inf, /expr)
>
> print,abs(lin-rin)
>
> end
>
>
>
> This outputs a number ( i am not sure if it is correct) but it also says floating point error.
>
> Thanks in advance
>
> Gompie
Hi Gompie,
I didn't go through the whole of your code, but you define T as an integer.
Try rerunning your code substituting this to the line where you define T:
T=300.0
Doing this, the value of c changes from 3.48598e-009 (T is integer) to 0.00225099 (T is floating point).
Maybe also change the definitions of llimit and ulimit to:
llimit=c2*2380.0/T
ulimit=c2*2940.0/T
This might help. But I have not tried the integration.
Cheers,
h
|
|
|
|
| Re: Integration qpint [message #81060 is a reply to message #81059] |
Mon, 13 August 2012 02:00 |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
On Sunday, August 12, 2012 1:42:57 PM UTC-4, Gompie wrote:
> Hi
>
> Can anyone help. I am using qpint1d to integrate but i get
>
> "Program caused arithmetic error: Floating overflow:"
>
>
>
> The code is
>
> C1 = 1.1910427E-12
>
> c2=1.43883
>
> T=300
>
> c=c1*(c2^(-4))*T^4
>
> llimit=c2*2380/T
>
> ulimit=c2*2940/T
>
> print,"X=",llimit,ulimit
>
> lin= c*qpint1d('x^3/(EXP(X)-1)', llimit, +inf, /expr)
>
> rin=c*qpint1d('x^3/(EXP(X)-1)', ulimit, +inf, /expr)
>
> print,abs(lin-rin)
>
> end
>
>
>
> This outputs a number ( i am not sure if it is correct) but it also says floating point error.
Check the status variable. It's not an error. You can verify this yourself by
IDL reports *any* overflow or underflow that occurs. That's not necessarily bad.
What do you think happens when you try to evaluate EXP(+inf)? Or EXP(+3000)? Your user function generates an overflow, that's what happens. IDL is warning you that it is treating EXP(+3000) equivalent to +infinity. That is what happens when IDL cannot evaluate the full range of your user function with full precision. EXP(+3000) simply cannot be represented on a modern computer.
If you don't want overflow errors, then rewrite your expression to use EXP(-X) instead. Of course, then you will get underflow errors. Life sucks, doesn't it? :-)
If you really don't want overflow errors for this function, then you must avoid using +inf as your upper bound.
Craig Markwardt
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
