PHYSICS · Multi-part · Fermi University Physics · Mathematical Methods for Physics · Computational Physics · Research software practice, verification and validation
Scientific Computing for Physicists, Ch Is the equation solved right? Testing, debugging and code verification
Problem
An explicit diffusion code solves ∂tu=D∂x2u on [0,1] with D=0.5m2s−1, diffusion number α=DΔt/Δx2=0.4 held fixed, and zero-flux (insulating) boundaries at both ends imposed with mirror ghost nodes. It starts from u(x,0)=1+cos(πx) and runs to t=0.08s; the exact solution u=1+e−Dπ2tcos(πx) has zero slope at both ends, so it satisfies the boundary conditions exactly. Maximum errors against it are measured as e(N=50)=1.23×10−4 and e(N=100)=3.06×10−5, where N is the number of grid intervals.
((a)) Using the two measured errors and the fact that this scheme's error scales as Δx2∝N−2, predict the maximum error at N=200. A unit test at N=200 compares the computed field with the exact solution; set its tolerance at twice the prediction, and say whether a measured error of 7.66×10−6 passes it.
((b)) Because the boundaries are insulating, ∫01udx is conserved exactly in time. A property test integrates u with the trapezium rule on the N+1=401 grid points at N=400 and checks this integral against its initial value to relative tolerance 10(N+1)ϵmach/2, with ϵmach=2−52. Compute the tolerance, and say whether a measured relative change of 4.4×10−16 passes it.
((c)) Errors are also measured at N=200 and N=400: e(200)=7.66×10−6 and e(400)=1.91×10−6. Using all four grids N=50,100,200,400, compute the observed order p^ from each consecutive pair via p^=ln[e(N)/e(2N)]/ln2, and say whether the code is verified against the formal order pformal=2.