Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 13

Optimization and Data Fitting

These additional exercises and physics projects accompany Chapter 13 of Computational Physics: From Equations to Simulation. They explore scaling, stationary points, search reliability, measurement errors, and the information available in a data set. The physics projects introduce new models and independent checks. Use the companion codes or suitable numerical libraries in your preferred language.

Level guide: Foundation uses a supplied routine with a direct check; Intermediate adds repeated runs or diagnostic analysis; Advanced combines methods or requires more statistical or physical interpretation. Extensions may be more demanding than the core task.

Chapter 13 codes Chapter 13 resources Chapter 13 updates & errata Chapter 13 overview

Additional exercises

  1. Objective scaling and stopping criteria

    Level: Foundation

    Minimize fs(x) = s(x − 0.37)2 on [−1,1] for s = 10−6, 1, 106. The minimizer is unchanged by this scaling. Compare a golden-section search with a local gradient-based routine starting at x = 0.9.

    Request an interval-width tolerance of 10−7 for the first method and a gradient-norm tolerance of 10−6 for the second. Report the position error, objective value, gradient, stopping reason, and objective-evaluation count. Explain why the same absolute gradient tolerance gives different position accuracy as s changes. Then rescale the objective or tolerance to obtain comparable position accuracy.

    Extension: repeat the golden-section experiment with gs(x) = 1 + s(x − 0.37)2 and requested widths from 10−6 down to 10−14. Compare actual location errors with reported bracket widths. Near the minimum, rounded values can become indistinguishable when the quadratic increment is of order machine precision. Estimate the position scale sqrt(eps/s) and test how it varies with s. This is a scale estimate, not a universal accuracy bound; distinguish objective rounding from eventual stagnation of the interval endpoints.

  2. A stationary point that is not a minimum

    Level: Foundation

    f(x,y) = x2/2 + y4/4 − y2/2.

    Plot contours on [−2,2]2 and run gradient descent from (0,0), (0.1,0.01), and (0.1,−0.01). Supply the analytic gradient (x, y3 − y). Compare the reported termination status with the local shape of the surface.

    Check the Hessian at the origin and at the two minima (0,1) and (0,−1), where f = −1/4. Explain why starting exactly at the saddle can satisfy a gradient stopping rule immediately. Design a small perturbation check that would reveal the problem.

  3. Several equally good solutions

    Level: Intermediate

    f(x,y) = (x2 + y − 11)2 + (x + y2 − 7)2.

    This is the Himmelblau test function. Search [−6,6]2 using particle swarm or differential evolution. Use 20 recorded seeds and an objective budget of 5000 evaluations per run. Count all initialization and search evaluations. Keep population size and other settings fixed within each comparison. With SciPy differential_evolution, set polish=False for this comparison or count polishing within the same total budget. Set the generation limit to respect the actual population size and initialization cost.

    The function is nonnegative and has f(3,2) = 0. Classify a run as successful when its best value is below 10−6. Report the success fraction, median best value, and distinct solution clusters. Check other near-zero candidates by substituting them into both squared terms. Explain why different coordinates need not mean that one run has failed. As an extension, compare two methods at the same evaluation budget and report any local-refinement cost separately.

  4. When measurement errors vary across the data

    Level: Intermediate

    Generate 40 equally spaced x values on [0,5] and independent observations yi = 2 − 0.4xi + sigmaizi, with sigmai = 0.05 + 0.1xi and standard normal zi. Fit the same data using ordinary and weighted least squares. Use weights 1/sigmai2.

    Plot both raw and standardized residuals. Treat the supplied standard deviations as known when estimating the weighted-fit covariance. Repeat for 200 seeds and compare the bias and spread of the intercept and slope estimates. Does weighting reduce the observed parameter spread? Explain why it need not improve every individual noisy realization. Compare the usual homoscedastic OLS standard errors and the weighted-fit standard errors with the empirical spread across seeds. The usual OLS covariance assumes a common noise variance, which this data model violates. As an extension, compare a heteroscedasticity-robust OLS covariance estimate.

  5. A good curve with poorly determined parameters

    Level: Advanced

    y(t) = A exp(−t/τ) + B.

    Generate 40 independent noisy measurements with A = 2, τ = 1, B = 0.5, and Gaussian noise of standard deviation 0.02. Compare two observation windows: [0,0.2] and [0,5]. Fit all three parameters with τ > 0 from several initial guesses.

    Report residual sums, parameter standard errors, and parameter correlations where the covariance estimate is reliable. Examine the singular values of the residual Jacobian. Use the short-time approximation y(t) ≈ A + B − (A/τ)t to explain which parameter combinations are measured in the short window. Compare predicted curves and parameter estimates across several seeds. Report fits that hit evaluation limits, drift to extreme parameter values, or reach any imposed bounds; give the bounds explicitly. Nearly fixed A + B and A/τ allow A and τ to grow while B becomes very negative, with little change in the short-window curve. What additional measurements would help?

  6. Bootstrap uncertainty in a derived quantity

    Level: Advanced

    Generate 25 independent x values uniformly on [−1,1], then generate y = 0.5 + x + epsilon with independent Gaussian noise of standard deviation 0.2. Fit y = a + bx and estimate its zero crossing xc = −a/b.

    Use 2000 pairs-bootstrap resamples, keeping each (xi,yi) pair together. Refit every resample and calculate its zero crossing. Plot the distribution and report its standard deviation and 2.5th and 97.5th percentiles. Compare with uncertainty propagation using the full covariance of (a,b), including the intercept–slope covariance.

    Repeat with true slope b = 0.1. Explain why a slope near zero can produce extreme zero-crossing estimates and why the new true crossing lies outside the observed interval. Report singular fits and near-zero slopes rather than silently removing inconvenient resamples. Increase the resample count to check whether each uncertainty summary is stable.

    Advanced extension: construct a Fieller confidence set for −a/b using the joint coefficient covariance. Allow an unbounded or disconnected set rather than forcing a finite interval. Compare its interpretation with a propagated standard error and bootstrap percentiles when the slope is poorly separated from zero. See E. C. Fieller, Some Problems in Interval Estimation (1954).

  7. A closed-tour benchmark for heuristic methods

    Level: Intermediate; comparing all three methods is an advanced extension.

    Place N = 12 cities at equally spaced angles on a unit circle: (xj,yj) = (cos(2πj/N),sin(2πj/N)), with j = 0,...,N − 1. Randomly permute the city labels with a recorded seed so that label order does not reveal the perimeter tour. Find the shortest closed tour using ant colony optimization, simulated annealing, or a genetic algorithm.

    Represent each candidate as a permutation. For simulated annealing, use swaps or reversals of tour segments; for a genetic algorithm, use permutation-preserving crossover and mutation. The continuous or bit-string genetic representation must be adapted to this task. Start with one method, then compare the others over 20 seeds with at most 10000 complete-tour length evaluations per run. Count initialization and all search costs.

    Lopt = 2N sin(π/N),   Lopt ≈ 6.211657082 for N = 12.

    Derive this known optimum from the regular polygon: every tour edge is at least the nearest-neighbor chord length, and perimeter order attains that bound. Include the last-to-first edge in every tour. Check that each city occurs exactly once, plot the best tour, and report relative excess length (L − Lopt)/Lopt. Count a run as successful when this excess is below 10−8. Compare success rates and evaluation counts, treating rotations and reversals as the same tour.

