Scientific computing, from the ground up

From real‑world mathematical models to engineering solutions

The power of scientific computing and simulation.

Harish Jayaraj
Scientific Computing and Software Development Engineer

Every one of these was tested before it existed

  • A bridge deck, in wind that had not blown yet
  • A processor, at a temperature it had never reached
  • A wing, in air it had never flown through
  • A car, in a crash that never happened

Not one of them was built first and checked afterwards. They were run inside a computer until somebody was willing to sign for them.

Not one part of this was measured

Simulated surface currents across the North Atlantic, drawn as dense red and blue swirls, with the Gulf Stream leaving the American coast as one strong ribbon.
North Atlantic surface currents, from a global ocean model. NASA / JPL-Caltech

Somebody wrote down how sea water moves, and let a computer step it forward.
Not one of those eddies was drawn. They are what the arithmetic did.

The loop every one of them went through

real world model discretise solve decide refine, and go round again

Two of those five boxes are where the difficulty actually lives.

Box two

Model

Decide what to keep and what to throw away. Keep too little and the answer is wrong. Keep too much and nobody can solve it.

Box four

Solve

Almost nothing worth modelling has a formula. Getting a number out of an equation you cannot solve is a craft of its own.

The next ninety minutes walk that loop end to end, once, on a problem small enough to hold in your head.

Whatever you are modelling, it ends up here

m \mathbf{a} = \sum \mathbf{F} m \dfrac{d\mathbf{v}}{dt} = \sum \mathbf{F} \dfrac{d\mathbf{u}}{dt} = \mathbf{f}(\mathbf{u})

A state, and a rule for how fast it changes.

Four questions, four different fields

  • How hot does this chip get before it throttles?
  • Will this circuit settle, or ring forever?
  • How many people are infected by Friday?
  • Where will the satellite be in six hours?

Different departments. Different journals. Different people, who mostly do not talk to each other.

Write each one down and they look unrelated

chip temperature

$$mc\,\frac{dT}{dt} = P - hA\,(T - T_\infty)$$

rlc circuit

$$L\,\frac{d^{2}q}{dt^{2}} + R\,\frac{dq}{dt} + \frac{q}{C} = 0$$

epidemic

$$\frac{dS}{dt} = -\beta S I, \qquad \frac{dI}{dt} = \beta S I - \gamma I$$

orbit

$$\frac{d^{2}\mathbf{r}}{dt^{2}} = -GM\,\frac{\mathbf{r}}{\lvert\mathbf{r}\rvert^{3}}$$

The one move that makes all four the same

A second derivative is awkward to solve numerically. So stop tracking position on its own, and track position and velocity together as a single state.

\dfrac{d^{2}\mathbf{r}}{dt^{2}} = -GM\dfrac{\mathbf{r}}{\lvert\mathbf{r}\rvert^{3}} \dfrac{d}{dt}\begin{bmatrix}\mathbf{r}\\ \mathbf{v}\end{bmatrix} = \begin{bmatrix}\mathbf{v}\\ -GM\,\mathbf{r}/\lvert\mathbf{r}\rvert^{3}\end{bmatrix} \dfrac{d\mathbf{u}}{dt} = \mathbf{f}(\mathbf{u})

Newton's law for an orbit.
The second derivative is the problem: our methods take first ones.

Top row: the definition \(d\mathbf{r}/dt = \mathbf{v}\).
Bottom row: the original equation, because \(d\mathbf{v}/dt\) is \(d^{2}\mathbf{r}/dt^{2}\).

Nothing was approximated, only rewritten.
Name the state \(\mathbf{u}\) and the rule \(\mathbf{f}\).

All four models do exactly this. They land on the same line.

The shape of every model in this talk

$$\frac{d\mathbf{u}}{dt} = \mathbf{f}(\mathbf{u})$$

u, the state

Every number you need to know where the system is right now. Temperature. Charge and current. Position and velocity.

f, the rule

Given the state, how fast is it changing. This is the only place your physics appears.

Which matters more than it sounds

the solver, written once time stepping step size control linear algebra error estimates event detection output and plots f(u) your physics

Swap the red box and you have changed field. Everything around it is written once, tested once, and reused by everybody.

