Error in Mathematica

  • #1
EngWiPy
1,367
61
Hello,

I have the following code in Mathematica, and it gives the following error:

Code:
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 12 recursive bisections in x near {x} = {0.000156769}. NIntegrate obtained 0.21447008480474927` and 5.637666621985554`*^-13 for the integral and error estimates. >>

How to solve this error?

Code:
yp = 10^(5/10);
GSS = 10;

For[yQdB = -10, yQdB <= 10, yQdB++;
yQ = 10^(yQdB/10);
Print[yQ];
A1 = NIntegrate[
   1/x^2*Exp[(GSS*(x + GSS*yQ))/(x*yp) - x]*
    ExpIntegralE[1, (GSS*(x + GSS*yQ))/(x*yp)], {x, LL, Infinity},
   PrecisionGoal -> 12, MaxRecursion -> 12];
]
 

Answers and Replies

  • #2
36,026
12,926
You have a loop, so a natural step to narrow down the error would be to test individual iterations of the loop. Does every value of yQdB produce the problem? Did you try a lower precision goal or a deeper recursion? Note that the precision is close to the target.
 
  • #3
EngWiPy
1,367
61
No actually it's not for every value of yQdB. What do you mean by "Note that the precision is close to the target."?
 
  • #4
36,026
12,926
NIntegrate obtained 0.21447008480474927` and 5.637666621985554`*^-13
That is a relative uncertainty of about 3*10-12.
 
  • #5
EngWiPy
1,367
61
So, do I need to increase the precision?
 
  • #6
EngWiPy
1,367
61
Another thing, which is that the lower limit of the integral is theoretically zero. However, since the program doesn't compute it by saying that the integral has evaluated to overflow, indeterminate, or infinity, I replaced it but ##LL=10^{-5}##. However, when I compare the results with Monte-Carlo simulations (the above formula is the numerical result), there is some discrepancy. How to make the numerical evaluation as accurate as possible?
 
  • #7
36,026
12,926
So, do I need to increase the precision?
The opposite. Mathematica is unable to require the precision you ask for with the number of iterations you allow. Allow more iterations or reduce the precision you ask for and it should work.
However, when I compare the results with Monte-Carlo simulations (the above formula is the numerical result), there is some discrepancy. How to make the numerical evaluation as accurate as possible?
Why do you think the problem is with Mathematica?
 
  • #8
EngWiPy
1,367
61
Fair point. I'm not sure if the problem is with Mathematica. I posed a question in the calculus forums to double check if there is something fundamentally wrong with the math, but I haven't gotten a conclusive answer.
 
  • #9
36,026
12,926
I posed a question in the calculus forums to double check if there is something fundamentally wrong with the math
Well, we can't tell. You wrote down an integral. Mathematica evaluates it. Is this the right integral you want to calculate? No idea.
 
  • #10
