import json import math from pathlib import Path import numpy as np def pade_target(s, z): p = sum(math.factorial(2*s-k)*math.factorial(s-1) / (math.factorial(2*s)*math.factorial(k)*math.factorial(s-1-k)) * z**k for k in range(s)) q = sum(math.factorial(2*s-k)*math.factorial(s+1) / (math.factorial(2*s)*math.factorial(k)*math.factorial(s+1-k)) * (-z)**k for k in range(s+2)) return p / q def make_coeffs(s=2, nodes=None): if nodes is None: nodes = np.linspace(0.0, 1.0, s + 1) q = np.zeros(2*s + 1) for k in range(s + 2): q[k] = (math.factorial(2*s-k)*math.factorial(s+1) / (math.factorial(2*s)*math.factorial(k)*math.factorial(s+1-k)) * (-1)**k) pi = np.zeros((s+1, 2*s + 1)) for j, c in enumerate(nodes): for r in range(2*s + 1): pi[j, r] = sum(q[k] * c**(r-k) / math.factorial(r-k) for k in range(min(s+1, r)+1)) # B represents Pi_j and z Pi_j in P_{2s+1}. basis = [np.r_[pi[j], 0.0] for j in range(s+1)] basis += [np.r_[0.0, pi[j]] for j in range(s+1)] B = np.column_stack(basis) aa = np.zeros((s, s+1)); hh = np.zeros_like(aa) for i in range(1, s+1): # (Pi_i-Q)/z: constant difference is zero; represent through degree 2s+1. rhs = np.r_[pi[i, 1:] - q[1:], 0.0, 0.0] sol = np.linalg.solve(B, rhs) aa[i-1] = sol[:s+1] hh[i-1] = sol[s+1:] return np.asarray(nodes), aa, hh, q def step_ph(u, t, h, fun, jac, A, H, nodes, tol=1e-11, maxit=100): s = len(nodes)-1; u = np.asarray(u, float) U = np.array([u + nodes[i]*h*fun(t, u) for i in range(1, s+1)]) for it in range(maxit): states = np.vstack([u, U]) fs = np.array([fun(t + nodes[j]*h, states[j]) for j in range(s+1)]) gs = np.array([jac(t + nodes[j]*h, states[j]) @ fs[j] for j in range(s+1)]) new = np.array([u + h*(A[i] @ fs) + h*h*(H[i] @ gs) for i in range(s)]) err = np.max(np.abs(new-U)) / max(1.0, np.max(np.abs(new))) U = 0.7*new + 0.3*U if err < tol: return U[-1], it+1 raise RuntimeError('nonlinear stage solve failed') def ph_integrate(u, t0, t1, h, fun, jac, A, H, nodes): n = int(round((t1-t0)/h)); x = np.asarray(u, float); t=t0; its=0 for _ in range(n): x, k = step_ph(x, t, h, fun, jac, A, H, nodes) its += k; t += h return x, its def rk4(u, t, h, f): k1=f(t,u); k2=f(t+h/2,u+h*k1/2); k3=f(t+h/2,u+h*k2/2); k4=f(t+h,u+h*k3) return u+h*(k1+2*k2+2*k3+k4)/6 def main(): np.random.seed(3) nodes,A,H,q = make_coeffs(2) zs = np.array([-1., -5., -10., -100., -1000.]) pade_err = [abs(pade_target(2,z) - (1+z/4)/(1-3*z/4+z*z/4-z**3/24)) for z in zs] order_rows=[] for h in [0.2,0.1,0.05,0.025]: f=lambda t,x: np.array([-x[0]]) j=lambda t,x: np.array([[-1.]]) y,_=ph_integrate([1.],0.,1.,h,f,j,A,H,nodes) order_rows.append((h, float(abs(y[0]-math.exp(-1))))) orders=[math.log(a[1]/b[1],2) for a,b in zip(order_rows[:-1],order_rows[1:])] stability=[] for z in [-1,-2,-3,-5,-10,-50,-100]: r=pade_target(2,z) stability.append({'z':z,'R':float(r),'euler_abs':abs(1+z)}) lam=80. def f(t,x): return np.array([-lam*x[0] - 0.2*x[1]**3, -x[1] + 0.1*x[0]]) def jac(t,x): return np.array([[-lam, -0.6*x[1]**2],[0.1,-1.]]) x0=np.array([1.,1.]); exact_ref,_=ph_integrate(x0,0.,1.,0.002,f,jac,A,H,nodes) cmp=[] for h in [0.05,0.025,0.0125]: try: yp,it=ph_integrate(x0,0.,1.,h,f,jac,A,H,nodes) pe=float(np.linalg.norm(yp-exact_ref)) except Exception: pe=float('inf'); it=-1 y=x0.copy(); ok=True for k in range(int(round(1/h))): y=rk4(y,k*h,h,f) if not np.all(np.isfinite(y)): ok=False cmp.append({'h':h,'ph_error':pe,'ph_iterations':it, 'rk4_error':float(np.linalg.norm(y-exact_ref)) if ok else float('inf')}) # Quantitative math checks: Hermite moments and A-stability grid. moment_residual=0.0 for i,c in enumerate(nodes[1:]): for m in range(1,5): lhs=sum(A[i,j]*nodes[j]**(m-1)/math.factorial(m-1) + H[i,j]*(m-1)*nodes[j]**(m-2)/math.factorial(m-1) if m>=2 else A[i,j]*nodes[j]**(m-1)/math.factorial(m-1) for j in range(3)) # Integral of t^(m-1) from 0 to c, with derivative-data term target=c**m/math.factorial(m) moment_residual=max(moment_residual,abs(lhs-target)) grid=[] for re in np.linspace(-20,0,41): for im in np.linspace(-20,20,81): grid.append(abs(pade_target(2,re+1j*im))) a_stab_max=max(grid) # One-step amplification from the actual anchored stages, across node choices. node_sweep=[] for nd in [np.array([0.,.5,1.]), np.array([0.,.25,1.]), np.array([0.,.75,1.])]: nn,aa,hh,_=make_coeffs(2,nd) vals=[] for z in [-1.,-10.,-100.]: # Exact solve of the scalar linear stage equations, avoiding # convergence assumptions in this stability-only diagnostic. c=nn; m=len(c)-1 M=np.eye(m)-z*aa[ :,1:]-z*z*hh[:,1:] rhs=np.ones(m)+z*aa[:,0]+z*z*hh[:,0] yy=np.linalg.solve(M,rhs)[-1] vals.append(float(yy)) node_sweep.append({'nodes':nd.tolist(),'R':vals}) # Predicted explicit-Euler boundary is z=-2; measure the transition. euler_boundary=[] for z in [-1.9,-2.0,-2.1]: euler_boundary.append({'z':z,'abs_R':abs(1+z)}) asymptotic_slope=(math.log(abs(pade_target(2,-1000.)/pade_target(2,-100.)),10)) out={'s':2,'nodes':nodes.tolist(),'A':A.tolist(),'H':H.tolist(), 'pade_negative_axis_errors':pade_err,'order_errors':order_rows, 'observed_orders':orders,'stability_sweep':stability, 'node_independence':node_sweep,'euler_boundary':euler_boundary, 'observed_log10_decay_slope_100_to_1000':asymptotic_slope, 'max_moment_residual':moment_residual,'max_left_half_plane_abs_R':float(a_stab_max), 'stiff_comparison':cmp} Path('results.json').write_text(json.dumps(out,indent=2)); print(json.dumps(out,indent=2)) if __name__=='__main__': main()