How to get the Linf norm of the solution

[zhyxue0]

In order to satisfy a CFL-type stability condition, I can get the maximum value of u by

np.abs(u.vector().array() ).max()

But this is just the maximum of the polynomial coefficients, which is very different from the pointwise maximum.

In order to get the Linf_error

Linf_error = np.abs( u.vector().array() - u_exact.vector().array() ).max()

Same problem here, but worse getting the wrong order.

norms.errornorm(u, u_exact, 'Lp' ) can only get the error of Lp (p>0), so how to get the Linf norm of the solution ?

[wence-]

tl;dr: it's not as trivial as you might think.

You can get an estimate of the Linf norm if you can put your function in a bernstein space (fiat at least implements bernstein elements, but I don't think they're hooked up). These basis functions are bezier curves which have the nice property that their maximum is bounded by the convex hull of the control points. So once you have the function in that space, pointwise maxing over the dofs gives you something that approximates the Linf norm.

[dham]

I believe Bernstein is wired in for intervals, triangles and tets. You can probably tensor product elements for quads, wedges and hexes too.

[wence-]

I guess the other thing you can do, if you know the basis functions, is just solve for the max on each element. But I don't think we have much facility to help with that.

[colinjcotter]

I think probably we are making too much fuss about this. There must be results in the literature for convergence of infinity norm evaluated over quadrature points to the true infinity norm for sufficiently smooth functions, so it might be reasonable to provide this. I just can't find these results because searching for L-infinity and approximation just provides a deluge of papers about L-infinity estimates for approximation error, not approximation of the L-infinity norm itself.

[rckirby]

Depending on the polynomial basis, the max of the polynomial may be many times the max of the nodal values . The Gauss-Lovatto points have a Lebesgue constant that’s not too big… One approach would be to use GLL and multiplayer the max value by the Lebesgue constant to be safe. This doesn’t help finding the max-norm error, just the max-norm of a polynomial.

[zhyxue0]

Maybe we can take a certain maximum value in the numerical sense. For example, for 1D divided into N cells, we can take the values of 10 equal points( or more) in each interval, and then get 10*N point values.

import numpy as np
def Max_Pts_Func(mesh_coord):
    x = np.unique(mesh_coord)
    ref_pts = np.linspace(0,1, 10)
    pts = np.zeros( (len(x)-1, len(ref_pts) ) )

    for i in range(len(ref_pts)):
        pts[:,i] = 0.5*(x[1:] + x[0:-1]) + 0.5*ref_pts[i]*(x[1:] - x[0:-1])
    pts = pts.flatten()

    return pts

...
...
mesh = IntervalMesh(Nx, Lx)
mesh_coord = mesh.coordinates.dat.data_ro
...
...
pts = Max_Pts_Func(mesh_coord)
u_max = np.isfinite(u.at(pts, dont_raise=True)).max()
Linf_error= np.abs( np.isfinite(u.at(pts, dont_raise=True) - u_exact.at(pts, dont_raise=True))).max()

But this work efficiency is low.
Number of points to test for 2D or 3D:
2D: (10*N)^2

3D: (10*N)^3