JorisL
491
189
Also it helps to give the integral you are trying to evaluate in latex with proper symbols.
That makes it easier to check for typos in the code (don't ask me how many times I made that kind of mistake :S) or whether there is some mathematical reason for the error message.
 
  • #11
EngWiPy
1,367
61
I think the mathematical derivations were accurate. I'm not sure about the limits though. The integral is the following

[tex]\int_0^{\infty}\frac{1}{x^2}\exp\left(\frac{A}{x}-x\right)E_1\left(B+\frac{A}{x}\right)\,dx[/tex]

More details are found here.
 
  • #12
JorisL
491
189
You mentioned in post #3 that the error doesn't happen for every yQdB, could you narrow it down some more?
With these kinds of expressions you sometimes overlook undefined cases.

Edit;
I also think your expression in mathematica has ##A=A(x)## meaning you have an extra subtlety.
 
Last edited:
  • #13
EngWiPy
1,367
61
Really? Do you mean I should use square brackets all the time in Mathematica?
 
  • #14
JorisL
491
189
I have rewritten the expression in your code to make it more clear what happens.
$$
\frac{1}{x^2}\exp\left[ \frac{G_{SS}\left(x+G_{SS} \cdot y_Q\right)}{xy_p}-x\right]
E_1\left[\frac{G_{SS}\left(x+G_{SS}y_Q\right)}{xy_p}\right]=\frac{1}{x^2}\exp\left[ \frac{G_{SS}}{y_p}+\frac{G_{SS}^2 \cdot y_Q}{xy_p}-x\right]E_1\left[\frac{G_{SS}}{y_p}+\frac{G_{SS}^2y_Q}{y_px}\right]$$

If you now define ##A=\frac{G_{SS}^2y_Q}{y_p}##, ##B=\frac{G_{SS}}{y_p}## and ##C=\frac{G_{SS}}{y_p}## this results in

$$
\frac{1}{x^2}\exp\left[ \frac{G_{SS}\left(x+G_{SS} \cdot y_Q\right)}{xy_p}-x\right]
E_1\left[\frac{G_{SS}\left(x+G_{SS}y_Q\right)}{xy_p}\right] = \frac{1}{x^2}\exp\left[C +\frac{A}{y_p}-x\right]E_1\left[\frac{G_{SS}}{y_p}+\frac{A}{x}\right]
$$

So either you have different A's as you can get the right form in the exponential by setting ##A^\prime(x)=\frac{G_{SS}\left(x+G_{SS} \cdot y_Q\right)}{y_p}## leading to

$$\frac{1}{x^2}\exp\left[\frac{A^\prime(x)}{x}-x\right]E_1\left[\frac{G_{SS}}{y_p}+\frac{A}{x}\right]$$

So you have to check that if really want this or the integral in post #11.

Edit;
I used square brackets and parentheses for clarity but these expressions are not code, I don't like reading the code directly most of the time as it can obscure things which seem obvious in retrospect.

Edit 2;
As Tom pointed out I made a mistake, it's now fixed
 
Last edited:
  • #15
Tom.G
Science Advisor
Gold Member
4,463
3,231
I have rewritten the expression...

1x2exp&#x2061;[GSS(x+GSS&#x22C5;yQ)xyp&#x2212;x]E1[GSS(x+GSSyQ)xyp]=1x2exp&#x2061;[C+Ayp&#x2212;x]E1[GSSyp+Ax]" style="font-size: 106%; position: relative;" tabindex="0" class="mjx-chtml MathJax_CHTML" id="MathJax-Element-8-Frame">1x2exp[C+Ayp−x]E1[GSSyp+Ax]

As you you Defined A, above, should the exp term in the rewrite have A/x rather than A/yp ?

EDIT: well, that didn't copy thru as intended!
 
  • #16
JorisL
491
189
You're absolutely right, thanks for pointing that out.

Best way to copy formulas is to use quote or reply since that copies the LaTeX codes.
 
  • #17
EngWiPy
1,367
61
I have rewritten the expression in your code to make it more clear what happens.
$$
\frac{1}{x^2}\exp\left[ \frac{G_{SS}\left(x+G_{SS} \cdot y_Q\right)}{xy_p}-x\right]
E_1\left[\frac{G_{SS}\left(x+G_{SS}y_Q\right)}{xy_p}\right]=\frac{1}{x^2}\exp\left[ \frac{G_{SS}}{y_p}+\frac{G_{SS}^2 \cdot y_Q}{xy_p}-x\right]E_1\left[\frac{G_{SS}}{y_p}+\frac{G_{SS}^2y_Q}{y_px}\right]$$

If you now define ##A=\frac{G_{SS}^2y_Q}{y_p}##, ##B=\frac{G_{SS}}{y_p}## and ##C=\frac{G_{SS}}{y_p}## this results in

$$
\frac{1}{x^2}\exp\left[ \frac{G_{SS}\left(x+G_{SS} \cdot y_Q\right)}{xy_p}-x\right]
E_1\left[\frac{G_{SS}\left(x+G_{SS}y_Q\right)}{xy_p}\right] = \frac{1}{x^2}\exp\left[C +\frac{A}{y_p}-x\right]E_1\left[\frac{G_{SS}}{y_p}+\frac{A}{x}\right]
$$

So either you have different A's as you can get the right form in the exponential by setting ##A^\prime(x)=\frac{G_{SS}\left(x+G_{SS} \cdot y_Q\right)}{y_p}## leading to

$$\frac{1}{x^2}\exp\left[\frac{A^\prime(x)}{x}-x\right]E_1\left[\frac{G_{SS}}{y_p}+\frac{A}{x}\right]$$

So you have to check that if really want this or the integral in post #11.

Edit;
I used square brackets and parentheses for clarity but these expressions are not code, I don't like reading the code directly most of the time as it can obscure things which seem obvious in retrospect.

Edit 2;
As Tom pointed out I made a mistake, it's now fixed

Yes, you are right. But if you saw the Mathematica code, it's detailed. I used the constants A and B in the mathematical integral for convenience.
 
  • #18
JorisL
491
189
So my question is, does the variable A in the mathematical integral really represent a constant?

Because that is not reflected in your mathematica code as I've shown.

Is there a unique value of yQdB for which the evaluation fails? Or are there several values that cause trouble?

@Tom.G
Copying equations works best through reply or quote. They are typeset using LaTeX through MathJax.
You seem to have copied the html-code representing the equation which is quite odd.
 
  • #19
EngWiPy
1,367
61
In the Mathematica code there is no A. As I said, I used A and B for mathematical convenience, and your expression is more accurate mathematically. But again Mathematica code is still correct as I wrote it detailed; I didn't write it in terms of A and B.

The problem occurs at almost all values of yQdB where the lower limit of the integral is 0. But when I use a lower limit that is close to zero (like in the following code), I will have the error occurring at certain values.

Code:
yp = 10^(5/10);
GSS = 1;
Cons = 10^-7;


For[yQdB = -10, yQdB <= 15, yQdB++;
yQ = 10^(yQdB/10);
A1 = 0.5*((GSS^2)*yQ )/yp*
   NIntegrate[
    1/x^2*Exp[(GSS*(x + GSS*yQ))/(x*yp) - x]*
     ExpIntegralE[1, (GSS*(x + GSS*yQ))/(x*yp)], {x, Cons, Infinity},
    PrecisionGoal -> 5, MaxRecursion -> 12];
Print[A1]
]

The above code is working fine. But when I change GSS to 100, the error appears again, but at certain values. I'm not sure why?
 
  • #20
36,215
8,197
yp = 10^(5/10);
Does Mathematica evaluate this as the square root of 10? Many programming languages would evaluate 5/10 as 0, using integer division rather than floating point division, thereby setting yp to 1.

Out of curiosity, why are you raising 10 to the 5/10 power?
 
  • #21
JorisL
491
189
In the Mathematica code there is no A. As I said, I used A and B for mathematical convenience, and your expression is more accurate mathematically. But again Mathematica code is still correct as I wrote it detailed; I didn't write it in terms of A and B.

Even if you don't have it explicitly, I showed in post #14 that you're mathematica code doesn't compute an integral of the form you gave.

Try putting this in your loop

Code:
Print["A1 = ", 0.5*((GSS^2)*yQ)/yp, "*NIntegrate[  1/x^2*Exp[ ( ", GSS, "*(x +", GSS*yQ, "))/(x*", yp, ") - x]* ExpIntegralE[1, (", GSS, "*(x + ", GSS*yQ, "))/(x*", yp, ")]"]

It will show you which integral is being evaluated (without the domain and options as they don't change).
That way you'll know what the exact integral is that doesn't work.
 
  • #22
EngWiPy
1,367
61
Does Mathematica evaluate this as the square root of 10? Many programming languages would evaluate 5/10 as 0, using integer division rather than floating point division, thereby setting yp to 1.

Out of curiosity, why are you raising 10 to the 5/10 power?

Because usually we choose the value of yp in dB, while in computation we use it in linear scale. So, yp in dB is 5 while in linear scale it's 10^(5/10).
 
  • #23
EngWiPy
1,367
61
Even if you don't have it explicitly, I showed in post #14 that you're mathematica code doesn't compute an integral of the form you gave.

Try putting this in your loop

Code:
Print["A1 = ", 0.5*((GSS^2)*yQ)/yp, "*NIntegrate[  1/x^2*Exp[ ( ", GSS, "*(x +", GSS*yQ, "))/(x*", yp, ") - x]* ExpIntegralE[1, (", GSS, "*(x + ", GSS*yQ, "))/(x*", yp, ")]"]

It will show you which integral is being evaluated (without the domain and options as they don't change).
That way you'll know what the exact integral is that doesn't work.

In mathematical form, the integrad is

[tex]\frac{1}{x^2}\exp\left[ \frac{G_{SS}\left(x+G_{SS} \cdot y_Q\right)}{xy_p}-x\right]
E_1\left[\frac{G_{SS}\left(x+G_{SS}y_Q\right)}{xy_p}\right]
[/tex]

In Mathematica it's

Code:
1/x^2*Exp[(GSS*(x + GSS*yQ))/(x*yp) - x]*
     ExpIntegralE[1, (GSS*(x + GSS*yQ))/(x*yp)]

They both look the same to me!! Am I missing something here?

I copied your line, but didn't know how to use, or what I'm looking for. :oldconfused:
 
Last edited:
  • #24
JorisL
491
189
I think the mathematical derivations were accurate. I'm not sure about the limits though. The integral is the following

[tex]\int_0^{\infty}\frac{1}{x^2}\exp\left(\frac{A}{x}-x\right)E_1\left(B+\frac{A}{x}\right)\,dx[/tex]

More details are found here.

If you compare the two you see that in the code ##B=0## while ##A## explicitly contains ##x## (it's linear in x).
Meaning it isn't the same integral. I'm sure your code is right but the mathematical discussion in the other thread isn't really relevant.

<snip>
I copied your line, but didn't know how to use, or what I'm looking for. :oldconfused:

You put it inside your For loop above the NIntegrate-function.
It will print your command without executing it i.e. show you what the integrand looks like.
 
  • #25
EngWiPy
1,367
61
But again you referred to the integral which was written in terms of A and B. I just wrote it this way here to make it easier for the reader to read. The integral as it's in post #23 is the same as the integral in Mathematica. This is the integral in my calculations. Not the one I wrote here which you just quoted.

Your line gave me this for the first value of yQdB

Code:
A1 = 199.054*NIntegrate[  1/x^2*Exp[ ( 100*(x +10 10^(1/10)))/(x*Sqrt[10]) - x]* ExpIntegralE[1, (100*(x + 10 10^(1/10)))/(x*Sqrt[10])]
 
  • #26
JorisL
491
189
All I'm saying is that they aren't the same. And you gave that expression when asked to reproduce the integral you where trying to evaluate.
This is confusing (at best).

Apparently I made a mistake when copying, it isn't important though. (The constant factor before NIntegrate shouldn't be there according to post #1, it shows up in #19 though)

Am I correct with saying ##GSS=100## and ##yQdB=1## since that would mean there's a problem with " 10 10^(1/10)" in the output.
I cannot reproduce the prefactor either (everything is off by a factor of 10)

Edit;
Can you post the exact code you're using right now once more?
That way I can safely assume the same constants are used and compare output.
 
  • #27
EngWiPy
1,367
61
All I'm saying is that they aren't the same. And you gave that expression when asked to reproduce the integral you where trying to evaluate.
This is confusing (at best).

Apparently I made a mistake when copying, it isn't important though. (The constant factor before NIntegrate shouldn't be there according to post #1, it shows up in #19 though)

Am I correct with saying ##GSS=100## and ##yQdB=1## since that would mean there's a problem with " 10 10^(1/10)" in the output.
I cannot reproduce the prefactor either (everything is off by a factor of 10)

Edit;
Can you post the exact code you're using right now once more?
That way I can safely assume the same constants are used and compare output.

I agree, it's confusing. I made a mistake. The code I'm using right now, including your line is

Code:
yp = 10^(5/10);
GSS = 100;
Cons = 10^-7;


For[yQdB = -10, yQdB <= 15, yQdB++;
yQ = 10^(yQdB/10);
A1 = 0.5*((GSS^2)*yQ )/yp*
   NIntegrate[
    1/x^2*Exp[(GSS*(x + GSS*yQ))/(x*yp) - x]*
     ExpIntegralE[1, (GSS*(x + GSS*yQ))/(x*yp)], {x, Cons, Infinity},
    PrecisionGoal -> 5, MaxRecursion -> 12];
Print["A1 = ", 0.5*((GSS^2)*yQ)/yp,
  "\[Times] NIntegrate[  1/x^2\[Times]Exp[ ( ", GSS, "\[Times](x +",
  GSS*yQ, "))/(x\[Times]", yp, ") - x]\[Times] ExpIntegralE[1, (",
  GSS, "\[Times](x + ", GSS*yQ, "))/(x\[Times]", yp, ")]"]
]
 
  • #28
Tom.G
Science Advisor
Gold Member
4,463
3,231
Ref Post #18:
Copying equations works best through reply or quote.

Thanks, my bad. The confusion arose because the "PREVIEW..." button does not show the nicely formatted result; so I thought I'd try a different approach.
Next time I'll know.
 
  • #29
EngWiPy
1,367
61
Any hint?
 
  • #30
JorisL
491
189
I'll try to look at it tommorow, I recall getting "Null" results when I saved the results in a table.
 
  • #31
EngWiPy
1,367
61
OK, thanks. I'm trying to do the mathematical derivations all over again. But I suspect it's just a technical issue with Mathematica.
 
  • #32
Hepth
Gold Member
449
39
Code:
For[yQdB = -10, yQdB <= 15, yQdB++;
 yQ = 10^(yQdB/10);
 INTEGRAND = 
  Simplify[1/x^2*Exp[(GSS*(x + GSS*yQ))/(x*yp) - x]*
    Exp[-(GSS*(x + GSS*yQ))/(x*yp) T]/T];
 A1 = 0.5*((GSS^2)*yQ)/yp*
   NIntegrate[INTEGRAND, {T, 1, \[Infinity]}, {x, Cons, Infinity}, 
    PrecisionGoal -> 5, MaxRecursion -> 12, WorkingPrecision -> 20];
 Print["A1 = ", 0.5*((GSS^2)*yQ)/yp, 
  "\[Times] NIntegrate[  1/x^2\[Times]Exp[ ( ", GSS, "\[Times](x +", 
  GSS*yQ, "))/(x\[Times]", yp, ") - x]\[Times] ExpIntegralE[1, (", 
  GSS, "\[Times](x + ", GSS*yQ, "))/(x\[Times]", yp, ")]"]]

Try that.
Basically its evaluating E1 and the Exponential separately, and one is hitting $MaxNumber and the other $MinNumber.

To fix this I just took the definiton of the E1 function (its in its help file under "more information") and then let it combine the two exponentials before integration.

Now of course you have a 2-dimensional numerical integration, which takes a bit longer, but it seems to get the results. It complains about precision, but the results seem ok. You can increase working precision to try to get that error to go away, but it makes it slower and slower.

Does that help?
 
  • #33
EngWiPy
1,367
61
Thanks. But I get results in term of x. How to get numerical values? (I'm no expert in Mathematica)
 
  • #34
Hepth
Gold Member
449
39
Code:
yp = 10^(5/10);
GSS = 100;
Cons = 10^-7;
results = {};
For[yQdB = -10, yQdB <= 15, yQdB++;
 yQ = 10^(yQdB/10);
 INTEGRAND = 
  Simplify[1/x^2*Exp[(GSS*(x + GSS*yQ))/(x*yp) - x]*
    Exp[-(GSS*(x + GSS*yQ))/(x*yp) T]/T];
 A1 = 0.5*((GSS^2)*yQ)/yp*
   NIntegrate[INTEGRAND, {T, 1, \[Infinity]}, {x, Cons, Infinity}, 
    PrecisionGoal -> 5, MaxRecursion -> 12, WorkingPrecision -> 20];
 results = Join[results, {yQdB, A1}];
 Print[{yQdB, A1}];]

Sorry, I left in your old code showing the equation.

This should print you real results, and at the end of the run you can just look at "results".

Though I'm not certain if this is giving the correct result to be honest, as once yQdB passes 5 it changes to something low. Maybe thats how it should be? I don't know.
 
  • #35
EngWiPy
1,367
61
It's difficult to know if it's correct or not, as A1 is just one value from the final result. The final result is in the form of A1+A2-A3, where all values have similar integrad form approximately. I increased the working precision to 40 but I still get the "error" message!! When I combined all the results, I got some negative values. This shouldn't happen, as the final value must be between 0 and 0.5.
 

Suggested for: Error in Mathematica

  • Last Post
Replies
8
Views
3K
  • Last Post
Replies
16
Views
940
  • Last Post
Replies
3
Views
3K
Replies
3
Views
6K
Replies
2
Views
3K
Replies
6
Views
15K
Replies
3
Views
2K
Replies
3
Views
5K
  • Last Post
Replies
1
Views
5K
  • Last Post
Replies
6
Views
10K
Top