Additional physics projects

The first two projects use optimization to estimate equilibrium states or energies. The remaining three infer physical parameters from simulated measurements. Use the stated values for an initial calculation, then explore the suggested extensions.

1. Equilibrium of coupled nonlinear springs

Level: Intermediate

Objective: find equilibrium by minimizing energy and verify it through force balance.

Two particles move along a line, with displacements x1 and x2 from their reference positions. The first is attached to a wall by a linear spring. A nonlinear spring connects the particles, and a constant force F acts on the second particle. In reduced units, let

U = k1x12/2 + k2(x2 − x1)2/2 + β(x2 − x1)4/4 − Fx2.

Use k1 = 10, k2 = 5, β = 20, and initially F = 1. Minimize from (0,0), (1,1), and (−1,2). Derive the gradient, compare methods with and without supplied derivatives, and report displacements, energy, force residuals, and evaluation counts.

Validation: writing d = x2 − x1, equilibrium requires k1x1 = F and k2d + βd3 = F. Independently solve the scalar equation. The default gives x1 = 0.1 and x2 ≈ 0.277594729. Show that the Hessian is positive definite everywhere for these coefficients, so the energy is strictly convex and has a unique minimum. All successfully converged minimizations should agree.

Possible extension: vary F from 0 to 5 and plot the force–displacement response. Compare with the linear limit β = 0 and use each equilibrium as the next starting guess.

