# NIntegrate gives the wrong result

I have tried to calculate an integral of a function numerically using NIntegrate, but sadly it seems like I do not get the correct result from it.

The function $$\phi(x,y)$$ is defined using a function $$A(c)$$ which is given by an parametric numerical integral. The definition is

A[c_] := Evaluate[2 NIntegrate[(2 ArcCos[Sqrt[1 - c]/Sqrt[(c + Cos[π x]) Sec[(πx)/2]^2]])/π,
{x, -((2 ArcCos[Sqrt[1 - c]/Sqrt[1 + c]])/π), (2 ArcCos[Sqrt[1 - c]/Sqrt[1 + c]])/π}]]

r[x_, y_] := 1/Sqrt[π] Sqrt[A[1 - 2/(Cos[π/2 x]^-2 + Cos[π/2 y]^-2)]]

a = Sqrt[π]/2;
ϕ[x_, y_] := 2 ArcSin[1/Sqrt[2] a r[x/a, y/a]];


The function $$\phi(x,y)$$ defined in the triangle with corner points of $$(0,0), (a,a), (a,-a)$$ for $$a=\frac{\pi}{2}$$. By ploting the function $$\phi(x,y)$$ we can see it is pretty simillar to the function $$\frac{\pi}{2a}x$$:

Plot3D[{ϕ[x, y], π/(2 a) x}, {x, 0, a}, {y, -x, x}]


Therefore I expect the integral $$\intop_{0}^{a}\intop_{-x}^{x}\left(\frac{\partial\phi}{\partial x}\right)^{2}dydx$$ to be close to $$a^2\cdot(\frac{\pi}{2a})^2=\frac{\pi^2}{4}\approx2.467$$. But when I'm integrating numerically over $$\left(\partial_{x}\phi\right)^{2}$$ I get something very different:

In[69]:= NIntegrate[(\!$$\*SubscriptBox[\(∂$$, $$x$$]$$ϕ[x, y]$$\))^2, {x, 0, a}, {y, -x, x}]

During evaluation of In[69]:= NIntegrate::nlim: x = Indeterminate is not a valid limit of integration.

During evaluation of In[69]:= General::ivar: Sqrt[π]/2 is not a valid variable.

During evaluation of In[69]:= NIntegrate::nlim: x = Indeterminate is not a valid limit of integration.

During evaluation of In[69]:= General::ivar: Sqrt[π]/2 is not a valid variable.

During evaluation of In[69]:= NIntegrate::nlim: x = Indeterminate is not a valid limit of integration.

During evaluation of In[69]:= General::stop: Further output of NIntegrate::nlim will be suppressed during this calculation.

During evaluation of In[69]:= General::ivar: Sqrt[π]/2 is not a valid variable.

During evaluation of In[69]:= General::stop: Further output of General::ivar will be suppressed during this calculation.

During evaluation of In[69]:= NIntegrate::write: Tag Times in (Sqrt[π] x)/2 is Protected.

During evaluation of In[69]:= NIntegrate::write: Tag Times in (Sqrt[π] x)/2 is Protected.

During evaluation of In[69]:= NIntegrate::write: Tag Times in (Sqrt[π] x)/2 is Protected.

During evaluation of In[69]:= General::stop: Further output of NIntegrate::write will be suppressed during this calculation.

During evaluation of In[69]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

During evaluation of In[69]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 18 recursive bisections in x near {x,y} = {0.93974199467253037576508489792104228399693965911865234375000000000,0.713879}. NIntegrate obtained 0.605531 +0.0100204 I and 0.00006594026997516764 for the integral and error estimates.

Out[69]= 0.605531 + 0.0100204 I


I get many warnings which are probably due to the definition of $A(c)$, but in other cases I did this I got the correct results despite these warnings. In this case I get that the integral fail to converge, and indeed the result I get here is $0.605531$ (plus an imagenary contribution) is probablt too far from the approximately expected $2.467$. I tried many integration methods using the Methods option, but non of them helped. Monte-Carlo methods also gave results of about $$0.605$$.

In addition I tried to calculate the integral by definin $$A(c)$$ using A[c_?NumericQ]. In this case the integral did not give results at all.

Does anyone know what the problem is? Am I wrong and the integral is so different from the one over $$\frac{\pi}{2a}$$ despite the similar apperence of the functions?

Thank you!

• It looks like there's a \[Pi]x that should be \[Pi]*x, but that might be a copy-and-paste problem. That's at least one possible coding issue, here. Also, I don't think you you can take the derivative of a numerical function that's defined like A[x_?NumericQ] := ... Commented Dec 16, 2022 at 21:55
• @march Thank you for mentioning, the $\pi$ is indeed just a copy-and-paste problem. It did not give a problem in the code. About the ?NumericQ, you are right. I just tried both with and without the ?NumericQ to see if it changes the result, but it doesn't.
– Roy
Commented Dec 16, 2022 at 22:06
• I am trying to run your code as we speak, but it's running for a long time. It hasn't kicked back any errors yet, though. In any case, the point is that A is defined in terms of a numerical integral, and I don't think that D is equipped to handle that situation. Commented Dec 16, 2022 at 22:09
• Thank you very much @march. It seems like I messed up, and another thing I wrote in the notebook interacted with the calculation for reason I'm not sure of. The results I wrote are actually the results I get if I do not use ?NumericQ. If I do, the integral just does not give any results, as you mentioned. I updated the question so it would be clear.
– Roy
Commented Dec 16, 2022 at 23:19

