Numerical Lab: Integrate, Differentiate, Solve ODEs
Three classic numerical methods, each laying out its working the way a lecturer would on the board. Integrate approximates the area under a curve with Simpson's rule, which fits a parabola through each pair of strips; its table lists every point, its weight and the sum. Differentiate estimates the slope at a point from the function's values just either side of it, a central difference, and the second derivative from the same three values. Solve an ODE steps y′ = f(x, y) forward with the classical fourth-order Runge–Kutta method and tabulates the four slopes of every step.
Each answer is checked by running again with half the step: the change between the two runs estimates the error, and combining the pair (Richardson extrapolation) is usually closer still.
The working
Step h = (b − a) / n = (3.141593 − 0) / 8 = 0.3926991. Simpson's rule fits a parabola through each pair of strips, which weights the points 1, 4, 2, 4, …, 2, 4, 1:
| i | xᵢ | f(xᵢ) | weight | weight × f(xᵢ) |
|---|---|---|---|---|
| 0 | 0 | 0 | 1 | 0 |
| 1 | 0.3926991 | 0.3826834 | 4 | 1.530734 |
| 2 | 0.7853982 | 0.7071068 | 2 | 1.414214 |
| 3 | 1.178097 | 0.9238795 | 4 | 3.695518 |
| 4 | 1.570796 | 1 | 2 | 2 |
| 5 | 1.963495 | 0.9238795 | 4 | 3.695518 |
| 6 | 2.356194 | 0.7071068 | 2 | 1.414214 |
| 7 | 2.748894 | 0.3826834 | 4 | 1.530734 |
| 8 | 3.141593 | 1.224647e-16 | 1 | 1.224647e-16 |
Sum = 15.28093, and h / 3 × sum = 0.1308997 × 15.28093 = 2.000269.
Check: with n = 16 the rule gives 2.000017. Its error shrinks about 16× each time the step halves, so that answer is off by about |2.000017 − 2.000269| / 15 = 0.000016839, and 2.000017 + (2.000017 − 2.000269) / 15 = 2 is closer still (Richardson extrapolation).
Type a function of x with + − * / ^ and brackets; 2x means 2 × x. Also pi, e and sin cos tan asin acos atan sinh cosh tanh sqrt abs exp ln log floor ceil (angles in radians; log is base 10, ln natural). Number fields take pi/2 too.
The working
f(x₀ + h) = f(1.01) = 2.745601
f(x₀) = f(1) = 2.718282
f(x₀ − h) = f(0.99) = 2.691234
f′(x₀) ≈ (f(x₀ + h) − f(x₀ − h)) / 2h = (2.745601 − 2.691234) / 0.02 = 2.718327
f″(x₀) ≈ (f(x₀ + h) − 2 f(x₀) + f(x₀ − h)) / h² = 2.718304
Check: with h / 2 the slope comes out 2.718293. The central difference's error shrinks about 4× when h halves, so (4 × 2.718293 − 2.718327) / 3 = 2.718282 is usually far closer (Richardson extrapolation). A very small h loses digits to rounding instead; around 0.001 to 0.01 suits most smooth functions.
Type a function of x with + − * / ^ and brackets; 2x means 2 × x. Also pi, e and sin cos tan asin acos atan sinh cosh tanh sqrt abs exp ln log floor ceil (angles in radians; log is base 10, ln natural). Number fields take pi/2 too.
The working
Step h = (x − x₀) / n = 0.1. Each step samples the slope four times, k₁ = f(x, y), k₂ = f(x + h/2, y + h k₁/2), k₃ = f(x + h/2, y + h k₂/2) and k₄ = f(x + h, y + h k₃), then moves on with y + h (k₁ + 2k₂ + 2k₃ + k₄) / 6:
| step | x | y | k₁ | k₂ | k₃ | k₄ |
|---|---|---|---|---|---|---|
| 0 | 0 | 1 | 1 | 1.05 | 1.0525 | 1.10525 |
| 1 | 0.1 | 1.105171 | 1.105171 | 1.160429 | 1.163192 | 1.22149 |
| 2 | 0.2 | 1.221403 | 1.221403 | 1.282473 | 1.285526 | 1.349955 |
| 3 | 0.3 | 1.349858 | 1.349858 | 1.417351 | 1.420726 | 1.491931 |
| 4 | 0.4 | 1.491824 | 1.491824 | 1.566415 | 1.570145 | 1.648839 |
| 5 | 0.5 | 1.648721 | 1.648721 | 1.731157 | 1.735278 | 1.822248 |
| 6 | 0.6 | 1.822118 | 1.822118 | 1.913224 | 1.917779 | 2.013896 |
| 7 | 0.7 | 2.013752 | 2.013752 | 2.114439 | 2.119474 | 2.225699 |
| 8 | 0.8 | 2.22554 | 2.22554 | 2.336817 | 2.34238 | 2.459778 |
| 9 | 0.9 | 2.459601 | 2.459601 | 2.582581 | 2.58873 | 2.718474 |
| 10 | 1 | 2.71828 |
Check: with 20 steps y(1) comes out 2.718282. The method's error shrinks about 16× each time the step halves, so that answer is off by about |2.718282 − 2.71828| / 15 = 1.299e-7, and 2.718282 + (2.718282 − 2.71828) / 15 = 2.718282 is closer still.
Type a function of x and y (write x*y, not xy) with + − * / ^ and brackets; 2x means 2 × x. Also pi, e and sin cos tan asin acos atan sinh cosh tanh sqrt abs exp ln log floor ceil (angles in radians; log is base 10, ln natural). Number fields take pi/2 too.
The other accessories