Advanced extension: set k2 = −5, keeping β = 20. This models an effective bistable coupling. Find all roots of −5d + 20d3 = F and classify their stability. There are three distinct equilibria for |F| < 0.962250449. Follow a stable local minimum while slowly sweeping F from −1.2 to 1.2 and back. Use the previous solution as the starting point and monitor curvature. Compare the resulting metastable branches and snap-through with a search that selects the global minimum at every force. Explain why local branch following can show hysteresis while global minimization does not give the same loop.

2. A variational estimate for an anharmonic oscillator

Level: Advanced

Objective: optimize a trial wavefunction and assess the limits of its flexibility.

In units with hbar = m = ω = 1, consider H = −(1/2)d2/dx2 + x2/2 + λx4, with λ ≥ 0. Use the normalized Gaussian trial function

ψa(x) = (πa2)−1/4 exp(−x2/(2a2)),   a > 0.

Evaluate its kinetic and potential energy integrals numerically, then minimize their sum over a > 0. Start with λ = 0 and 0.1; repeat for 0.5 and 1. If using a finite integration interval, enlarge it and refine the quadrature until the optimized energy is stable. Compare the changes in width and energy with increasing anharmonicity.

Validation: derive the Gaussian expectation E(a) = 1/(4a2) + a2/4 + 3λa4/4 and use it as an independent check. At λ = 0, the optimum is a = 1, E = 1/2. At λ = 0.1, expect a ≈ 0.90491376 and E ≈ 0.56030737 within this trial family. This optimized trial energy is an upper bound to the exact ground-state energy.

Possible extension: compare with an independent finite-difference eigenvalue calculation. For λ = 0.1, a converged numerical reference is E0 ≈ 0.55914633, about 0.21% below the Gaussian estimate. Establish mesh and domain convergence rather than treating a single discretized eigenvalue as exact. Separate the error of the trial family from quadrature and minimization errors.

3. Measuring damping and frequency from a noisy signal

Level: Advanced

Objective: estimate an oscillation frequency and decay rate while testing sensitivity to initial guesses and sampling.

x(t) = A exp(−γt) cos(ωt + φ) + c.

Generate 81 equally spaced observations on [0,8] with A = 1.2, γ = 0.15, ω = 2π(0.8), φ = 0.3, and c = 0.1, in consistent reduced units. Add independent Gaussian noise with standard deviation 0.05. Record the seed and fit all five parameters.

Use A > 0, γ ≥ 0, and 0 < ω < 10. Try starting angular frequencies 2, 4, 6, 8, with starting values A = 1, γ = 0.1, φ = 0, c = 0. Reduce fitted phases to a common interval when comparing solutions. Plot the fitted signal and residuals and compare local fits with a global search followed by local refinement.

Validation: first fit noise-free data and recover a zero-residual solution with the stated parameters, allowing phase equivalence modulo 2π. For noisy data, report parameter standard errors and distinguish angular frequency ω from frequency ω/(2π). Compare residual structure across different fitted frequency basins.