One solver simulates a ball,
an amplifier and a chemical plant.

Learn the machinery once and it pays out for the rest of your career, in whichever field you end up in.

One problem, small enough to hold in your head

u₀ θ R
$$R = \frac{u_0^{2}\sin 2\theta}{g}$$

You have had this since school.
Maximum range at 45 degrees.

The ball, marked at equal time intervals

Horizontal speed never changes.

same Δt between marks, same Δx every time

Do you believe that about a real cricket ball?

What the formula quietly left out

AIR

Air resistance, in one term

The wind tunnel gives you a size. Common sense gives you a direction:
drag opposes motion. Multiply them together.

\lvert\mathbf{F}_D\rvert = \tfrac{1}{2}\rho C_D A \lvert\mathbf{v}\rvert^{2} \mathbf{F}_D = \tfrac{1}{2}\rho C_D A \lvert\mathbf{v}\rvert^{2} \left(-\dfrac{\mathbf{v}}{\lvert\mathbf{v}\rvert}\right) \mathbf{F}_D = -\tfrac{1}{2}\rho C_D A \lvert\mathbf{v}\rvert\mathbf{v}

One power of the speed cancels. It is \(\lvert\mathbf{v}\rvert\mathbf{v}\), a number times a vector, not a vector squared.

Put it in the sum, divide by the mass, group the constants

m\dfrac{d\mathbf{v}}{dt} = -mg\hat{\jmath} - \tfrac{1}{2}\rho C_D A\lvert\mathbf{v}\rvert\mathbf{v} \dfrac{d\mathbf{v}}{dt} = -g\hat{\jmath} - \dfrac{\rho C_D A}{2m}\lvert\mathbf{v}\rvert\mathbf{v} \dfrac{d\mathbf{v}}{dt} = -g\hat{\jmath} - k\lvert\mathbf{v}\rvert\mathbf{v} \dfrac{d\mathbf{u}}{dt} = \mathbf{f}(\mathbf{u})

Gravity's mass cancels.
A feather and a coin fall alike in vacuum.

Drag's mass does not.
It survives inside \(k = \rho C_D A / 2m\).

Now write it out in components

$$\frac{dv_x}{dt} = -k\,\underbrace{\sqrt{v_x^{2}+v_y^{2}}}_{\text{holds } v_y}\,v_x$$ $$\frac{dv_y}{dt} = -g - k\,\underbrace{\sqrt{v_x^{2}+v_y^{2}}}_{\text{holds } v_x}\,v_y$$

The horizontal equation now contains the vertical speed, and the vertical equation contains the horizontal one. They cannot be solved separately any more.

And there is no closed form.

There is no formula. So stop looking for one.

Give up knowing it
everywhere.

Settle for knowing it at \(t_0,\; t_1,\; t_2,\; \ldots\) separated by a step \(h\).

Exactness traded for existence. Every method in this act is that trade.

You already met this in your first calculus course

\mathbf{u}(t+h) = \mathbf{u}(t) + h\,\mathbf{u}'(t) + \dfrac{h^{2}}{2}\mathbf{u}''(\xi) \mathbf{u}(t+h) \approx \mathbf{u}(t) + h\,\mathbf{u}'(t) \mathbf{u}_{n+1} = \mathbf{u}_{n} + h\,\mathbf{f}(\mathbf{u}_{n})

Taylor's theorem, exactly true.
The last term is everything you are about to throw away.

Throw it away.
You have just made your first approximation, and you know its size.

