University of Wales Trinity Saint David
ACCB5017 - Games Architecture and Mechanics

Runge-Kutta 4 Integration

Overview

This lab builds a small SFML program, RK4Demo.cpp, that runs the same spring-mass oscillator twice side by side. The left copy is advanced with the semi-implicit Euler update from the Numerical Integration Fundamentals lab and the right copy with fourth-order Runge-Kutta (RK4). Both use the same deliberately coarse fixed timestep. A grey ring drawn from the exact analytic solution marks where each mass should be, and a live readout reports how far each one has strayed from it. You will build the window, the oscillator state and the Euler baseline, watch it drift away from the exact answer, then write RK4 and compare the two.

Learning Objectives

By the end of this lab you will be able to:

Concept

Semi-implicit Euler is cheap and well behaved, but every step treats the rate of change as if it stayed constant for the whole step. On a path that keeps curving, such as a mass swinging on a spring, each step lands a little off the true curve, and the small errors add up. Run at a coarse step size, the simulated oscillator slowly falls out of phase with the real motion.

Fourth-order Runge-Kutta (RK4) keeps the idea of advancing in fixed steps but looks at the curve more carefully before committing to one. It samples the derivative four times within a single step. k1 is the derivative at the start of the step. k2 is the derivative at the midpoint, found by moving the state half a step along k1. k3 is another midpoint sample, this time moving half a step along k2. k4 is the derivative at the end of the step, found by moving a whole step along k3. The step then advances by a weighted average of the four, with the two midpoint samples counting double. Because the error of a single step falls with roughly the fifth power of the step size, RK4 tracks a smooth curve far more closely than a single-slope update at the same step size. The price is four derivative evaluations per step instead of one.

The state is held in a small Oscillator struct of displacement and velocity. A derivative() function returns how fast each of those is changing, and both integrators are written as functions that take an Oscillator& and a step size dt. The exact answer to compare against is a cosine, computed with std::cos.

Tasks

Task 1 - Opening the window

Create a new Console App (C++) project in Visual Studio the same way as in the Environment Setup lab, naming it RK4Demo, and link SFML through vcpkg as before. Visual Studio generates a RK4Demo.cpp containing a Hello World program. Replace its whole content with the window below.

Success condition: build and run the project. A white 900 by 700 window titled Semi-Implicit Euler vs RK4 opens and closes when you click its close button.

RK4Demo.cpp
#include <SFML/Graphics.hpp>

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        window.clear(sf::Color::White);
        window.display();
    }
}

Task 2 - Describing the oscillator as a state and a derivative

A mass on a spring obeys Hooke's law, F = -k * x, so its acceleration is -(k / m) * x. That gives two first-order equations, one for how displacement changes and one for how velocity changes.

dxdt=v,dvdt=−km x\frac{dx}{dt} = v, \qquad \frac{dv}{dt} = -\frac{k}{m}\,x

Above main() in RK4Demo.cpp, add an Oscillator struct holding displacement and velocity, two constants springConstant and mass, and a derivative() function. Given a state, derivative() returns a second Oscillator whose fields hold the two rates of change from the equations above.

The project still builds, but nothing new appears on screen yet. Nothing calls derivative() until a later task.

RK4Demo.cpp
#include <SFML/Graphics.hpp>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        window.clear(sf::Color::White);
        window.display();
    }
}

Task 3 - Drawing the mass on its spring

Extend main() so the state is visible. Add an Oscillator called euler starting at initialDisplacement with zero velocity, an anchor point at the top of the window, and a restLength that is the distance from the anchor to the mass when the spring is at rest. Inside the loop, position a sf::CircleShape at restLength + euler.displacement below the anchor and draw a grey line from the anchor to it.

Success condition: build and run the project. An orange ball hangs 80 pixels below its rest position on a grey line from the top of the window. It does not move yet.

RK4Demo.cpp
#include <SFML/Graphics.hpp>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    const float initialDisplacement = 80.0f;
    Oscillator euler{ initialDisplacement, 0.0f };

    const sf::Vector2f anchor(450.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape ball(14.0f);
    ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    ball.setFillColor(sf::Color(200, 120, 30));

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        window.clear(sf::Color::White);

        sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + euler.displacement);

        sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
        rod[0].position = anchor;
        rod[0].color = sf::Color(120, 120, 120);
        rod[1].position = ballPosition;
        rod[1].color = sf::Color(120, 120, 120);

        ball.setPosition(ballPosition);
        window.draw(rod);
        window.draw(ball);

        window.display();
    }
}

