# How to prevent numerical instability during NIntegrate

I'm not entirely new to Mathematica but I use it only very sporadically. I'm foremost a C/C++ programmer by profession.

I recently had a few weeks of spare time and launched Mathematica to do some computation on the GGX BRDF model. I wanted to compute the Fresnel Albedo function for all values of roughnesses and incident angles. The Albedo is a directional-hemispherical integral.

I got pretty much everything running the way I wanted but at very small roughness values, I start getting instabilities in the Albedo computations.

To illustrate, I boiled the code with the instability issue down to this:

(* Utility function for creating a unit 3D vector in the upper \
hemisphere from a polar coordinates with radius 1 *)
PolarToNormal[theta_, phi_] := {Sin[theta] Cos[phi], Sin[theta] Sin[phi], Cos[theta]};

(* Utility function for creating a unit 3D vector in the upper hemisphere from a cos *)
CosToNormal[cos_] := {Sqrt[1 - cos^2], 0, cos};

(* Utility function for computing the cosine of the half vector \
between an incident and an exitant vectors *)
cosM[NdotI_, theta_, phi_] := Normalize[
CosToNormal[NdotI] + PolarToNormal[theta, phi]][[3]];

(* The GGX normal distribution function *)
GGXD[alpha_, NdotI_, theta_, phi_] :=
alpha^2/(Pi (cosM[NdotI, theta, phi]^2 (alpha^2 - 1) + 1)^2);

(* Plot of the exitant hemispherical integral of GGX D for incident cosines from 0 to 1 *)
Plot[NIntegrate[
GGXD[0.00005, cosI, theta, phi] Sin[theta] Cos[theta], \
{phi, -Pi, Pi}, {theta, 0, Pi/2}], {cosI, 0.000000001, 1}]


And this is the plot I get:

How can I debug this?

How can I find which part of the integration process is causing this instability?

Are there any Mathematica parameters such as WorkingPrecision, PrecisionGoal, etc. that I can tweak to fix this?

BTW, Is there any way I can find out what is the current working precision on my installed Mathematica?

Under John Boyd's MAG philosophy (make-a-graph), we should examine the graph of a function when a numerical routine behaves oddly. The support for the integral mainly lies in two spikes:

Block[{cosI = 5/10},
Plot3D[GGXD[0.00005, cosI, theta, phi] Sin[theta] Cos[theta],
{phi, -Pi, Pi}, {theta, 0, Pi/2}, PlotPoints -> 25,
PlotRange -> All, PlotLabel -> N@ArcCos[cosI]]
]


These spikes might be missed by the default sampling. We ca add points to the iterators to ensure that they will be found. The occur where the denominator of GGXD[5/100000, cosI, theta, phi] has a minimum. (You can get the denominator with GGXD[5/100000, cosI, theta, phi] /. Abs -> RealAbs // Simplify // Denominator.)

Solve[
D[(1 - (399999999 (cosI + Cos[theta])^2)/(800000000 (1 + cosI Cos[theta] +
Sqrt[1 - cosI^2] Cos[phi] Sin[theta])))^2 /. phi -> Pi,
theta] == 0, theta]


Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

(*
{{theta -> -ArcCos[-cosI]}, {theta -> ArcCos[cosI]},
{theta -> -ArcCos[(400000001 cosI - 40000 Sqrt[-1 + cosI^2])/399999999]},
{theta -> ArcCos[(400000001 cosI - 40000 Sqrt[-1 + cosI^2])/399999999]},
{theta -> -ArcCos[(400000001 cosI + 40000 Sqrt[-1 + cosI^2])/399999999]},
{theta -> ArcCos[(400000001 cosI + 40000 Sqrt[-1 + cosI^2])/399999999]}}
*)


The second solution, theta == ArcCos[cosI], is the location of the spike near phi == {-Pi, Pi}. The width in the phi direction of the spikes is quite small, and adding two points near phi == {-Pi, Pi} ensures that the spike will be adequately sampled by NIntegrate.

Plot[NIntegrate[
GGXD[0.00005, cosI, theta, phi] Sin[theta] Cos[
theta], {phi, -Pi, -Pi + 10^-3, Pi - 10^-3, Pi}, {theta, 0,
ArcCos[cosI], Pi/2}], {cosI, 0.000000001, 1}, MaxRecursion -> 2]


• Thanks Michael. I see what you are doing. Very good debugging technique and I reused it on the same function but on the cos theta dimension. I have a question: The change you suggest will add more samples by overlapping the troublesome region. This will change the probability for the samples in this region. I guess Mathematica will detect the overlap and adjust the samples probabilities accordingly? Is that correct? – Yves Poissant Jun 17 '18 at 23:24
• @YvesPoissant Yes NIntegrate computes the integral by adjusting the rule to the density of the sampling in each subregion. (In fact, it almost certainly subdivides further some regions automatically; the process is called "adaptive sampling." It will increase sampling in regions in order to reduce the error in the parts that have the worst error.) – Michael E2 Jun 19 '18 at 2:01