Therefore I expect the integral $$\intop_{0}^{a}\intop_{-x}^{x}\left(\frac{\partial\phi}{\partial x}\right)^{2}dydx$$ to be close to $$a^2\cdot(\frac{\pi}{2a})^2=\frac{\pi^2}{4}\approx2.467$$

If functions $$f(x)$$ and $$g(x)$$ are close on some interval, the derivatives $$f'(x)$$ and $$g'(x)$$ are not necessary close there. In fact, in your case the code

A[c_] := Evaluate[ 2 NIntegrate[(2 ArcCos[
Sqrt[1 - c]/
Sqrt[(c +  Cos[\[Pi] x]) Sec[(\[Pi] x)/
2]^2]])/\[Pi], {x, -((2 ArcCos[
Sqrt[1 - c]/Sqrt[1 + c]])/\[Pi]), (2 ArcCos[
Sqrt[1 - c]/Sqrt[1 + c]])/\[Pi]}]]

r[x_, y_] := 1/Sqrt[\[Pi]] Sqrt[A[1 - 2/(Cos[\[Pi]/2 x]^-2 + Cos[\[Pi]/2 y]^-2)]]

a = Sqrt[\[Pi]]/2;\[Phi][x_, y_] := 2 ArcSin[1/Sqrt[2] a r[x/a, y/a]];

Plot3D[{Evaluate[D[\[Phi][x, y], x]], \[Pi]/(2 a) }, {x, 0, a}, {y, -x, x}]


results in

As we see, the plots are not close.

• Thank you @user64494, but I do not see how it makes sense. We can see from the plot of $\phi(x,y)$ that its derivative is around to the one of $\frac{\pi}{2a}x$. If we try ploting constant $y$ slices of $\phi(x,y)$, and then we can clearly see that its derivative is sometimes above and sometimes below the one $\frac{\pi}{2a}x$, but the function you ploted gives only smaller derivative. This is probably a numerical error. Following @Bob Hanlon's, by taking the derivative using ND and ploting we get a different function, that has values around $\frac{\pi}{2a}$.
– Roy
Commented Dec 17, 2022 at 10:50
• @Roy: I am out of MMA now so I'll reply later. Commented Dec 17, 2022 at 14:32
• @Roy: (i) In the Bob Hanlon's notations Plot3D[D[\[Phi][x, y], x], {x, 0, a}, {y, -x, x}] produces an empty plot. (ii) The plots in ny answer and the Bob Hanlon's numeric result 2.767 are in accordance. The derivatives differ by $\approx 0.5$. This implies the integrals differ by $\approx 0.5\cdot 1\cdot 1 /2=0.25$. Commented Dec 17, 2022 at 19:00
• But as Bob Hanlon specified, using ND instead of D does give a plot. Also taking the numerical derivative by finite differences gives a plot similar to the one ND gives. If the derivatives differ by $\approx0.5$, than the why would the integrals differ by about $0.25$? Did you mean to multiply by the area of the triangle? The area is $1/2\cdot a \cdot 2a=a^2$, and the integrals should differ by $~0.5a^2\approx0.4$. In addition, in the plot you showed the derivative is smaller than $\pi/2a$, but Bob Hanlon's result is larger than the one given from $\pi/2a$.
– Roy
Commented Dec 17, 2022 at 23:44
Clear["Global*"]


The argument of a function which uses a numeric technique should be restricted to numeric values. However, do not put Evaluate on the RHS of the definition. The Evaluate will try to evaluate immediately before there is a numeric value for the argument.

A[cv_?NumericQ] := Module[
{c = SetPrecision[cv, Max[15, Precision[cv]]]},
2 NIntegrate[(2 ArcCos[Sqrt[1 - c]/
Sqrt[(c + Cos[π x]) Sec[(π x)/2]^2]])/π,
{x, -((2 ArcCos[Sqrt[1 - c]/Sqrt[1 + c]])/π),
(2 ArcCos[Sqrt[1 - c]/Sqrt[1 + c]])/π},
WorkingPrecision -> Max[15, Precision[cv]]]]

r[x_, y_] := 1/Sqrt[π] *
Sqrt[A[1 - 2/(Cos[π/2 x]^-2 + Cos[π/2 y]^-2)]]

a = Sqrt[π]/2;

ϕ[x_, y_] := 2 ArcSin[1/Sqrt[2] a r[x/a, y/a]];


To take the derivative of a numeric function you will need to use a numeric technique, i.e., ND

Needs["NumericalCalculus"]

NIntegrate[
(ND[ϕ[x, y], x, x0, WorkingPrecision -> 12])^2,
{x0, 0, a}, {y, -x0, x0},
AccuracyGoal -> 10,
WorkingPrecision -> 12]

(* 2.76682450023 *)

• Thanks a lot @Bob Hanlon! This seems to be the solution. I still wonder what exactly Mathematica is doing when I'm asking for a derivative using "regular" derivative commands that gives such wrong results.
– Roy
Commented Dec 17, 2022 at 11:17
• I checked it more, and something is still different from what I expect. I tried to calculate the numeriacal derivative by myself using finite differences by d\[Phi][x_, y_] := (\[Phi][x + h, y] - \[Phi][x - h, y])/(2 h) for small $h$ (h = 0.00000001 for example). Integrating over the square of this function gives a different result of 2.3798, which is about 0.4 less. Changing the AccuracyGoal and WorkingPrecision does not change this result significantly. The ND` function is probably more accurate, but with such differences I'm not sure I can trust the results I get...
– Roy
Commented Dec 17, 2022 at 13:03