Task 4 - A fixed-timestep loop and an empty Euler step

Add a stepSemiImplicitEuler() function and call it from a fixed-timestep loop built the same way as in the Fixed-Timestep Game Loop lab. Real elapsed time is added to an accumulator, and a while loop spends it in whole fixedStep pieces. The step is set to 0.1f on purpose. It is coarse enough that integration error becomes easy to see.

Leave the body of stepSemiImplicitEuler() as the two TODO comments for now. The project builds and the ball still hangs still, because the step does not change the state yet.

RK4Demo.cpp
#include <SFML/Graphics.hpp>
#include <algorithm>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

void stepSemiImplicitEuler(Oscillator& s, float dt)
{
    // TODO: add the acceleration, derivative(s).velocity, scaled by dt to s.velocity
    // TODO: add the updated s.velocity scaled by dt to s.displacement
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    const float initialDisplacement = 80.0f;
    Oscillator euler{ initialDisplacement, 0.0f };

    // a deliberately coarse fixed step of 10 updates per second
    const float fixedStep = 0.1f;
    float accumulator = 0.0f;
    sf::Clock clock;

    const sf::Vector2f anchor(450.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape ball(14.0f);
    ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    ball.setFillColor(sf::Color(200, 120, 30));

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        accumulator += std::min(clock.restart().asSeconds(), 0.25f);

        while (accumulator >= fixedStep)
        {
            stepSemiImplicitEuler(euler, fixedStep);
            accumulator -= fixedStep;
        }

        window.clear(sf::Color::White);

        sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + euler.displacement);

        sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
        rod[0].position = anchor;
        rod[0].color = sf::Color(120, 120, 120);
        rod[1].position = ballPosition;
        rod[1].color = sf::Color(120, 120, 120);

        ball.setPosition(ballPosition);
        window.draw(rod);
        window.draw(ball);

        window.display();
    }
}

Task 5 - Completing the semi-implicit Euler step

Fill in the two TODO lines in stepSemiImplicitEuler() using the same update order as in the Numerical Integration Fundamentals lab.

  1. Add the acceleration, derivative(s).velocity, scaled by dt to s.velocity.
  2. Add the just-updated s.velocity scaled by dt to s.displacement.

Success condition: build and run the project. The orange ball swings up and down about its rest position roughly once per second, moving in ten visible jumps per second because the state only updates at the fixed step.

RK4Demo.cpp
#include <SFML/Graphics.hpp>
#include <algorithm>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

void stepSemiImplicitEuler(Oscillator& s, float dt)
{
    s.velocity += derivative(s).velocity * dt;
    s.displacement += s.velocity * dt;
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    const float initialDisplacement = 80.0f;
    Oscillator euler{ initialDisplacement, 0.0f };

    // a deliberately coarse fixed step of 10 updates per second
    const float fixedStep = 0.1f;
    float accumulator = 0.0f;
    sf::Clock clock;

    const sf::Vector2f anchor(450.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape ball(14.0f);
    ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    ball.setFillColor(sf::Color(200, 120, 30));

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        accumulator += std::min(clock.restart().asSeconds(), 0.25f);

        while (accumulator >= fixedStep)
        {
            stepSemiImplicitEuler(euler, fixedStep);
            accumulator -= fixedStep;
        }

        window.clear(sf::Color::White);

        sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + euler.displacement);

        sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
        rod[0].position = anchor;
        rod[0].color = sf::Color(120, 120, 120);
        rod[1].position = ballPosition;
        rod[1].color = sf::Color(120, 120, 120);

        ball.setPosition(ballPosition);
        window.draw(rod);
        window.draw(ball);

        window.display();
    }
}

Task 6 - Drawing the exact solution as a ring

An undamped mass on a spring has an exact solution that needs no stepping at all. Starting from rest at displacement x0, its position at time t is a cosine with angular frequency omega.

x(t)=x0cos⁡(ωt),ω=kmx(t) = x_0 \cos(\omega t), \qquad \omega = \sqrt{\frac{k}{m}}
  1. Add #include <cmath> and an omega constant computed with std::sqrt.
  2. Add an elapsed variable and increase it by fixedStep inside the while (accumulator >= fixedStep) loop, so it counts simulated time exactly.
  3. Each frame, compute exact as initialDisplacement * std::cos(omega * elapsed) and place a transparent, grey-outlined exactRing at that displacement below the anchor.

