# Need help to speed up the Gaussian quadrature

I am trying to compute the $$L_2$$-norm of the solution $$y(z,t)$$ of a PDE with Gauss quadrature using the following code, where $$z$$ is space position and $$t$$ is time, then to construct it as a function of $$t$$ by interpolating.

Get["NumericalDifferentialEquationAnalysis"];

np = 101; points = weights = Table[Null, {np}];
Do[points[[i]] = GaussianQuadratureWeights[np, 0, L][[i, 1]], {i, 1, np}]
Do[weights[[i]] = GaussianQuadratureWeights[np, 0, L][[i, 2]], {i, 1, np}]
GaussInt[f_(*integrand*), z_] := Sum[(f /. z -> points[[i]])*weights[[i]], {i, 1, np}]
GQL2norm = Table[{t, GaussInt[sol[z, t]^2, z]}, {t, 0, tm, 1}];
L2normInt = Interpolation[GQL2norm, InterpolationOrder -> 2, Method -> "Spline"];


Here, sol[z,t] is an InterpolatingFunction that was obtained from NDSolveValue.

The following PDE is for testing, although the computation using the above code for the specific example is indeed fast, it is really slow for my real problem.

tm = 10; L = 5;
sol = NDSolveValue[{\!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]$$u[z, t]$$\) == \!$$\*SubscriptBox[\(\[PartialD]$$, $$z, z$$]$$u[z, t]$$\), u[z, 0] == 0,
u[0, t] == Sin[t], u[5, t] == 0}, u, {t, 0, tm}, {z, 0, L}]


Plot the evolution of the $$L_2$$-norm:

Plot[L2normInt[t]^(1/2), {t, 0, tm}, PlotRange -> {{0, tm}, All}, Frame -> True]


Problem:

1. The code can give a correct result with a sufficient node points np, but it runs for an extremely long time. Please help me to improve it.

2. As we do not know the degree of polynomials for an InterpolatingFunction obtained from NDSolve-type solver, is there a rule of thumb to estimate the an adequate value of np?

From wikipedia:

An $$n$$-point Gaussian quadrature rule is a quadrature rule constructed to yield an exact result for polynomials of degree $$2n − 1$$ or less by a suitable choice of the nodes $$x_i$$ and weights $$w_i$$ for $$i = 1, ..., n.$$

Thank you for any suggestion!

• Your sample is incomplete. GaussianQuadratureWeights is from NumericalDifferentialEquationAnalysis, you should add the corresponding Needs[……] to the sample. Also, a sol should be added for testing. Commented Aug 6, 2022 at 3:47
• As to Gauss-Legendre quadrature, an efficient implemention can be found here: mathematica.stackexchange.com/a/6966/1871 Commented Aug 6, 2022 at 3:53
• Integrating a piecewise polynomial sol[z, t]^2 should be a simple business normally....I cannot say what would be the particular problem with your sol[t, z], though. Commented Aug 6, 2022 at 4:27
• @xzczd, thank you, I updated. The practical sol is too large to add here. Let me find a simple example Commented Aug 6, 2022 at 4:53
• @xzczd, I added an example for testing. Could you adapt the implement of Gauss quadrature in that link to fit an InterpolatingFunction obtained from NDSolveValue. Commented Aug 6, 2022 at 5:52

First of all, you should really put a bit more effort in learning the core language of Mathematica. (Look at those annoying loops! ) Then, if you insist on using Gauss quadrature and speed is concerned, you should use Gauss quadrature nodes from the very beginning i.e. use it as spatial grid of NDSolve:

