This script solves the equation of motion for the acceleration phase of the Bloodhound supersonic car using Chebfun.

The car is powered by a jet engine and a rocket and aims to reach a speed in excess of 1000 mph. More details on the project can be found at the website

Mathematically, the problem takes the form of an initial value ODE

$$ a(t)v' + b(t)v + cv^2 = d(t), \qquad t \in [0,50],~~ v(0) = 0 . $$

Here $a$ is a piecewise linear function and $b$ and $d$ are piecewise constants. The break occurs at the discontinuity caused by the rocket ignition.

One can vary the parameters of the run, such as the rocket start time or the mass of the car, to see the effect on the time taken to reach 1000 miles per hour.

Define the domain.

dom = [0 50];
RocketStart = 11;      % in seconds

Compute the mass.

MassOriginal = 6500;   % in kg
JetThrust = 80;        % in kN
RocketThrust = 115;    % in kN
JSFC = 0.0005*102;     % in kg/(kN*h) Jet Specific Fuel Consumption
RSFC = 102/220;        % Rocket Specific Fuel Consumption

The mass satisfies

$$ \mathrm{mass} = \mathrm{original} - \mathrm{(jet~fuel~burned)} - \mathrm{(rocket~fuel~burned)} $$

mass = @(t) MassOriginal - (JSFC*JetThrust*t) - ...   % in kg, t in seconds

Create a chebfun for mass to capture the piecewise linear behaviour.

cmass = chebfun(mass, dom, 'splitting', 'on');

Create chebfuns or anonymous functions for the forces.

thrust = @(t) 1000*(JetThrust + RocketThrust.*(t>RocketStart)); % in N

Create a chebfun for thrust to capture the piecewise constant behavior.

cthrust = chebfun(thrust, dom, 'splitting', 'on');


aerodrag = @(v) (175./289).*v.^2; % in N
surfacedrag = (2/5)*cmass*9.81;   % in N

Create a chebop to represent the differential equation. Newton's 2nd law (variable mass):

$$ \frac{d(mv)}{dt} = \mathrm{thrust} - \mathrm{drag}, $$ $$ m \frac{dv}{dt} + \frac{dm}{dt} v = \mathrm{thrust} - \mathrm{drag}. $$

N = chebop(@(t,v) cmass.*diff(v) + diff(cmass).*v - cthrust + ...
    aerodrag(v) + surfacedrag, dom);

Set up the right-hand side of the differential equation so that $ N(v) = \mathrm{rhs} $.

rhs = 0;

Assign boundary conditions to the chebop.

N.bc = @(t,v) v(0);

Define an initial solution.

N.init = chebfun(@(t) t, dom);

Solve the problem using backslash.

u = N\rhs;
u_mph = u/0.44704; % convert to mph for plotting
Warning: Newton iteration failed.
Please try supplying a better initial guess via the .init field
of the chebop. 

Find where the solution is equal to 1000 mph.

t1000 = find(u == 447); % number of seconds when we hit 1000 mph

Plot the solution.

hold on
plot([0 50],[1000 1000],'r', 'LineWidth', 2); % 1000 mph line
hold on;
text(20, 500, ['Time to 1000 mph = ' num2str(t1000, 3) ' seconds']);
title('Velocity of the Bloodhound Supersonic Car (Acceleration Phase)')
ylabel('velocity in mph')

Integrate to find the distance.

s = cumsum(u);
s_miles = s/1609;     % convert to miles to plot
plot(s_miles,'LineWidth',2); title('Distance travelled')
xlabel('Time in seconds'),   ylabel('Distance in miles')