Success condition: build and run the project. A grey ring appears exactly around the orange ball at the start. Both move, and the ring shows where a perfect simulation would have the mass.

RK4Demo.cpp
#include <SFML/Graphics.hpp>
#include <algorithm>
#include <cmath>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

void stepSemiImplicitEuler(Oscillator& s, float dt)
{
    s.velocity += derivative(s).velocity * dt;
    s.displacement += s.velocity * dt;
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    const float initialDisplacement = 80.0f;
    const float omega = std::sqrt(springConstant / mass);
    Oscillator euler{ initialDisplacement, 0.0f };

    // a deliberately coarse fixed step of 10 updates per second
    const float fixedStep = 0.1f;
    float accumulator = 0.0f;
    float elapsed = 0.0f;
    sf::Clock clock;

    const sf::Vector2f anchor(450.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape ball(14.0f);
    ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    ball.setFillColor(sf::Color(200, 120, 30));

    // an outline where the mass should be according to the exact solution
    sf::CircleShape exactRing(14.0f);
    exactRing.setOrigin(sf::Vector2f(14.0f, 14.0f));
    exactRing.setFillColor(sf::Color::Transparent);
    exactRing.setOutlineColor(sf::Color(90, 90, 90));
    exactRing.setOutlineThickness(2.0f);

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        accumulator += std::min(clock.restart().asSeconds(), 0.25f);

        while (accumulator >= fixedStep)
        {
            stepSemiImplicitEuler(euler, fixedStep);
            elapsed += fixedStep;
            accumulator -= fixedStep;
        }

        float exact = initialDisplacement * std::cos(omega * elapsed);

        window.clear(sf::Color::White);

        sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + euler.displacement);

        sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
        rod[0].position = anchor;
        rod[0].color = sf::Color(120, 120, 120);
        rod[1].position = ballPosition;
        rod[1].color = sf::Color(120, 120, 120);

        ball.setPosition(ballPosition);
        exactRing.setPosition(anchor + sf::Vector2f(0.0f, restLength + exact));
        window.draw(rod);
        window.draw(exactRing);
        window.draw(ball);

        window.display();
    }
}

Task 7 - Watching Euler drift

Run the program and watch the ball against the ring for about ten seconds. Then try one experiment.

  1. Watch the ring and the ball for the first few swings, then keep watching. Note when they stop overlapping.
  2. Close the program, change fixedStep from 0.1f to 0.05f and run again. Watch for the same amount of time.
  3. Set fixedStep back to 0.1f before moving on.

Success condition: with fixedStep at 0.1f, the ball starts inside the ring but by around five seconds it is often well outside it, sometimes at the opposite end of its swing. The oscillation has kept its size but is now out of phase with the exact motion. With fixedStep at 0.05f the drift is much slower and the ball stays close to the ring for longer. That gap is the integration error you will now try to remove without shrinking the step.

Task 8 - How one RK4 step samples the curve

Writing y for the whole state (displacement and velocity together) and f for the derivative() function, one RK4 step of size dt takes four samples and blends them.

k1=f(yn)k2=f ⁣(yn+Δt2 k1)k3=f ⁣(yn+Δt2 k2)k4=f ⁣(yn+Δt k3)\begin{aligned} k_1 &= f(y_n) \\ k_2 &= f\!\left(y_n + \tfrac{\Delta t}{2}\,k_1\right) \\ k_3 &= f\!\left(y_n + \tfrac{\Delta t}{2}\,k_2\right) \\ k_4 &= f\!\left(y_n + \Delta t\,k_3\right) \end{aligned}
yn+1=yn+Δt6(k1+2k2+2k3+k4)y_{n+1} = y_n + \frac{\Delta t}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)

The diagram below shows one such step on the oscillator's displacement curve, drawn with a much larger step than the program uses so the curvature is easy to see. Each coloured segment is the slope of one sample, placed where that sample was taken. The red hollow point is where k1 alone would have carried the state, using a single slope for the whole step. The navy point is where the weighted average of all four lands, and the gold ring is the exact answer.

Reading the diagram is enough for this task. There is nothing to build yet.

Task 9 - Setting up the RK4 step

Add two functions above main(). addScaled() returns a copy of a state moved along a derivative by some scale, which is the y_n + (dt / 2) * k part of each sample. It is given below in full. stepRK4() starts with k1 computed for you and four TODO comments for the remaining work, which you will write in the next task.