AbsoluteTiming[
points = 500; domain = {0, L};
{nodes, weights} = Most[NIntegrateGaussRuleData[points, MachinePrecision]];
midgrid = Rescale[nodes, {0, 1}, domain];
grid = Flatten[{domain[[1]], midgrid, domain[[-1]]}];

tm = 10; L = 5;
solmodified = NDSolveValue[{D[u[z, t], t] == D[u[z, t], z, z], u[z, 0] == 0,
u[0, t] == Sin[t], u[5, t] == 0}, u, {t, 0, tm}, {z, grid} // Flatten];

solmat = solmodified["ValuesOnGrid"];
timegrid = solmodified["Coordinates"][[2]];

int = -Subtract @@ domain weights . solmat[[2 ;; -2]]^2;

mynormint = ListInterpolation[int, timegrid]]


Though I've chosen an unnecessarily dense grid (points = 500), the timing is just about 0.1 second on my laptop, while yours takes about 3.5 second.

To understand my code, consider reading

If you don't bother to re-evaluate NDSolve and Gauss quadrature isn't necessary for you, then you should turn to trapezoid rule, because trapezoid rule doesn't require one to use a special grid, the default uniform grid of NDSolve is enough. Then again we can directly extract the solution list and integrate:

AbsoluteTiming[
points = sol["Coordinates"][[1]] // Length; domain = {0, L};
With[{h = -Subtract @@ domain/(points - 1)},
trap[value_] := h (Total[value] - 1/2 (value[[1]] + value[[-1]]))];
solmat = sol["ValuesOnGrid"];
timegrid = sol["Coordinates"][[2]];
int = solmat^2 // trap;

mynorminttrap = ListInterpolation[int, timegrid]]
(* Timing: 0.0014401 *)

ListPlot[{L2normInt, mynormint, mynorminttrap}, PlotMarkers -> "OpenMarkers"]


• thank you for the optimal code and useful links. Well, my concern is how to speed up the Gauss quadrature itself after one has obtained a solution from NDSolveValue, which did not use Gauss nodes as spatial grid. Also, the computation of NDSolve is slow. Could you just apply J.M.'s GaussLegendreQuadrature function to the sol from NDSolveValue directly? Sorry for the further request! Commented Aug 6, 2022 at 8:04
• @Nobody Then just use the GaussLegendreQuadrature as shown in J.M.'s post. Its usage is quite straightforward, which part are you having difficulty? Also, since the default grid of NDSolve is uniform grid, turning to e.g. trapezoid rule instead of Gauss quadrature and use my method above should be quite efficient. I've updated the answer to elaborate. Commented Aug 6, 2022 at 8:46
• @Nobody, if you're using a low interpolation order to begin with, then the extra effort and high order in Gaussian quadrature is wasted, and as xzczd says, trapezoidal might just be good enough for your purposes. Commented Aug 6, 2022 at 8:51

Here's an analytic solution. You have to know the order of the interpolation, which depends on the (spatial) order of the pde if using InterpolationOrder -> Automatic in NDSolve.

intOrder2 = (CoefficientList[
InterpolatingPolynomial[{{{z1}, y1, p1, pp1}, {{z2}, y2, p2,
pp2}}, z]^2, z]/Range[11]) . (z^Range[11]) /. {{z ->
z2}, {z -> z1}} // Apply@Subtract // Simplify
(* can use Integrate[..., {z, z1, z2}] instead
-(1/27720)(10860 y1^2 + 10860 y2^2 +
y1 (6000 y2 + (-3732 p1 +
1812 p2 + (281 pp1 + 181 pp2) (z1 - z2)) (z1 - z2)) +
y2 (-1812 p1 + 3732 p2 + (181 pp1 + 281 pp2) (z1 - z2)) (z1 -
z2) + (416 p1^2 - 532 p1 p2 + 416 p2^2 -
p1 (69 pp1 + 52 pp2) (z1 - z2) +
p2 (52 pp1 + 69 pp2) (z1 - z2) + (3 pp1^2 + 5 pp1 pp2 +
3 pp2^2) (z1 - z2)^2) (z1 - z2)^2) (z1 - z2)
*)

zvals = Transpose@
ConstantArray[sol["Coordinates"][[1]],
Length[sol["Coordinates"][[2]]]];
yvals = Derivative[0, 0][sol]["ValuesOnGrid"];
pvals = Derivative[1, 0][sol]["ValuesOnGrid"];
ppvals = Derivative[2, 0][sol]["ValuesOnGrid"];
normIFN = intOrder2 /.
Thread[{z1, y1, p1, pp1, z2, y2, p2, pp2} ->
Flatten[{Most /@ {zvals, yvals, pvals, ppvals},
Rest /@ {zvals, yvals, pvals, ppvals}}, 1]
] //
Transpose //
Map@Total //
Transpose[{sol["Coordinates"][[2]], #}] & //
Interpolation

ListPlot[normIFN]
`