Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 11

Solving Nonlinear Equations

These exercises supplement the problems in the book. They emphasize method selection, misleading convergence, and the possibility of missing a solution. The projects use different physical models from those assigned in Chapter 11. For each calculation, plot the function on its stated domain, report the root and residual, and verify the physical meaning of the result. The parameters below are reproducible starting points; vary them after the initial calculation.

Chapter 11 codes Chapter 11 resources Chapter 11 updates & errata Chapter 11 overview

Additional exercises

  1. Root accuracy and the scale of a residual

    Solve fs(x) = s(x − 0.37) = 0 on [0,1] for s = 10−6, 1, 106. Use bisection with a final bracket width of 10−7 and compare the estimated position with the exact root x = 0.37. Record |fs(x)| for each scale.

    Then try a residual stopping rule of 10−8 on the same bracket. Does it guarantee comparable accuracy in x for all three scales? Explain which tolerance expresses a useful physical error after x and f have been assigned units.

  2. When false position barely moves

    Find the root of f(x) = x10 − 1 on [0,2]. Compare ordinary false position and bisection using a requested residual of 10−8 for false position and a final bracket width of 10−8 for bisection. Limit each method to 200 iterations. Record the current bracket, root estimate, and residual after 5, 20, 100, and 200 steps where applicable.

    The exact root is x = 1. Which endpoint stays fixed under false position? Compare with a library bracketed solver if available, taking care to compare the actual root errors rather than unlike stopping criteria. As a further experiment, implement the Illinois modification, which reduces the stored function value at an endpoint that remains fixed, and compare its progress with ordinary false position.

  3. Two crossings inside one scan cell

    Scan f(x) = (x − 0.3734)(x − 0.3917) on [0,1] with 10 and then 100 equal subintervals. Refine every detected sign-change bracket to width 10−8. The exact roots are 0.3734 and 0.3917.

    Explain why the coarse scan reports no sign change despite two roots. Plot the function around the two crossings and identify the smallest scan spacing needed to separate them reliably. Repeat with 1000 subintervals as a check.

  4. A Newton iteration that cycles

    Apply ordinary Newton iteration to f(x) = x3 − 2x + 2, with f′(x) = 3x2 − 2 and initial guess x0 = 0. Print the first ten iterates instead of waiting for the 200-step limit. Verify algebraically that the iterates alternate between 0 and 1. Then start at x0 = −2 and request |f(x)| ≤ 10−10.

    Independently bisect [−2,−1] to width 10−10. Compare both converged roots with x ≈ −1.769292354239. What does the first run show about the need to monitor convergence and set an iteration limit?

  5. A sign change across a pole

    Consider f(x) = tan(x) − 2 on [1.2,1.8]. Its endpoint signs differ, but the interval contains the pole at π/2 and no root. Try bisection with bracket-width tolerance 10−8, recording |f(x)| as the bracket shrinks, as well as the final position. Explain why the usual sign-change guarantee does not apply.

    Now solve on the continuous bracket [1,1.2]. The true root there is arctan(2). Exclude a margin of 0.01 on both sides of π/2 when plotting the two continuous pieces.

  6. Companion code versus a library routine

    Solve f(x) = e−x − x on [0,1] with the companion bisection code and a library routine such as MATLAB fzero, SciPy root_scalar, Boost.Math, or Netlib ZEROIN. Set a requested position tolerance of approximately 10−10 where the routine permits it. The reference root is x ≈ 0.567143290410. For an independent check, this value equals the principal Lambert W function W0(1).

    Compare root error, residual, status information, and function evaluations if reported. Check the library documentation for its precise stopping rule. Which implementation would you use in a larger calculation, and which makes the bracket updates easier to inspect?

Physics projects

The first three projects have one physical root for their default parameters. The final two have several roots and require care near a known solution or a singularity. Keep each variable within its physical domain and check the relevant balance or equilibrium equation after solving.

1. Temperature of a heated surface

Objective: find the steady temperature of a surface that loses heat by both convection and thermal radiation.

P = HA(T − Ta) + εσA(T4 − Ta4).

Here P is supplied power, A is area, H is the convection coefficient, T and Ta are absolute temperatures, ε is emissivity, and σ is the Stefan–Boltzmann constant. Use A = 0.020 m2, H = 10 W m−2 K−1, ε = 0.8, Ta = 295 K, P = 20 W, and σ = 5.670374419 × 10−8 W m−2 K−4. Find the physical root on [295,500] K using a bracket-width tolerance of 10−6 K. Then try Newton’s method from the convection-only estimate T0 = Ta + P/(HA) = 395 K. The derivative of the heat loss is HA + 4εσAT3. Request a residual below 10−8 W and compare the two results. Vary P and plot the equilibrium temperature.

Validation: substitution into the heat balance should reproduce P. The default gives T ≈ 356.270872 K. At P = 0, the root is exactly T = Ta; with ε = 0, it is T = Ta + P/(HA). The heat loss increases with T for T ≥ Ta, so the default physical root is unique.

Possible extension: vary emissivity and convection strength to identify the range where radiation carries most of the heat loss.

2. How deeply does a sphere float?

Objective: determine the submerged depth of a homogeneous sphere from the balance of buoyancy and weight.

πd2(R − d/3) = (ρs/ρf)4πR3/3,    0 ≤ d ≤ 2R.

