""" Scientific ML Studio -- reference code (written by hand, not generated by the Studio) Licence : Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) https://creativecommons.org/licenses/by-nc/4.0/ (c) 2026 Scientific ML Studio. All rights are reserved to the Scientific ML Studio platform except the permissions granted by the licence above. You may share and adapt this code for NON-COMMERCIAL purposes, provided you credit "Scientific ML Studio" with a link to https://scimlstudio.com/, and indicate any changes you make. Commercial use, paid client work and production deployment need prior written permission, which you can ask for through the site's Contact page. ------------------------------------------------------------------------------- Reference solution for the "Heat in a plate" worked problem. u_t = alpha (u_xx + u_yy) on the unit square, 0 < t <= 1, alpha = 0.05 u = 0 on x = 0, y = 0, y = 1; du/dx = 0 on x = 1 u(x, y, 0) = exp(-50 [(x - 0.35)^2 + (y - 0.5)^2]) Two independent solutions are provided: 1. fem_solve() Finite elements. Quadratic (Q2) Lagrange elements on a uniform N x N mesh of squares, built as tensor products of 1-D quadratic elements, so the 2-D matrices are Kronecker products of 1-D ones. Time stepping: Crank-Nicolson with two implicit-Euler half steps at the start (Rannacher start-up). 2. series() Exact eigenfunction series. The initial condition is a product f(x) g(y) and the boundary conditions are separable, so u = 4 X(x, t) Y(y, t) with sin((k + 1/2) pi x) in x and sin(n pi y) in y. Licence: CC BY-NC 4.0, (c) 2026 Scientific ML Studio. """ from __future__ import annotations import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla ALPHA = 0.05 X0, Y0, K = 0.35, 0.5, 50.0 def ic(x, y): return np.exp(-K * ((x - X0) ** 2 + (y - Y0) ** 2)) # ---------------------------------------------------------------- 1-D P2 pieces _G = np.array([-np.sqrt(3 / 5), 0.0, np.sqrt(3 / 5)]) # 3-point Gauss, enough for P2*P2 _W = np.array([5 / 9, 8 / 9, 5 / 9]) def _shape(xi): """Quadratic Lagrange basis on [-1, 1] and its derivative.""" xi = np.asarray(xi, float) phi = np.stack([0.5 * xi * (xi - 1), 1 - xi ** 2, 0.5 * xi * (xi + 1)], axis=-1) dphi = np.stack([xi - 0.5, -2 * xi, xi + 0.5], axis=-1) return phi, dphi def matrices_1d(n_el): """Consistent mass and stiffness on [0, 1], n_el quadratic elements, 2 n_el + 1 nodes.""" h = 1.0 / n_el phi, dphi = _shape(_G) # (3 gauss, 3 shape) m_loc = (phi.T * _W) @ phi * (h / 2) k_loc = (dphi.T * _W) @ dphi * (2 / h) nn = 2 * n_el + 1 rows, cols, mv, kv = [], [], [], [] for e in range(n_el): idx = np.array([2 * e, 2 * e + 1, 2 * e + 2]) rows += list(np.repeat(idx, 3)); cols += list(np.tile(idx, 3)) mv += list(m_loc.ravel()); kv += list(k_loc.ravel()) M = sp.csr_matrix((mv, (rows, cols)), shape=(nn, nn)) Kk = sp.csr_matrix((kv, (rows, cols)), shape=(nn, nn)) return M, Kk def eval_basis_1d(pts, n_el): """For points in [0,1]: (first node index, 3 basis values) for the element containing each.""" pts = np.clip(np.asarray(pts, float), 0.0, 1.0) e = np.minimum((pts * n_el).astype(int), n_el - 1) xi = 2 * (pts * n_el - e) - 1 phi, _ = _shape(xi) return 2 * e, phi # --------------------------------------------------------------------- FEM solve def fem_solve(n_el=64, n_steps=1000, t_end=1.0): """Return a dict with the nodal solution at every step. Shape of U: (n_steps+1, nx, ny).""" M1, K1 = matrices_1d(n_el) nn = 2 * n_el + 1 M = sp.kron(M1, M1, format="csr") K2 = (sp.kron(K1, M1) + sp.kron(M1, K1)).tocsr() ix = np.arange(1, nn) # x = 0 fixed; x = 1 free (natural zero-flux condition) iy = np.arange(1, nn - 1) # y = 0 and y = 1 fixed free = (ix[:, None] * nn + iy[None, :]).ravel() Mf = M[free][:, free].tocsc() Kf = K2[free][:, free].tocsc() xs = np.linspace(0, 1, nn) X, Y = np.meshgrid(xs, xs, indexing="ij") u = ic(X, Y).ravel()[free] # nodal interpolation of the initial condition dt = t_end / n_steps # Rannacher start-up: two implicit-Euler steps of size dt/2, then Crank-Nicolson with step dt. # Both use the same left-hand matrix M + (dt/2) alpha K, so one factorisation serves both. lu = spla.splu((Mf + 0.5 * dt * ALPHA * Kf).tocsc()) hist = np.zeros((n_steps + 1, nn * nn)) full = np.zeros(nn * nn); full[free] = u; hist[0] = full # step 1 uses two BE half-steps (dt/2 each, matrix M + dt/2 a K) for _ in range(2): u = lu.solve(Mf @ u) full = np.zeros(nn * nn); full[free] = u; hist[1] = full rhs_mat = (Mf - 0.5 * dt * ALPHA * Kf).tocsr() for n in range(2, n_steps + 1): u = lu.solve(rhs_mat @ u) full = np.zeros(nn * nn); full[free] = u; hist[n] = full return dict(n_el=n_el, n_steps=n_steps, dt=dt, nn=nn, U=hist.reshape(n_steps + 1, nn, nn)) def fem_eval(sol, x, y, t): """Evaluate the FEM solution at arbitrary points (x, y, t): quadratic in space, linear in time.""" n_el, nn, U, dt = sol["n_el"], sol["nn"], sol["U"], sol["dt"] x, y, t = (np.asarray(a, float).ravel() for a in (x, y, t)) ix0, px = eval_basis_1d(x, n_el) iy0, py = eval_basis_1d(y, n_el) s = np.clip(t / dt, 0, sol["n_steps"] - 1e-12) n0 = s.astype(int); w = s - n0 out = np.zeros_like(x) for a in range(3): for b in range(3): coef = px[:, a] * py[:, b] u0 = U[n0, ix0 + a, iy0 + b] u1 = U[n0 + 1, ix0 + a, iy0 + b] out += coef * ((1 - w) * u0 + w * u1) return out # ------------------------------------------------------------------ exact series def series(x, y, t, n_modes=400): """Eigenfunction series. u = 4 X(x,t) Y(y,t) because the initial condition is a product.""" x, y, t = (np.asarray(a, float).ravel() for a in (x, y, t)) gx, gw = np.polynomial.legendre.leggauss(600) s = 0.5 * (gx + 1); w = 0.5 * gw # quadrature on [0, 1] k = np.arange(n_modes) a = (k + 0.5) * np.pi n = np.arange(1, n_modes + 1) * np.pi fx = np.exp(-K * (s - X0) ** 2) fy = np.exp(-K * (s - Y0) ** 2) Ix = (np.sin(np.outer(a, s)) * (fx * w)).sum(axis=1) # int f(x) sin(a_k x) dx Iy = (np.sin(np.outer(n, s)) * (fy * w)).sum(axis=1) X = (Ix[None, :] * np.sin(np.outer(x, a)) * np.exp(-ALPHA * np.outer(t, a ** 2))).sum(axis=1) Y = (Iy[None, :] * np.sin(np.outer(y, n)) * np.exp(-ALPHA * np.outer(t, n ** 2))).sum(axis=1) return 4.0 * X * Y