And \(\mathbf{u}'\) is exactly what \(\mathbf{f}\) hands you.
That is Euler's method.

One step, drawn honestly

start true path Euler follows the tangent local error

The tangent is right where you are standing, and increasingly wrong everywhere else.

Drop one term of Taylor, and one step costs you \(O(h^{2})\).

Take \(T/h\) of them and the error accumulates to \(O(h)\). Halve the step, halve the error.

That is a bad exchange rate. Ten times the accuracy costs ten times the work.

Seven Euler steps at \(h = 0.6\,\mathrm{s}\), against the truth

Euler apex 23.6 m true apex 18.4 m

Every step fires off along the tangent and ignores the curve bending beneath it,
so the errors all lean the same way.

The defect, once you have seen the picture

Euler trusts the slope it measured at the start of the step.

By the time the step ends, that slope is stale.

So measure it four times instead of once.

Each probe decides where the next one is taken

tₙ tₙ + h/2 tₙ + h true solution k₁ k₂ k₃ k₄ uₙ₊₁ off by O(h⁵)

Two probes overshoot, two undershoot. The step takes their weighted average, and the four errors very nearly cancel.

RK4

$$\begin{aligned} \mathbf{k}_1 &= \mathbf{f}(\mathbf{u}_n) \\[2pt] \mathbf{k}_2 &= \mathbf{f}\!\left(\mathbf{u}_n + \tfrac{h}{2}\mathbf{k}_1\right) \\[2pt] \mathbf{k}_3 &= \mathbf{f}\!\left(\mathbf{u}_n + \tfrac{h}{2}\mathbf{k}_2\right) \\[2pt] \mathbf{k}_4 &= \mathbf{f}\!\left(\mathbf{u}_n + h\,\mathbf{k}_3\right) \end{aligned}$$
$$\mathbf{u}_{n+1} = \mathbf{u}_n + \frac{h}{6}\bigl(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4\bigr)$$

Four evaluations of \(\mathbf{f}\) per step instead of one.

Global error \(O(h^{4})\). Halve the step and the error falls by sixteen.

Four times the work buys thousands of times the accuracy. In computing, the algorithm is worth more than the hardware.

Do not take that on trust. Measure it.

h (s) Euler error ratio RK4 error ratio
0.16002.09 m—7.0×10⁻⁶ m—
0.08001.04 m2.014.2×10⁻⁷ m16.6
0.04000.52 m2.002.6×10⁻⁸ m16.3
0.02000.26 m2.001.6×10⁻⁹ m16.1
0.01000.13 m2.001.0×10⁻¹⁰ m15.9

2.00, row after row. 16, row after row. Not claimed, observed.

One more problem, and it is not integration

ground last point above first point below interpolated landing 0.06 m out 1.44 m if you stop at the last point

The solver marches in fixed steps. It has no idea the ground is there, so it steps straight past it.

Finding the landing is root finding: solve \(y(t^{*}) = 0\) between the two points that straddle it.

$$\tau = \frac{y_n}{y_n - y_{n+1}}, \qquad R \approx x_n + \tau\,(x_{n+1}-x_n)$$

Worth the care. The integrator was accurate to \(10^{-10}\) m. Stopping at the last point throws away ten orders of magnitude in the final step.

Everything in the last act, as code

Four functions. About thirty lines.

  • f   the physics, and nothing else
  • rk4_step   one step, the four probes
  • fly   the loop, and where it lands
  • sweep   ask it the question

Nothing here is a library call you have to take on faith.

01   the physics

import numpy as np g, rho = 9.81, 1.2 def f(u, k): # u = [x, y, vx, vy] x, y, vx, vy = u speed = np.hypot(vx, vy) return np.array([vx, vy, -k * speed * vx, -g - k * speed * vy])

The first two rows are the definition of velocity. The last two are Newton's law with drag.

This is the only file you rewrite when you change problem. Everything after this slide is the same for a ball, an orbit or a circuit.

02   one step

def rk4_step(u, h, k): k1 = f(u, k) k2 = f(u + h*k1/2, k) k3 = f(u + h*k2/2, k) k4 = f(u + h*k3, k) return u + h/6 * (k1 + 2*k2 + 2*k3 + k4)
[0] [1] [2] [3] return = [ x , y , vx , vy ]

The same four numbers back, one step later. Index 1 is y, which matters on the next slide.

The four probes from the diagram, in the order you watched them being taken.

Each line feeds the next. k2 is measured at the point k1 pointed to.

Six lines. That is the whole of the method that carries spacecraft across the solar system.

03   the loop, and where it lands

def fly(speed, angle, k, h=0.001): th = np.radians(angle) u = np.array([0.0, 0.0, speed*np.cos(th), speed*np.sin(th)]) while True: nxt = rk4_step(u, h, k) if nxt[1] < 0: tau = u[1] / (u[1] - nxt[1]) return u[0] + tau * (nxt[0] - u[0]) u = nxt

Degrees in, radians out.

The starting state: origin, speed split into vx and vy.

Work out a candidate step, but do not commit to it.

nxt[1] is \(y\). Negative means the step overshot the ground.

How far along that segment \(y\) reaches zero.

The same fraction applied to \(x\). That is the range.

[0] [1] [2] [3] x y vx vy u = [ 61.136, 1.940, 12.723, -16.540] nxt = [ 64.258, -2.421, 12.251, -18.332]

The two steps that straddle the ground. nxt[1] has gone negative, so the landing is somewhere between them.

04   ask it the question

m, d, Cd = 0.156, 0.072, 0.45 # a cricket ball A = np.pi * d**2 / 4 k = rho * Cd * A / (2 * m) angles = np.arange(10, 86, 1.0) ranges = [fly(30.0, a, k) for a in angles] best = angles[np.argmax(ranges)] print(best) # 42.0

Seventy six complete flight simulations, to answer one question with one unknown.

And the answer is 42 degrees, not 45.

Hold that thought. Twenty unknowns instead of one, and this loop is why supercomputers exist.

Before you believe any of it

# with no air, the solver must reproduce # the formula they gave you in school sim = fly(30.0, 45.0, k=0.0) exact = 30.0**2 * np.sin(np.radians(90)) / 9.81 print(sim, exact) # 91.743 91.743

Turn the air off and the hard problem
collapses into one you already know.

If those two numbers disagree, something is broken, and every result after it is worthless.

This is the habit worth stealing. Find the case where you already know the answer, and check it first.

Now run it. Everything here is computed as you drag.

  • Cricket ball
  • Football
  • Table-tennis ball
  • Shuttlecock
range, with air— range, no air— lost to the air— drag length 1/k—

Every point on this curve is a full flight simulation

  • Cricket ball
  • Football
  • Table-tennis ball
  • Shuttlecock
best angle, with air— range there— best angle, no air— range the air cost—

Push the speed up and watch the red peak walk away from 45.

So what is the answer?

The best angle is not 45. And it does not stay put.

cricket ball · 156 g · 72 mm · Cₔ = 0.45

launch speed best angle range there range at 45° what 45° costs
20 m/s43°33.6 m33.5 m0.04 m
30 m/s42°62.9 m62.6 m0.27 m
50 m/s39°118.4 m116.6 m1.78 m
92 m/s35°202.5 m194.3 m8.16 m
120 m/s32°241.7 m228.6 m13.09 m

At throwing speeds, aiming at 45 costs you four centimetres.
At bowling speeds and above, it costs you metres.

Why the peak moved at all

Drag charges rent on every second you spend in the air.

A steep launch buys height, and height buys time,
and time is now something you pay for.

So the winning shot is flatter and quicker, and gets flatter the harder you throw.

Worth noticing. The peak is flat. An optimum you do not have to hit precisely is a good one, and engineers actively look for that.

Count what that answer actually cost

  • 1 unknown, the launch angle
  • 76 complete flight simulations
  • ~4000 RK4 steps in each one
  • ~1.2 million evaluations of f

All of that, and it finished before you noticed. That is the good case.

Now make it a real design problem

A wing has twenty numbers you are allowed to change, not one.

Sweeping one unknown took 76 runs.
Sweeping twenty the same way takes more runs than there are atoms in the room.

So you stop sweeping, and start using the gradient.

And our \(\mathbf{f}\) was four numbers

a thrown ball x y vx vy 4 numbers → air over a car hundreds of millions of cells, several unknowns each

The loop does not change. Only the size of \(\mathbf{u}\) changes, and with it the machine you need.

This is the air a single wing throws away

A light aircraft flying low through coloured smoke, which is being wound into one enormous tight spiral trailing from its wingtip.
Wingtip vortex, made visible with coloured smoke. NASA Langley

Every wing does this, all the time, and it is the reason aircraft have to be spaced out on approach.
Predicting that spiral is the same job as predicting the ball, with a great many more unknowns.

Same pipeline, different \(\mathbf{f}\)

aerodynamics

Every aircraft and every car shape is flown thousands of times inside a computer before one is built.

crash and structures

A car is crashed a few hundred times in simulation for every one crashed on a test track.

circuits

No chip is fabricated before it has been solved, node by node, step by step, in exactly this way.

weather and climate

The forecast you read this morning is the same loop, run over the atmosphere, several times a day.

None of this is new

Artillery firing tables are exactly this computation.

Angle against range, corrected for muzzle speed, air density and wind.
Produced before electronic computers existed.

One trajectory, worked by hand, took about a day.
A single table needed thousands of them.

The job title was "computer". It described a person. ENIAC, the first general purpose electronic computer, was built to do this arithmetic faster.

What the work is done with

  • Python and NumPy for thinking, prototyping and plotting
  • C++, Fortran, Julia when the loop has to be fast
  • MPI and CUDA when one machine is no longer enough
  • OpenFOAM, deal.II, PETSc when somebody has already written the hard part

You used the first one today, and hand-wrote the part the last one usually hides.

What the job actually looks like

  • Decide what physics to keep, and defend the choice
  • Turn it into something a machine can step forward
  • Prove the answer is converged, not merely produced
  • Make it fast enough to be useful to somebody waiting
  • Explain the result to an engineer who will sign for it

You did the first three today, on a cricket ball.

The formula you memorised is a special case of an equation nobody can solve.

Scientific computing is how engineers work in the general case,
which is every case that has ever been built.

What you did in the last ninety minutes

  • Took a question with no analytical answer
  • Built a model from one law, and named what you threw away
  • Built a solver, and measured its accuracy instead of trusting it
  • Found an answer nobody in the room could have guessed

That is the whole working method. The cricket ball was incidental.

Take home

A different system. The same method.

A pendulum, length \(L = 1\,\mathrm{m}\), released from rest at an angle \(\theta_0\). You were given this in school:

$$T = 2\pi\sqrt{L/g}$$

Look at what is missing from the right hand side. The amplitude. That formula claims a big swing and a small swing take exactly the same time.

The honest model has no such claim in it:

$$\frac{d^{2}\theta}{dt^{2}} = -\frac{g}{L}\sin\theta$$

The school formula comes from replacing \(\sin\theta\) with \(\theta\), which is only true for small angles. Put \(\sin\theta\) back and the closed form is gone.

The question. At what amplitude is the true period one per cent longer than the formula says?

Four things, one page

  • One check. At \(\theta_0 = 1^{\circ}\) your period must match \(2\pi\sqrt{L/g}\) to three decimals. Run this first.
  • One plot. Period against amplitude, from \(1^{\circ}\) to \(120^{\circ}\), with the school formula drawn as a flat line.
  • One number. The amplitude where your curve first sits one per cent above that line.
  • One paragraph. The period gets longer, not shorter. Explain why, using \(\sin\theta < \theta\).

If the check fails, stop. Nothing else means anything until it passes.

If you want more. Push the amplitude towards \(180^{\circ}\). The period grows without limit. Say what the pendulum is physically doing, and what that does to your step size.

Two things you already know how to do

import numpy as np g, L = 9.81, 1.0 def f(u): # u = [theta, omega] return np.array([u[1], -(g/L)*np.sin(u[0])]) def rk4_step(u, h): # slide 29, minus the k ... def quarter(theta0, h=1e-4): u = np.array([np.radians(theta0), 0.0]) t = 0.0 while True: nxt = rk4_step(u, h) if nxt[0] <= 0: # crossed the bottom tau = u[0] / (u[0] - nxt[0]) return t + tau*h u = nxt; t += h period = 4 * quarter(theta0)

The state is two numbers, angle and angular velocity. Exactly the reduction from the start of the talk.

Finding the period is event detection. Release from rest, integrate until \(\theta\) crosses zero, and interpolate the crossing. That is a quarter of a period, and it is the landing problem again.

Nothing new is required. Change four lines of what you already have and the pendulum falls out.

Send me what you build

Harish Jayaraj
Scientific Computing and Software Development Engineer

codeynamics@gmail.com
linkedin.com/in/harish-jayaraj-p

Questions.

Codeynamics

Slide index, click to jump, Esc to close

Keys

→ / Spacenext build or slide
←back
↓ / ↑skip a whole slide
Home / Endfirst / last
Ffull screen
Oslide index
?this panel