Here d is the height of the submerged spherical cap measured from the bottom of the sphere, R its radius, and ρs and ρf are the sphere and fluid densities. Use R = 0.100 m, ρf = 1000 kg/m3, and ρs = 650 kg/m3. Find d in [0,0.200] m to a bracket width of 10−9 m. Plot d/R versus ρs/ρf between 0.1 and 0.9. In terms of y = d/R, the balance is y3 − 3y2 + 4ρs/ρf = 0. Compare all three polynomial roots with the physical range 0 ≤ y ≤ 2.

Validation: the default gives d/R ≈ 1.202779395. Substitute the result into the two volumes and compare the buoyant force with the sphere's weight. At density ratio 0.5, symmetry gives exactly d = R. The submerged volume grows monotonically over the allowed interval, so there is one physical root for each density ratio strictly between 0 and 1.

Possible extension: examine how changing the fluid density, as in fresh versus salt water, changes the submerged fraction for a sphere of fixed density.

3. Resonant frequency with frequency-dependent damping

Objective: locate the frequency of maximum response of a driven oscillator when the damping depends weakly on frequency.

Let q = ω/ω0. In this model, the response amplitude is proportional to 1/sqrt(D(q)), with

D(q) = (1 − q2)2 + 4g(q)2q2,   g(q) = g0(1 + βq2).

Choose g0 = 0.08 and β = 0.30. Find the nonzero solution of D′(q) = 0 on [0.8,1.1] with bracket-width tolerance 10−9. An explicit derivative is

D′(q) = −4q(1 − q2) + 8g(q)2q + 8g(q)g′(q)q2,   g′(q) = 2g0βq.

Plot the amplitude around the root and check that the root gives a maximum, not merely a stationary point. The endpoint q = 0 is another stationary point of D and must not be selected as the resonance. Since every term in D′(q) contains q, factor out this known root before solving for a nonzero resonance.

Validation: for β = 0, the nonzero resonance is exactly q = sqrt(1 − 2g02) when this is positive. For the default frequency-dependent damping, q ≈ 0.984413316. Check that the amplitude decreases on both sides and that direct maximization of the sampled response gives a compatible frequency.

Possible extension: vary β and g0 to see how the peak moves and when a nonzero resonance peak ceases to exist.

4. Spontaneous magnetization in a mean-field model

Objective: locate the possible magnetizations below the Curie temperature and distinguish mathematical roots from stable physical states.

m = tanh((m + h)/τ),   τ = T/Tc.

Here m is reduced magnetization, h is reduced applied field, and τ is reduced temperature. Start with τ = 0.9 and h = 0. The root m = 0 is known. Find the positive nonzero root on [0.1,1] to a bracket width of 10−9 and obtain its negative counterpart by symmetry. Try Newton’s method with m0 = 0 and then m0 = 1, using f′(m) = 1 − sech2((m+h)/τ)/τ and a residual tolerance of 10−10. Follow the positive root for 0.20 ≤ τ ≤ 0.99 and plot it against τ. At lower temperatures the root becomes indistinguishable from 1 in ordinary double precision.

Next set τ = 0.8 and h = −0.02. Scan [−1,1] with 1000 equal cells and refine each sign-change bracket to width 10−9. A stationary state is locally stable when τ/(1 − m2) − 1 > 0. Use this criterion and the sign of the applied field to label the stable, metastable, and unstable solutions.

Validation: the default positive root is m ≈ 0.525429513. In zero field, inversion gives τ = m/artanh(m); in particular, m = 0.5 corresponds to τ ≈ 0.910239227. Near τ = 1, compare with m ≈ sqrt(3(1 − τ)); at τ = 0.99, this gives 0.173205 versus the numerical root 0.172511. The three applied-field roots are approximately −0.739499930, 0.101398707, and 0.673074790. The negative state is stable, the positive outer state is metastable, and the middle state is unstable.

Possible extension: sweep h in both directions, following a locally stable root from its previous value, and plot the hysteresis loop. The spinodal field at τ = 0.8 has magnitude about 0.062244.

5. Collinear Earth–Moon Lagrange points

Objective: locate L1, L2, and L3 along the Earth–Moon line while distinguishing equilibrium points from sign changes across singularities.

In a rotating frame, take the Earth–Moon distance, total mass, and angular velocity as units. Earth is at x = −μ, the Moon at x = 1 − μ, and μ = 0.012150585. Equilibrium on the axis requires F(x) = 0, where

F(x) = x − (1 − μ)(x + μ)/|x + μ|3 − μ(x − 1 + μ)/|x − 1 + μ|3.

First scan [−2,2] with 4000 equal cells. It gives five sign changes, including two across the poles at the masses. Examine |F(x)| near these poles. Then bracket the three actual roots separately on [−2,−0.1], [0.1,0.9], and [1.05,2], using a final bracket width of 10−10. Convert distances to kilometers with 384400 km per distance unit.

Validation: the roots are L1 ≈ 0.836915129, L2 ≈ 1.155682163, and L3 ≈ −1.005062646. L1 lies about 58019 km from the Moon and L2 about 64515 km beyond it. The approximate Hill radius (μ/3)1/3 ≈ 0.159401 lies between the Moon-to-L1 and Moon-to-L2 distances. The approximation L3 ≈ −(1+5μ/12) differs by about 10−7 in these units. Check the residual only inside each continuous interval.

Possible extension: repeat with μ = 3.0 × 10−6 for the Sun–Earth system. Adjust the scan spacing near Earth to distinguish L1 and L2 from its pole.

← Back to Chapter 11