The project still builds and behaves exactly as in Task 6. Nothing calls stepRK4() yet.

RK4Demo.cpp
#include <SFML/Graphics.hpp>
#include <algorithm>
#include <cmath>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

void stepSemiImplicitEuler(Oscillator& s, float dt)
{
    s.velocity += derivative(s).velocity * dt;
    s.displacement += s.velocity * dt;
}

Oscillator addScaled(const Oscillator& s, const Oscillator& d, float scale)
{
    return { s.displacement + d.displacement * scale, s.velocity + d.velocity * scale };
}

void stepRK4(Oscillator& s, float dt)
{
    Oscillator k1 = derivative(s);
    // TODO: k2 is the derivative at s moved half a step along k1
    // TODO: k3 is the derivative at s moved half a step along k2
    // TODO: k4 is the derivative at s moved a whole step along k3
    // TODO: add dt / 6 times (k1 + 2 * k2 + 2 * k3 + k4) to s.displacement and to s.velocity
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    const float initialDisplacement = 80.0f;
    const float omega = std::sqrt(springConstant / mass);
    Oscillator euler{ initialDisplacement, 0.0f };

    // a deliberately coarse fixed step of 10 updates per second
    const float fixedStep = 0.1f;
    float accumulator = 0.0f;
    float elapsed = 0.0f;
    sf::Clock clock;

    const sf::Vector2f anchor(450.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape ball(14.0f);
    ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    ball.setFillColor(sf::Color(200, 120, 30));

    // an outline where the mass should be according to the exact solution
    sf::CircleShape exactRing(14.0f);
    exactRing.setOrigin(sf::Vector2f(14.0f, 14.0f));
    exactRing.setFillColor(sf::Color::Transparent);
    exactRing.setOutlineColor(sf::Color(90, 90, 90));
    exactRing.setOutlineThickness(2.0f);

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        accumulator += std::min(clock.restart().asSeconds(), 0.25f);

        while (accumulator >= fixedStep)
        {
            stepSemiImplicitEuler(euler, fixedStep);
            elapsed += fixedStep;
            accumulator -= fixedStep;
        }

        float exact = initialDisplacement * std::cos(omega * elapsed);

        window.clear(sf::Color::White);

        sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + euler.displacement);

        sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
        rod[0].position = anchor;
        rod[0].color = sf::Color(120, 120, 120);
        rod[1].position = ballPosition;
        rod[1].color = sf::Color(120, 120, 120);

        ball.setPosition(ballPosition);
        exactRing.setPosition(anchor + sf::Vector2f(0.0f, restLength + exact));
        window.draw(rod);
        window.draw(exactRing);
        window.draw(ball);

        window.display();
    }
}

Task 10 - Completing stepRK4

Fill in the four TODO comments in stepRK4().

  1. Compute k2 as derivative(addScaled(s, k1, dt * 0.5f)).
  2. Compute k3 the same way but moving half a step along k2.
  3. Compute k4 by moving a whole dt along k3.
  4. Add (dt / 6.0f) times k1 + 2 * k2 + 2 * k3 + k4 to s.displacement, and the matching sum of the .velocity fields to s.velocity.

The finished function is shown below to check your work against. The project still builds and behaves as before, because nothing calls stepRK4() until the next task.

RK4Demo.cpp
#include <SFML/Graphics.hpp>
#include <algorithm>
#include <cmath>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

void stepSemiImplicitEuler(Oscillator& s, float dt)
{
    s.velocity += derivative(s).velocity * dt;
    s.displacement += s.velocity * dt;
}

Oscillator addScaled(const Oscillator& s, const Oscillator& d, float scale)
{
    return { s.displacement + d.displacement * scale, s.velocity + d.velocity * scale };
}

