Solving Equations and Calculus by Computer: Quick Revision

A brief note on Solving Equations and Calculus by Computer. Read it before a test. Try each question at the end before you open its answer.

Roots of an equation

• A root of $f$ is a number $x$ with $f(x) = 0$. Many equations have no formula for the root, so the computer searches for it.
• If $f$ is continuous and changes sign between $a$ and $b$, a root lies between them. The pair is a bracket.
• Bisection halves the bracket at each step and keeps the half with the change of sign. It always works on a bracket, but each new correct digit takes about three steps.
• Newton's method follows the tangent down to the axis:
  • $x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}$
• Near a simple root, Newton doubles the number of correct digits at each step. It needs $f'(x)$, and it can fail at a flat tangent, in a cycle, or from a far start.
• brentq(f, a, b) from scipy.optimize solves on a bracket, safely and quickly.

Polynomials, fits and systems

• A polynomial is an array of its coefficients, from the highest power down: $x^3 - 4x$ is [1, 0, -4, 0].
• np.polyval(p, x) works out a value, np.roots(p) solves, and np.polyder(p) differentiates. Roots can be complex.
• np.polyfit(x, y, 1) gives the slope and then the intercept of the best straight line. Look at the residuals too.
• curve_fit(model, x, y, p0=...) fits any model; the variable is the first argument of model. The uncertainties are np.sqrt(np.diag(cov)).
• Linear equations are $Ax = b$. np.linalg.solve(A, b) solves them, and np.allclose(A @ x, b) checks the answer.
• @ is the matrix product; * works place by place. A zero determinant means no single solution.

Slopes and areas from numbers

• The forward difference has an error proportional to $h$:
  • $f'(x) \approx \dfrac{f(x + h) - f(x)}{h}$
• The central difference has an error proportional to $h^2$:
  • $f'(x) \approx \dfrac{f(x + h) - f(x - h)}{2h}$
• A tiny $h$ brings rounding errors. For the central difference, $h$ near $10^{-5}$ works well.
• An area is a sum of narrow pieces. Rectangles have an error proportional to the width; trapezoids, to the width squared.
• Distance is the area under a speed graph. Work is the area under a force graph.

Motion in small steps

• A differential equation gives a rate of change. Euler's method moves forward in small steps of $dt$:
  • $y_{\text{new}} = y + \text{rate} \times dt$
• Euler is first order: half the step gives half the error.
• With air drag $-kv$, a falling body comes closer and closer to the speed $g/k$.
• solve_ivp(rate, (t0, t1), start) chooses its own steps. rate(t, y) takes the time first and returns a list of rates.
• An event function with terminal = True stops the solution when something happens, such as a landing.

A worked example

Find the cube root of 20 with Newton's method, starting at 3.

1. Write it as a root. $f(x) = x^3 - 20$, so $f'(x) = 3x^2$.
2. Step 1. $x_1 = 3 - \dfrac{27 - 20}{27} = 2.740741$.
3. Step 2. $f(x_1) = 0.587512$ and $f'(x_1) = 22.53498$, so $x_2 = 2.714670$.
4. Step 3 gives 2.714418, whose cube is 20.0000005. The errors were about 0.29, 0.026, 0.00025 and then 0.00000002.

The number of correct digits roughly doubles at each step. Bisection would need about 25 steps for the same accuracy.

Practice with new numbers

A program to try

A coach reads a sprinter's speed from a video once every second. The program finds the distance from the area under the speed graph, and the acceleration from its slope.

Mistakes to avoid

• No change of sign. brentq on $x^2 + 1$ from 0 to 2 stops with ValueError: f(a) and f(b) must have different signs. Draw the graph to find a bracket first.
• Newton with no limit on the steps. It can cycle or run away. Stop after, say, 50 steps, and check that $f(x)$ is close to 0.
• Coefficients in the wrong order. np.roots wants the highest power first, and a 0 for every missing power.
• *` for a matrix product.** A * x works place by place. Use A @ x`.
• A singular matrix. np.linalg.solve stops with LinAlgError: Singular matrix when the equations do not fix one answer.
• The arguments of a rate function in the wrong order. solve_ivp calls rate(t, y), with the time first.
• Trusting a method with no test. Try a case that you can solve by hand first.

Check yourself

1. A continuous function has $f(1) = -2$ and $f(2) = 3$. What does this tell you, and which point does bisection try next?

2. Why does Newton's method fail if it reaches a point where $f'(x) = 0$?

3. Write the array of coefficients for $2x^3 - 5x + 7$.

4. You halve $h$ in a forward difference. What happens to its error? And to the error of a central difference?

5. The force on a stretched cord is measured at six extensions. Should you find the work with quad or with np.trapezoid?

6. Euler's method with $dt = 0.1$ s gives a final position that is 0.8 m wrong. About how wrong is it with $dt = 0.05$ s?