数値実験編：外挿，defect correction，特異性の比較

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Extrapolation, defect correction, and singularity experiments.

This script accompanies the blog article:
外挿法と誤差展開（第8回）：数値実験編
外挿，defect correction，特異性の比較

The script has three parts.

Part A.
  A small 1D periodic model problem for defect correction.
  We solve
      -u'' + u = f
  on [0,1] with periodic boundary condition.
  The exact solution is u(x)=sin(2*pi*x).
  A second-order finite difference solution is improved by one defect
  correction step using a fourth-order residual.

Part B.
  A 2D P1 finite element code for the Poisson equation.
  We compare pointwise Richardson extrapolation on:
    - a smooth square-domain problem,
    - an L-shaped-domain singular problem.

Part C.
  A simple residual/defect indicator for the 2D P1 FEM solution.
  This is not yet the full Rannacher-type FEM defect correction.
  It is included to visualise where the defect concentrates,
  especially for the L-shaped-domain singular problem.

Requirements:
  numpy, scipy, matplotlib
"""

import math
import numpy as np
from scipy.sparse import csr_matrix, lil_matrix, diags
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import matplotlib.tri as mtri


# ============================================================
# Part A. 1D defect correction toy model
# ============================================================

def periodic_matrix_second_order(N):
    """
    Matrix for -u'' + u on a periodic grid with N points.
    The second-order approximation is used for -u''.
    """
    h = 1.0 / N

    main = (2.0 / h**2 + 1.0) * np.ones(N)
    off = (-1.0 / h**2) * np.ones(N - 1)

    A = diags([main, off, off], [0, -1, 1], shape=(N, N), format="lil")
    A[0, N - 1] = -1.0 / h**2
    A[N - 1, 0] = -1.0 / h**2

    return A.tocsr()


def periodic_matrix_fourth_order(N):
    """
    Matrix for -u'' + u on a periodic grid with N points.
    The fourth-order approximation is used for -u''.

    Formula:
        -u''(x_i) ≈
        (30 u_i -16(u_{i-1}+u_{i+1}) + (u_{i-2}+u_{i+2}))/(12 h^2).
    """
    h = 1.0 / N

    main = (30.0 / (12.0 * h**2) + 1.0) * np.ones(N)
    off1 = (-16.0 / (12.0 * h**2)) * np.ones(N - 1)
    off2 = (1.0 / (12.0 * h**2)) * np.ones(N - 2)

    A = diags(
        [main, off1, off1, off2, off2],
        [0, -1, 1, -2, 2],
        shape=(N, N),
        format="lil",
    )

    # periodic wrap for distance 1
    A[0, N - 1] = -16.0 / (12.0 * h**2)
    A[N - 1, 0] = -16.0 / (12.0 * h**2)

    # periodic wrap for distance 2
    A[0, N - 2] = 1.0 / (12.0 * h**2)
    A[1, N - 1] = 1.0 / (12.0 * h**2)
    A[N - 2, 0] = 1.0 / (12.0 * h**2)
    A[N - 1, 1] = 1.0 / (12.0 * h**2)

    return A.tocsr()


def defect_correction_1d(N):
    """
    One defect correction step for the periodic model problem.

    Exact solution:
        u(x)=sin(2*pi*x).

    Equation:
        -u'' + u = f.
    """
    x = np.arange(N) / N
    exact = np.sin(2.0 * math.pi * x)
    f = (1.0 + (2.0 * math.pi) ** 2) * exact

    A2 = periodic_matrix_second_order(N)
    A4 = periodic_matrix_fourth_order(N)

    # Low-order solution.
    u0 = spsolve(A2, f)

    # Defect measured by the high-order operator.
    defect = A4 @ u0 - f

    # Correct using the stable low-order operator.
    correction = spsolve(A2, -defect)
    u1 = u0 + correction

    err0 = np.linalg.norm(u0 - exact) / math.sqrt(N)
    err1 = np.linalg.norm(u1 - exact) / math.sqrt(N)

    return err0, err1


def run_defect_correction_1d():
    print("=" * 80)
    print("Part A. 1D defect correction toy model")
    print("Problem: -u'' + u = f, periodic, exact u=sin(2*pi*x)")
    print("=" * 80)
    print(f"{'N':>6} | {'err low-order':>14} {'rate':>7} | {'err corrected':>14} {'rate':>7}")
    print("-" * 80)

    old0 = None
    old1 = None

    for N in [16, 32, 64, 128, 256]:
        err0, err1 = defect_correction_1d(N)

        if old0 is None:
            r0 = "   ---"
            r1 = "   ---"
        else:
            r0 = f"{math.log(old0 / err0, 2.0):7.3f}"
            r1 = f"{math.log(old1 / err1, 2.0):7.3f}"

        print(f"{N:6d} | {err0:14.6e} {r0} | {err1:14.6e} {r1}")

        old0 = err0
        old1 = err1

    print()
    print("Expected behaviour:")
    print("  low-order solution       : about second order")
    print("  defect-corrected solution: about fourth order")
    print()


# ============================================================
# Part B. 2D P1 FEM for the Poisson equation
# ============================================================

def generate_square_mesh(N):
    """
    Uniform triangular mesh of [0,1]^2.

    N is the number of subdivisions in each coordinate direction.
    """
    nodes = []
    index = {}

    for j in range(N + 1):
        for i in range(N + 1):
            x = i / N
            y = j / N
            index[(i, j)] = len(nodes)
            nodes.append((x, y))

    triangles = []

    for j in range(N):
        for i in range(N):
            v00 = index[(i, j)]
            v10 = index[(i + 1, j)]
            v01 = index[(i, j + 1)]
            v11 = index[(i + 1, j + 1)]

            # Split each square along the same diagonal.
            triangles.append((v00, v10, v11))
            triangles.append((v00, v11, v01))

    return np.array(nodes, dtype=float), np.array(triangles, dtype=int)


def generate_lshape_mesh(N):
    """
    Uniform triangular mesh of the L-shaped domain

        Omega = (-1,1)^2 \\ ([0,1] x [-1,0]).

    The reentrant corner is at the origin.
    N should be even so that the origin and common test points are grid points.
    """
    nodes = []
    index = {}

    for j in range(N + 1):
        y = -1.0 + 2.0 * j / N
        for i in range(N + 1):
            x = -1.0 + 2.0 * i / N
            index[(i, j)] = len(nodes)
            nodes.append((x, y))

    triangles = []

    for j in range(N):
        for i in range(N):
            # cell centre
            xc = -1.0 + 2.0 * (i + 0.5) / N
            yc = -1.0 + 2.0 * (j + 0.5) / N

            # remove the fourth quadrant
            if (xc > 0.0) and (yc < 0.0):
                continue

            v00 = index[(i, j)]
            v10 = index[(i + 1, j)]
            v01 = index[(i, j + 1)]
            v11 = index[(i + 1, j + 1)]

            triangles.append((v00, v10, v11))
            triangles.append((v00, v11, v01))

    # Remove unused nodes.
    used = sorted(set(int(v) for tri in triangles for v in tri))
    new_id = {old: new for new, old in enumerate(used)}
    new_nodes = np.array([nodes[i] for i in used], dtype=float)
    new_tris = np.array([[new_id[int(v)] for v in tri] for tri in triangles], dtype=int)

    return new_nodes, new_tris


def boundary_nodes_from_triangles(triangles):
    """
    Return boundary nodes by detecting edges that belong to only one triangle.
    """
    edge_count = {}

    for tri in triangles:
        for a, b in [(tri[0], tri[1]), (tri[1], tri[2]), (tri[2], tri[0])]:
            edge = tuple(sorted((int(a), int(b))))
            edge_count[edge] = edge_count.get(edge, 0) + 1

    bnodes = set()
    bedges = []

    for edge, count in edge_count.items():
        if count == 1:
            bedges.append(edge)
            bnodes.update(edge)

    return np.array(sorted(bnodes), dtype=int), bedges


def assemble_poisson(nodes, triangles, f):
    """
    Assemble the P1 stiffness matrix and load vector for

        -Delta u = f.

    The load vector uses centroid quadrature.
    """
    nnode = nodes.shape[0]
    rows = []
    cols = []
    data = []
    rhs = np.zeros(nnode)

    for tri in triangles:
        coords = nodes[tri, :]
        x0, y0 = coords[0]
        x1, y1 = coords[1]
        x2, y2 = coords[2]

        detJ = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)
        area = abs(detJ) / 2.0

        # Gradients of barycentric basis functions.
        b = np.array([y1 - y2, y2 - y0, y0 - y1]) / detJ
        c = np.array([x2 - x1, x0 - x2, x1 - x0]) / detJ
        grads = np.vstack([b, c]).T

        localK = area * (grads @ grads.T)

        xc = (x0 + x1 + x2) / 3.0
        yc = (y0 + y1 + y2) / 3.0
        localF = area * f(xc, yc) / 3.0 * np.ones(3)

        for a in range(3):
            A = int(tri[a])
            rhs[A] += localF[a]
            for bidx in range(3):
                B = int(tri[bidx])
                rows.append(A)
                cols.append(B)
                data.append(localK[a, bidx])

    K = csr_matrix((data, (rows, cols)), shape=(nnode, nnode))
    return K, rhs


def solve_dirichlet_poisson(nodes, triangles, f, g):
    """
    Solve the P1 FEM system with Dirichlet boundary value g.
    """
    K, rhs = assemble_poisson(nodes, triangles, f)
    bnodes, _ = boundary_nodes_from_triangles(triangles)

    all_nodes = np.arange(nodes.shape[0])
    is_boundary = np.zeros(nodes.shape[0], dtype=bool)
    is_boundary[bnodes] = True

    interior = all_nodes[~is_boundary]
    boundary = all_nodes[is_boundary]

    u = np.zeros(nodes.shape[0])
    u[boundary] = np.array([g(nodes[i, 0], nodes[i, 1]) for i in boundary])

    rhs_i = rhs[interior] - K[interior, :][:, boundary] @ u[boundary]
    K_ii = K[interior, :][:, interior]

    u[interior] = spsolve(K_ii, rhs_i)

    return u


def find_node(nodes, point, tol=1.0e-12):
    """
    Find a node with coordinate equal to point.
    """
    p = np.array(point, dtype=float)
    dist = np.linalg.norm(nodes - p[None, :], axis=1)
    idx = int(np.argmin(dist))

    if dist[idx] > tol:
        raise ValueError(f"Point {point} is not a mesh node. nearest distance={dist[idx]}")

    return idx


def smooth_exact(x, y):
    return math.sin(math.pi * x) * math.sin(math.pi * y)


def smooth_rhs(x, y):
    return 2.0 * math.pi**2 * smooth_exact(x, y)


def lshape_angle(x, y):
    """
    Angle in the L-shaped domain with the removed fourth quadrant.

    Near the reentrant corner, the domain corresponds to
        0 < theta < 3*pi/2.
    """
    theta = math.atan2(y, x)
    if theta < 0.0:
        theta += 2.0 * math.pi
    return theta


def singular_exact_factory(alpha):
    def exact(x, y):
        r = math.hypot(x, y)
        if r == 0.0:
            return 0.0
        theta = lshape_angle(x, y)
        return (r ** alpha) * math.sin(alpha * theta)
    return exact


def zero_rhs(x, y):
    return 0.0


def run_point_extrapolation_smooth():
    """
    Smooth square-domain test.

    exact u = sin(pi x) sin(pi y).
    We evaluate at P=(1/2,1/2).
    """
    print("=" * 80)
    print("Part B1. Smooth square-domain test")
    print("u=sin(pi x) sin(pi y), P=(1/2,1/2)")
    print("=" * 80)

    point = (0.5, 0.5)
    exact_value = smooth_exact(*point)

    values = {}

    for N in [8, 16, 32, 64]:
        nodes, tris = generate_square_mesh(N)
        uh = solve_dirichlet_poisson(nodes, tris, smooth_rhs, smooth_exact)
        idx = find_node(nodes, point)
        values[N] = uh[idx]

    print(f"{'N':>6} | {'P1 value':>14} {'error':>12} | {'Richardson':>14} {'error':>12}")
    print("-" * 80)

    for N in [8, 16, 32]:
        val_h = values[N]
        val_h2 = values[2 * N]
        ext = (4.0 * val_h2 - val_h) / 3.0

        err = abs(val_h - exact_value)
        err_ext = abs(ext - exact_value)

        print(f"{N:6d} | {val_h:14.10f} {err:12.4e} | {ext:14.10f} {err_ext:12.4e}")

    print()


def run_point_extrapolation_lshape():
    """
    L-shaped-domain singular test.

    exact u = r^alpha sin(alpha theta), alpha=2/3.
    We evaluate at P=(-1/2,1/2).
    """
    print("=" * 80)
    print("Part B2. L-shaped-domain singular test")
    print("u=r^alpha sin(alpha theta), alpha=2/3, P=(-1/2,1/2)")
    print("=" * 80)

    alpha = 2.0 / 3.0
    exact = singular_exact_factory(alpha)
    point = (-0.5, 0.5)
    exact_value = exact(*point)

    values = {}

    for N in [16, 32, 64, 128]:
        nodes, tris = generate_lshape_mesh(N)
        uh = solve_dirichlet_poisson(nodes, tris, zero_rhs, exact)
        idx = find_node(nodes, point)
        values[N] = uh[idx]

    q = 2.0 * alpha
    factor = 2.0 ** q

    print(f"{'N':>6} | {'P1 value':>14} {'error':>12} | {'std ext':>14} {'error':>12} | {'singular ext':>14} {'error':>12}")
    print("-" * 112)

    for N in [16, 32, 64]:
        val_h = values[N]
        val_h2 = values[2 * N]

        std_ext = (4.0 * val_h2 - val_h) / 3.0
        sing_ext = (factor * val_h2 - val_h) / (factor - 1.0)

        err = abs(val_h - exact_value)
        err_std = abs(std_ext - exact_value)
        err_sing = abs(sing_ext - exact_value)

        print(
            f"{N:6d} | "
            f"{val_h:14.10f} {err:12.4e} | "
            f"{std_ext:14.10f} {err_std:12.4e} | "
            f"{sing_ext:14.10f} {err_sing:12.4e}"
        )

    print()
    print("The standard h^2 extrapolation is not designed for h^(2 alpha) pollution.")
    print("The singular extrapolation uses the exponent q=2 alpha.")
    print()


# ============================================================
# Part C. residual / defect indicator for 2D P1 FEM
# ============================================================

def element_gradients(nodes, triangles, uh):
    """
    Compute constant gradients of a P1 function on each triangle.
    """
    grads = np.zeros((len(triangles), 2))
    areas = np.zeros(len(triangles))

    for it, tri in enumerate(triangles):
        coords = nodes[tri, :]
        x0, y0 = coords[0]
        x1, y1 = coords[1]
        x2, y2 = coords[2]

        detJ = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)
        area = abs(detJ) / 2.0

        b = np.array([y1 - y2, y2 - y0, y0 - y1]) / detJ
        c = np.array([x2 - x1, x0 - x2, x1 - x0]) / detJ

        grad = np.array([np.dot(uh[tri], b), np.dot(uh[tri], c)])

        grads[it] = grad
        areas[it] = area

    return grads, areas


def p1_residual_indicator(nodes, triangles, uh, f):
    """
    A simple residual/defect indicator for P1 FEM.

    Since the P1 function is linear on each element, Delta uh = 0
    inside each element. Thus the element residual is approximately f.
    The edge part uses jumps of the normal derivative across interior edges.

    This is not a full a posteriori estimator implementation.
    It is a lightweight diagnostic to visualise where the defect concentrates.
    """
    grads, areas = element_gradients(nodes, triangles, uh)

    eta2 = np.zeros(len(triangles))

    # element residual part: h_T^2 ||f||^2_T, approximated at centroid
    for it, tri in enumerate(triangles):
        coords = nodes[tri, :]
        xc = float(np.mean(coords[:, 0]))
        yc = float(np.mean(coords[:, 1]))
        hT = max(
            np.linalg.norm(coords[0] - coords[1]),
            np.linalg.norm(coords[1] - coords[2]),
            np.linalg.norm(coords[2] - coords[0]),
        )
        eta2[it] += (hT**2) * (f(xc, yc) ** 2) * areas[it]

    # edge map
    edge_to_tri = {}
    for it, tri in enumerate(triangles):
        for local_edge in [(0, 1), (1, 2), (2, 0)]:
            a = int(tri[local_edge[0]])
            b = int(tri[local_edge[1]])
            edge = tuple(sorted((a, b)))
            edge_to_tri.setdefault(edge, []).append(it)

    for edge, elems in edge_to_tri.items():
        if len(elems) != 2:
            continue

        i, j = elems
        p0 = nodes[edge[0]]
        p1 = nodes[edge[1]]
        tangent = p1 - p0
        length = np.linalg.norm(tangent)

        # A unit normal. The sign is irrelevant because we square the jump.
        normal = np.array([tangent[1], -tangent[0]]) / length

        jump = np.dot(grads[i] - grads[j], normal)
        contribution = 0.5 * length * jump**2

        eta2[i] += contribution / 2.0
        eta2[j] += contribution / 2.0

    return np.sqrt(eta2)


def plot_lshape_defect_indicator(N=64, filename="lshape_defect_indicator.png"):
    """
    Plot the residual/defect indicator for the singular L-shaped problem.
    """
    alpha = 2.0 / 3.0
    exact = singular_exact_factory(alpha)

    nodes, tris = generate_lshape_mesh(N)
    uh = solve_dirichlet_poisson(nodes, tris, zero_rhs, exact)

    eta = p1_residual_indicator(nodes, tris, uh, zero_rhs)

    triang = mtri.Triangulation(nodes[:, 0], nodes[:, 1], tris)

    plt.figure(figsize=(6, 5))
    plt.tripcolor(triang, facecolors=eta, shading="flat")
    plt.colorbar(label="residual / defect indicator")
    plt.gca().set_aspect("equal")
    plt.title("L-shaped domain: defect indicator")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.tight_layout()
    plt.savefig(filename, dpi=200)
    plt.close()

    print("=" * 80)
    print("Part C. Defect indicator plot")
    print(f"Wrote {filename}")
    print("The defect indicator should be large near the reentrant corner.")
    print("=" * 80)
    print()


def main():
    run_defect_correction_1d()
    run_point_extrapolation_smooth()
    run_point_extrapolation_lshape()
    plot_lshape_defect_indicator(N=64)


if __name__ == "__main__":
    main()