Google Colabで始めるFreeFEM入門：Poisson方程式，誤差評価，log-logプロット

%%bash
set -e

apt-get update -y
apt-get install -y wget ca-certificates \
  libgsl-dev libhdf5-dev liblapack-dev libopenmpi-dev freeglut3-dev

source /etc/os-release
echo "Ubuntu version: $VERSION_ID"

if [ "$VERSION_ID" = "24.04" ]; then
  URL="https://github.com/FreeFem/FreeFem-sources/releases/download/v4.15/FreeFEM-4.15-amd64-ubuntu24.04.deb"
elif [ "$VERSION_ID" = "22.04" ]; then
  URL="https://github.com/FreeFem/FreeFem-sources/releases/download/v4.15/FreeFEM-4.15-amd64-ubuntu22.04.deb"
else
  echo "Unsupported Ubuntu version in this Colab runtime: $VERSION_ID"
  exit 1
fi

wget -q "$URL" -O /tmp/freefem.deb
dpkg -i /tmp/freefem.deb || apt-get -f install -y

echo "Installed binaries:"
which FreeFem++ || true
which FreeFem++-nw || true
ls -l /usr/local/bin/FreeFem* || true

%%bash
cat > poisson_error_demo.edp <<'EOF'
load "iovtk"

real pi2 = pi*pi;
int[int] Ns = [8, 16, 32, 64, 128];
ofstream ferr("errors.csv");
ferr << "N,h,absL2,relL2,absH1semi,relH1semi,absH1\n";

func sol  = sin(pi*x)*sin(pi*y);
func solx = pi*cos(pi*x)*sin(pi*y);
func soly = pi*sin(pi*x)*cos(pi*y);
func rhs  = 2*pi2*sin(pi*x)*sin(pi*y);

for (int k = 0; k < Ns.n; ++k) {
  int N = Ns[k];
  mesh Th = square(N, N);

  fespace Vh(Th, P1);
  Vh uh, v;

  problem Poisson(uh, v, solver=LU)
      = int2d(Th)( dx(uh)*dx(v) + dy(uh)*dy(v) )
      - int2d(Th)( rhs*v )
      + on(1, 2, 3, 4, uh=0);

  Poisson;

  real absL2 = sqrt(int2d(Th, qforder=10)((uh - sol)^2));
  real absH1semi = sqrt(int2d(Th, qforder=10)(
        (dx(uh) - solx)^2
      + (dy(uh) - soly)^2
  ));
  real absH1 = sqrt(int2d(Th, qforder=10)(
        (uh - sol)^2
      + (dx(uh) - solx)^2
      + (dy(uh) - soly)^2
  ));

  real normL2 = sqrt(int2d(Th, qforder=10)(sol^2));
  real seminormH1 = sqrt(int2d(Th, qforder=10)(solx^2 + soly^2));
  real relL2 = absL2 / normL2;
  real relH1semi = absH1semi / seminormH1;
  real h = sqrt(2.0) / N;

  cout << "N = " << N
       << ", h = " << h
       << ", absL2 = " << absL2
       << ", relL2 = " << relL2
       << ", absH1semi = " << absH1semi
       << ", relH1semi = " << relH1semi << endl;

  ferr << N << ","
       << h << ","
       << absL2 << ","
       << relL2 << ","
       << absH1semi << ","
       << relH1semi << ","
       << absH1 << "\n";
}

int Nvis = 16;
mesh ThVis = square(Nvis, Nvis);

fespace VhVis(ThVis, P1);
VhVis uhVis, vVis;

problem PoissonVis(uhVis, vVis, solver=LU)
    = int2d(ThVis)( dx(uhVis)*dx(vVis) + dy(uhVis)*dy(vVis) )
    - int2d(ThVis)( rhs*vVis )
    + on(1, 2, 3, 4, uhVis=0);

PoissonVis;

int[int] order = [1];
string dataName = "uh";
savevtk("uh_vis.vtu", ThVis, uhVis, dataname=dataName, order=order);
EOF

FreeFem++-nw poisson_error_demo.edp

import numpy as np
import matplotlib.pyplot as plt

# CSV を読み込みます
raw = np.genfromtxt("errors.csv", delimiter=",", names=True)

Ns = raw["N"]
hs = raw["h"]
abs_L2s = raw["absL2"]
rel_L2s = raw["relL2"]
abs_H1s = raw["absH1semi"]
rel_H1s = raw["relH1semi"]

def rates(errors, hs):
    return np.log(errors[1:] / errors[:-1]) / np.log(hs[1:] / hs[:-1])

print("L2 absolute rates  =", rates(abs_L2s, hs))
print("L2 relative rates  =", rates(rel_L2s, hs))
print("H1 absolute rates  =", rates(abs_H1s, hs))
print("H1 relative rates  =", rates(rel_H1s, hs))

plt.figure(figsize=(6, 5))
plt.loglog(hs, abs_L2s, "o-", label="absolute L2 error")
plt.loglog(hs, rel_L2s, "s--", label="relative L2 error")
plt.xlabel("h")
plt.ylabel("error")
plt.title("L2 errors on a log-log plot")
plt.grid(True, which="both")
plt.legend()
plt.show()

plt.figure(figsize=(6, 5))
plt.loglog(hs, abs_H1s, "o-", label="absolute H1 seminorm error")
plt.loglog(hs, rel_H1s, "s--", label="relative H1 seminorm error")
plt.xlabel("h")
plt.ylabel("error")
plt.title("H1 seminorm errors on a log-log plot")
plt.grid(True, which="both")
plt.legend()
plt.show()

import meshio
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

m = meshio.read("uh_vis.vtu")
points = m.points[:, :2]

triangles = None
for cell_block in m.cells:
    if cell_block.type == "triangle":
        triangles = cell_block.data
        break

if triangles is None:
    raise RuntimeError("triangle cells were not found in uh_vis.vtu")

if "uh" in m.point_data:
    values = np.asarray(m.point_data["uh"]).reshape(-1)
else:
    raise RuntimeError("point data 'uh' was not found in uh_vis.vtu")

plt.figure(figsize=(6, 6))
plt.triplot(points[:, 0], points[:, 1], triangles, linewidth=0.8)
plt.xlabel("x")
plt.ylabel("y")
plt.title("2D mesh")
plt.gca().set_aspect("equal")
plt.show()

plt.figure(figsize=(6, 5))
plt.tricontourf(points[:, 0], points[:, 1], triangles, values, levels=20)
plt.colorbar(label="u_h")
plt.xlabel("x")
plt.ylabel("y")
plt.title("2D contour of FEM solution")
plt.gca().set_aspect("equal")
plt.show()

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection="3d")
ax.plot_trisurf(points[:, 0], points[:, 1], values, triangles=triangles,
                linewidth=0.2, antialiased=True)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u_h")
ax.set_title("3D surface of FEM solution")
plt.show()