Need help to speed up the Gaussian quadrature

Question

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!

Answer 1

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[NIntegrate`GaussRuleData[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

https://mathematica.stackexchange.com/a/6966/1871 https://mathematica.stackexchange.com/a/240406/1871 https://mathematica.stackexchange.com/a/19043/1871 https://mathematica.stackexchange.com/a/28341/1871 https://mathematica.stackexchange.com/a/163273/1871

---

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"]

Answer 2

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]