void stepRK4(Oscillator& s, float dt)
{
    Oscillator k1 = derivative(s);
    Oscillator k2 = derivative(addScaled(s, k1, dt * 0.5f));
    Oscillator k3 = derivative(addScaled(s, k2, dt * 0.5f));
    Oscillator k4 = derivative(addScaled(s, k3, dt));

    s.displacement += (dt / 6.0f) * (k1.displacement + 2.0f * k2.displacement + 2.0f * k3.displacement + k4.displacement);
    s.velocity += (dt / 6.0f) * (k1.velocity + 2.0f * k2.velocity + 2.0f * k3.velocity + k4.velocity);
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    const float initialDisplacement = 80.0f;
    const float omega = std::sqrt(springConstant / mass);
    Oscillator euler{ initialDisplacement, 0.0f };

    // a deliberately coarse fixed step of 10 updates per second
    const float fixedStep = 0.1f;
    float accumulator = 0.0f;
    float elapsed = 0.0f;
    sf::Clock clock;

    const sf::Vector2f anchor(450.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape ball(14.0f);
    ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    ball.setFillColor(sf::Color(200, 120, 30));

    // an outline where the mass should be according to the exact solution
    sf::CircleShape exactRing(14.0f);
    exactRing.setOrigin(sf::Vector2f(14.0f, 14.0f));
    exactRing.setFillColor(sf::Color::Transparent);
    exactRing.setOutlineColor(sf::Color(90, 90, 90));
    exactRing.setOutlineThickness(2.0f);

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        accumulator += std::min(clock.restart().asSeconds(), 0.25f);

        while (accumulator >= fixedStep)
        {
            stepSemiImplicitEuler(euler, fixedStep);
            elapsed += fixedStep;
            accumulator -= fixedStep;
        }

        float exact = initialDisplacement * std::cos(omega * elapsed);

        window.clear(sf::Color::White);

        sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + euler.displacement);

        sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
        rod[0].position = anchor;
        rod[0].color = sf::Color(120, 120, 120);
        rod[1].position = ballPosition;
        rod[1].color = sf::Color(120, 120, 120);

        ball.setPosition(ballPosition);
        exactRing.setPosition(anchor + sf::Vector2f(0.0f, restLength + exact));
        window.draw(rod);
        window.draw(exactRing);
        window.draw(ball);

        window.display();
    }
}

Task 11 - Running RK4 beside Euler

Give the program a second oscillator that uses stepRK4(), and draw both.

  1. Add a second Oscillator named rk4 with the same starting state as euler, and a blue rk4Ball.
  2. Replace the single anchor with eulerAnchor at x = 280 and rk4Anchor at x = 680.
  3. Inside the fixed-step loop, call stepRK4(rk4, fixedStep) next to the Euler call, so both integrators take exactly the same steps.
  4. Move the drawing into a drawSide lambda that takes an anchor, a state and a ball, and call it once per side. Each side draws its own ring at the same exact position.

Success condition: build and run the project. Two oscillators swing side by side from the same start. On the left the orange Euler ball drifts out of its ring as before. On the right the blue RK4 ball stays inside its ring for the same ten seconds or more.

RK4Demo.cpp
#include <SFML/Graphics.hpp>
#include <algorithm>
#include <cmath>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

Oscillator addScaled(const Oscillator& s, const Oscillator& d, float scale)
{
    return { s.displacement + d.displacement * scale, s.velocity + d.velocity * scale };
}

void stepRK4(Oscillator& s, float dt)
{
    Oscillator k1 = derivative(s);
    Oscillator k2 = derivative(addScaled(s, k1, dt * 0.5f));
    Oscillator k3 = derivative(addScaled(s, k2, dt * 0.5f));
    Oscillator k4 = derivative(addScaled(s, k3, dt));

    s.displacement += (dt / 6.0f) * (k1.displacement + 2.0f * k2.displacement + 2.0f * k3.displacement + k4.displacement);
    s.velocity += (dt / 6.0f) * (k1.velocity + 2.0f * k2.velocity + 2.0f * k3.velocity + k4.velocity);
}