Possible extension: fit only the first two time units or retain every fifth observation. For every fifth observation starting at t = 0, the spacing is Δt = 0.5. Derive and identify the alias ω′ = 2π/Δt − ω ≈ 7.53982237, with φ′ = −φ. Show that the true and aliased parameters reproduce exactly the same model samples on this grid, with the same amplitude, damping, and offset. Both frequencies lie inside the stated bounds, so fitting these samples alone cannot distinguish them. Explain how changing the sampling or adding prior frequency information can resolve the ambiguity.

4. Temperature and calibration from a thermal spectrum

Level: Intermediate

Objective: infer temperature and an unknown instrument scale from spectral data.

Bλ(λ,T) = (2hc2/λ5) / [exp(hc/(λkBT)) − 1].

Model the measured signal as s(λ) = q Bλ(λ,T)/Bref, where q is a dimensionless calibration scale and Bref = Bλ(1000 nm,3000 K) is a fixed normalization. Generate 60 equally spaced wavelengths from 400 to 2000 nm with T = 3000 K, q = 0.8, and independent Gaussian errors of standard deviation 0.02 in the normalized signal.

Convert wavelengths to meters inside the model. Use h = 6.62607015 × 10−34 J s, c = 299792458 m/s, and kB = 1.380649 × 10−23 J/K. Fit T > 0 and q > 0 from several guesses, for example (2000,1), (4000,0.5), and (6000,0.2). Keep Bref fixed throughout the optimization. Improve parameter scaling by fitting T/1000 rather than T, or use a solver’s parameter-scale option, such as SciPy least_squares with suitable x_scale values.

Validation: recover the stated parameters from noise-free data. For noisy data, compare residuals and parameter covariance using the known noise level. Repeat with 60 observations restricted to [1500,2000] nm and compare the temperature–scale correlation and uncertainty with those from the full spectrum. Include several seeds before drawing conclusions about precision. At 400 nm, the true normalized signal is only about 0.0582, compared with noise standard deviation 0.02. Examine the model sensitivities to temperature and calibration, divided by the noise standard deviation, across wavelength. Explain why signal amplitude alone does not determine how much a measurement tells us about temperature.

Possible extension: add and fit a constant background. Investigate whether the chosen wavelength range provides enough information to separate temperature, calibration, and background.

5. Inferring a slit width from a diffraction pattern

Level: Intermediate

Objective: fit a diffraction model and connect the inferred width with the positions of intensity minima.

I(θ) = I0[sin(u)/u]2 + B,  u = πa sin(θ)/λ.

This is the far-field intensity pattern of a uniformly illuminated single slit. Define sin(u)/u = 1 at u = 0. Take known wavelength λ = 632.8 nm and true values a = 0.100 mm, I0 = 1, B = 0.05. Generate 101 equally spaced angles on [−0.025,0.025] rad, adding independent Gaussian noise of standard deviation 0.02 in intensity units.

Fit a, I0, and B with 0.05 ≤ a ≤ 0.20 mm and positive I0. Use several starting widths and consistent units. Plot the fitted pattern, its background level, and the residuals. Use the known noise to estimate parameter standard errors and repeat for several seeds.

Validation: the first minima satisfy a sin(θ) = ±λ. For the true width, their angles are approximately ±0.00632804 rad. Compare a width estimate based on the minima with the full-curve fit. Check noise-free recovery and handle the u = 0 limit explicitly; the specified angle grid includes θ = 0.

Possible extension: introduce a small unknown angular offset and fit it as a fourth parameter. Compare fits using the entire pattern and only the central peak, with attention to parameter correlations.

What to report

State the model, units, method, settings, tolerances, and random seeds. Include the plots needed to interpret the result, the stopping status, objective-evaluation count, and at least one independent check. For fits, show residuals and state whether reported uncertainties are standard errors or confidence intervals. The Chapter 13 resources explain covariance scaling, constraints, and bootstrap choices.

← Back to Chapter 13