void stepSemiImplicitEuler(Oscillator& s, float dt)
{
    s.velocity += derivative(s).velocity * dt;
    s.displacement += s.velocity * dt;
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    const float initialDisplacement = 80.0f;
    const float omega = std::sqrt(springConstant / mass);
    Oscillator euler{ initialDisplacement, 0.0f };
    Oscillator rk4{ initialDisplacement, 0.0f };

    // a deliberately coarse fixed step of 10 updates per second
    const float fixedStep = 0.1f;
    float accumulator = 0.0f;
    float elapsed = 0.0f;
    sf::Clock clock;

    const sf::Vector2f eulerAnchor(280.0f, 60.0f);
    const sf::Vector2f rk4Anchor(680.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape eulerBall(14.0f);
    eulerBall.setOrigin(sf::Vector2f(14.0f, 14.0f));
    eulerBall.setFillColor(sf::Color(200, 120, 30));

    sf::CircleShape rk4Ball(14.0f);
    rk4Ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    rk4Ball.setFillColor(sf::Color(50, 100, 200));

    // an outline where the mass should be according to the exact solution
    sf::CircleShape exactRing(14.0f);
    exactRing.setOrigin(sf::Vector2f(14.0f, 14.0f));
    exactRing.setFillColor(sf::Color::Transparent);
    exactRing.setOutlineColor(sf::Color(90, 90, 90));
    exactRing.setOutlineThickness(2.0f);

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        accumulator += std::min(clock.restart().asSeconds(), 0.25f);

        while (accumulator >= fixedStep)
        {
            stepSemiImplicitEuler(euler, fixedStep);
            stepRK4(rk4, fixedStep);
            elapsed += fixedStep;
            accumulator -= fixedStep;
        }

        float exact = initialDisplacement * std::cos(omega * elapsed);

        window.clear(sf::Color::White);

        auto drawSide = [&](const sf::Vector2f& anchor, const Oscillator& s, sf::CircleShape& ball)
        {
            sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + s.displacement);

            sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
            rod[0].position = anchor;
            rod[0].color = sf::Color(120, 120, 120);
            rod[1].position = ballPosition;
            rod[1].color = sf::Color(120, 120, 120);

            ball.setPosition(ballPosition);
            exactRing.setPosition(anchor + sf::Vector2f(0.0f, restLength + exact));
            window.draw(rod);
            window.draw(exactRing);
            window.draw(ball);
        };

        drawSide(eulerAnchor, euler, eulerBall);
        drawSide(rk4Anchor, rk4, rk4Ball);

        window.display();
    }
}

Task 12 - Reading the error as a number

Watching a ring is convincing, but a number is easier to compare. Add a readout to each side that shows the elapsed simulated time and the distance in pixels between the ball and the ring.

  1. Add <iostream>, <iomanip> and <sstream>, and load a font with font.openFromFile as in the Sourcing Fonts lab.
  2. Give drawSide a name and a text colour. Inside it, build a label with std::ostringstream containing the name, elapsed and std::abs(s.displacement - exact), using std::setprecision(2).
  3. Draw that label as a sf::Text above each side.

Success condition: build and run the project. Between t = 4.0 s and t = 6.0 s the Semi-Implicit Euler error swings between roughly 0 and 70 pixels as it goes in and out of phase, while the RK4 error stays between about 0.7 and 4 pixels.

RK4Demo.cpp
#include <SFML/Graphics.hpp>
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>

// one mass on a spring. displacement is measured from the rest position.
struct Oscillator
{
    float displacement;
    float velocity;
};

const float springConstant = 40.0f;
const float mass = 1.0f;

// the rate of change of each part of the state. the same struct is reused
// to hold the answer, so .displacement here means d(displacement)/dt and
// .velocity here means d(velocity)/dt, which is the acceleration.
Oscillator derivative(const Oscillator& s)
{
    return { s.velocity, -(springConstant / mass) * s.displacement };
}

Oscillator addScaled(const Oscillator& s, const Oscillator& d, float scale)
{
    return { s.displacement + d.displacement * scale, s.velocity + d.velocity * scale };
}

void stepRK4(Oscillator& s, float dt)
{
    Oscillator k1 = derivative(s);
    Oscillator k2 = derivative(addScaled(s, k1, dt * 0.5f));
    Oscillator k3 = derivative(addScaled(s, k2, dt * 0.5f));
    Oscillator k4 = derivative(addScaled(s, k3, dt));

    s.displacement += (dt / 6.0f) * (k1.displacement + 2.0f * k2.displacement + 2.0f * k3.displacement + k4.displacement);
    s.velocity += (dt / 6.0f) * (k1.velocity + 2.0f * k2.velocity + 2.0f * k3.velocity + k4.velocity);
}

void stepSemiImplicitEuler(Oscillator& s, float dt)
{
    s.velocity += derivative(s).velocity * dt;
    s.displacement += s.velocity * dt;
}

int main()
{
    sf::RenderWindow window(sf::VideoMode({ 900u, 700u }), "Semi-Implicit Euler vs RK4");

    // Arial from the Windows font folder, as in the Sourcing Fonts lab
    sf::Font font;
    if (!font.openFromFile("C:/Windows/Fonts/arial.ttf"))
    {
        std::cerr << "Failed to load font, readouts will not render." << std::endl;
    }

    const float initialDisplacement = 80.0f;
    const float omega = std::sqrt(springConstant / mass);
    Oscillator euler{ initialDisplacement, 0.0f };
    Oscillator rk4{ initialDisplacement, 0.0f };

    // a deliberately coarse fixed step of 10 updates per second
    const float fixedStep = 0.1f;
    float accumulator = 0.0f;
    float elapsed = 0.0f;
    sf::Clock clock;

    const sf::Vector2f eulerAnchor(280.0f, 60.0f);
    const sf::Vector2f rk4Anchor(680.0f, 60.0f);
    const float restLength = 140.0f;

    sf::CircleShape eulerBall(14.0f);
    eulerBall.setOrigin(sf::Vector2f(14.0f, 14.0f));
    eulerBall.setFillColor(sf::Color(200, 120, 30));

    sf::CircleShape rk4Ball(14.0f);
    rk4Ball.setOrigin(sf::Vector2f(14.0f, 14.0f));
    rk4Ball.setFillColor(sf::Color(50, 100, 200));

    // an outline where the mass should be according to the exact solution
    sf::CircleShape exactRing(14.0f);
    exactRing.setOrigin(sf::Vector2f(14.0f, 14.0f));
    exactRing.setFillColor(sf::Color::Transparent);
    exactRing.setOutlineColor(sf::Color(90, 90, 90));
    exactRing.setOutlineThickness(2.0f);

    while (window.isOpen())
    {
        while (const std::optional<sf::Event> event = window.pollEvent())
        {
            if (event->is<sf::Event::Closed>()) window.close();
        }

        accumulator += std::min(clock.restart().asSeconds(), 0.25f);

        while (accumulator >= fixedStep)
        {
            stepSemiImplicitEuler(euler, fixedStep);
            stepRK4(rk4, fixedStep);
            elapsed += fixedStep;
            accumulator -= fixedStep;
        }

        float exact = initialDisplacement * std::cos(omega * elapsed);

        window.clear(sf::Color::White);

        auto drawSide = [&](const sf::Vector2f& anchor, const Oscillator& s, sf::CircleShape& ball,
                            const char* name, sf::Color textColor)
        {
            sf::Vector2f ballPosition = anchor + sf::Vector2f(0.0f, restLength + s.displacement);

            sf::VertexArray rod(sf::PrimitiveType::LineStrip, 2);
            rod[0].position = anchor;
            rod[0].color = sf::Color(120, 120, 120);
            rod[1].position = ballPosition;
            rod[1].color = sf::Color(120, 120, 120);

            ball.setPosition(ballPosition);
            exactRing.setPosition(anchor + sf::Vector2f(0.0f, restLength + exact));
            window.draw(rod);
            window.draw(exactRing);
            window.draw(ball);

            // how far the filled ball is from where the exact solution says it should be
            std::ostringstream label;
            label << std::fixed << std::setprecision(2);
            label << name << "\nt = " << elapsed << " s\nError vs exact: " << std::abs(s.displacement - exact) << " px";
            sf::Text text(font, label.str(), 16);
            text.setFillColor(textColor);
            text.setPosition(sf::Vector2f(anchor.x - 140.0f, 20.0f));
            window.draw(text);
        };

        drawSide(eulerAnchor, euler, eulerBall, "Semi-Implicit Euler", sf::Color(170, 90, 10));
        drawSide(rk4Anchor, rk4, rk4Ball, "RK4", sf::Color(30, 60, 150));

        window.display();
    }
}

Task 13 - Halving the step and counting the cost

Change fixedStep from 0.1f to 0.05f, rebuild and run, and watch the two readouts between t = 4.0 s and t = 6.0 s again. Then answer the cost question from the code you wrote.

  1. Compare each integrator's error range with the one you saw at 0.1f.
  2. Count how many times derivative() is called per step in stepSemiImplicitEuler() and in stepRK4().

Success condition: with fixedStep at 0.05f the Euler error now peaks at roughly 25 pixels instead of 70, while the RK4 error stays below about 0.25 pixels, more than ten times smaller than before. Euler improves and RK4 improves much more. The price is that stepRK4() calls derivative() four times per step where stepSemiImplicitEuler() calls it once, so RK4 pays off most when the force is cheap to evaluate relative to the accuracy it buys.

Solution

Try the tasks above before revealing the solution.