Interpolation is widely used in mathematics and in the sciences – mathematics tends to concentrate on the theoretical side of things and the methods are then used by scientists with the confidence that what they are using works.

What is interpolation? Well an example would be if you are measuring the temperature outside on a particular day then you would take the temperature at certain times and record the temperature – but what about the times in between? We can’t record everything but we might need to know what the temperature was at one of these intermediate times – this is where interpolation comes in. Interpolation is a way of “joining the dots” within (not outside) a dataset so that estimates can be made about the behaviour. This can be done in hundreds of different ways using different methods each with their pros and cons.

In a previous post I talked about the limitations of Euler’s method – well here I’m going to talk about the limitations of polynomial interpolation and specifically the Lagrange representation and one way that the problem can be partially resolved. Don’t misunderstand me – the Lagrange representation is very useful but it, as with almost any numerical technique, has its own problems.

Let’s look at the function $f(x)=\dfrac{1}{1+25x^{2}}$ on the interval $[-1,1]$. This is what it looks like

Graph of f(x)

A graph of f(x) drawn using SAGE Math


Now given $n+1$ distinct points $x_{0}, x_{1}, \ldots , x_{n}$ in $[-1,1]$ and their corresponding $y$-values $y_{0}, y_{1}, \ldots , y_{n}$ then the Lagrange representation is defined as $$P_{n}=\sum_{j=0}^{n}y_{j}\ell_{j}$$ where $$\ell_{j}=\prod_{i=0,i\neq k}^{n}\dfrac{x-x_{i}}{x_{k}-x_{i}}$$

This representation can be proved to be unique and that $P(x_{i})=y_{i}$ but I don’t want to get sidetracked with this. So going back to the function $f(x)$ – first I’m going to choose $5$ equally-spaced points along the interval [-1,1]. The resulting interpolating polynomial looks as follows (I am not going to try to write an expression for the interpolating polynomials because they are very difficult to find explicitly and they don’t really tell us anything anyway)

Lagrange representation of f(x)

Lagrange representation of f(x) with 5 equally spaced points

The interpolating curve is shown in red and $f(x)$ is shown in blue. This is not bad, after all we only have $5$ points to work with so we might expect that as the number of points increases we get a more accurate picture – right? Well…no. As the number of points increases we have to increase the degree of the interpolating polynomial and if we choose $10$ equally spaced points on the interval $[-1,1]$ we get this picture

Lagrange representation of f(x)

Lagrange representation of f(x) with 10 equally spaced points

Things are starting to go a bit awry – maybe if we increase the number of points even further then things will start to settle down. Let’s look what happens when we have $16$ equally spaced points.

Lagrange representation of f(x)

Lagrange representation of f(x) with 16 equally spaced points

Maybe that wasn’t such a good idea. Things are starting to get so out of control that as the number of interpolation points (and consequently the degree of the interpolating polynomial) increases, the interpolating polynomial oscillates wildly between both extremely large positive and extremely large negative values near the edges of the interval $[-1,1]$. This behaviour gets so bad $n$ increases that the interpolating polynomial grows without bound near the edges of the interval – this is known as Runge’s phenomenon and makes the interpolating polynomial practically useless.

One way around this is to choose a different set of interpolation points. One of the problems is the equal spacing of the points – to resolve this in part we can choose Chebyshev nodes. Using these Chebyshev nodes we get a very different picture – the following diagram shows what things look like when $10$ points are chosen

Lagrange representation of f(x)

Lagrange representation of f(x) using 10 Chebyshev nodes

Now compare that to what we had before and we see something that is much better behaved. Has Runge’s phenomenon been eliminated? Mostly, yes – but in truth it can never be completely eliminated; the Chebyshev nodes, however, massively reduce the effects of Runge’s phenomenon. Runge’s phenomenon does not always occur but it is something that can go wrong from time-to-time, so as with all numerical methods you have to take care when applying the method to solve problems; there may be another more suitable method that needs to be applied.

There are many examples of differential equations that cannot be solved analytically – in fact, it is very rare for a differential equation to have an explicit solution. Euler’s Method is a way of numerically solving differential equations that are difficult or that can’t be solved analytically. Euler’s method, however, still has its limitations.

For a differential equation $y^{\prime}=f(x,y(x))$ with initial condition $y(x_{0})=y_{0}$ we can choose a step-length $h$ and approximate the solution to the differential equation by defining $x_{n}=x_{0}+nh$ and then for each $x_{n}$ finding a corresponding $y_{n}$ where $y_{n}=x_{n-1}+hf(x_{n-1},y_{n-1})$. This method works quite well in many cases and gives good approxiamtions to the actual solution to a differential equation, but there are some differential equations that are very sensitive to the choice of step-length $h$ as the following demonstrates.

Let’s look at the differential equation $y^{\prime}+110y=100$ with initial condition $y(0)=2$.

This differential equation has an exact solution given by $y=1+\mathrm{e}^{-100t}$ but this example is a very good example which demonstrates that Euler’s method cannot be used blindly. Let’s look at what happens for a few different step-lengths.

For the step-length $h=0.019$ step-length we get the following behaviour

Behaviour of numerical solution of an ODE

The red curve is the actual solution and the blue curve represents the behaviour of the numerical solution given by the Euler method – it is clear that the numerical solution converges to the actual solution so we should be very happy. However, look what happens when the step-length $h=0.021$ is chosen

Behaviour of numerical solution of an ODE

Again the actual solution is represented by the red line which on this diagram looks like a flat line because the blue curve gets bigger and bigger as you move along the $x$-axis. So a change of just $0.002$ in the step-length has completely changed the behaviour of the numerical solution. For a step-length $h=0.03$ the graph would look as follows

Behaviour of numerical solution of an ODEThe actual solution can barely be seen and the numerical solution gets out of control very quickly – this solution is completely useless – the scales on the $y$-axis are enormous and increasing the step-length only makes this worse. What has happened?

It can be shown by induction that for $n \in \mathbb{N}$ that $y_{n}=1+(1-100h)^{n}$. This converges only for $h<0.02$ and diverges for $h>0.02$. $h=0.02$ is a limiting case and gives an oscillating numerical solution that looks as follows

Behaviour of numerical solution of an ODE

For this particular example for $h<0.02$ and as the step-length gets closer to $0$ the solution will converge faster and for $h>0.02$ as the step-length increases the solution will diverge more rapidly.

So even though we have Euler’s method at our disposal for differential equations this example shows that care must be taken when dealing with numerical solutions because they may not always behave as you want them to. This differential equation is an example of a stiff equation – in other words, one that is very sensitive to the choice of step length. In general as the step-length increases the accuracy of the solution decreases but not all differential equations will be as sensitive to the step-length as this differential equation – but they do exist.

I mentioned vector fields in a previous post in the context of differential equations and over the last week or so I have been looking at them in a bit more detail. Vector fields sound quite complicated but they can be very simple. A vector field can be presented visually as a vector attached to each point in space. The space may be the $x$-$y$ plane, three-dimensional space, it could be a region of the $x$-$y$ plane or even a manifold. A typical vector field in 2-dimensions might look as follows.

A simple vector field

An image of a simple vector field drawn in SAGE Math

The arrows represent the vector that is attached to that particular point – the direction of the arrow gives the direction of the vector and the size of the arrow gives an idea of its relative magnitude. Graphical representations of vector fields can be a little misleading as it is tempting to think that only certain points have vectors attached to them – this is not the case; every point has a vector attached to it but if we were to try to show all of them the diagram would be too cluttered.

Vector fields are useful to model flows of liquids or gases; for example in weather prediction a vector field that changes over time could be used to model wind patterns. The vector attached to each point would tell you the direction and strength of the wind at that point and the vector field would evolve from one moment to the next. If you define a surface in a vector field then you can use integration to measure the flux across the surface – the physical interpretation of flux is as a measure of the amount of substance flowing across a surface. I came across this question in the book Advanced Engineering Mathematics Fifth Edition by Stroud and Booth and decided to give it a go.

Evaluate $\int_{S}\mathbf{F}.\mathrm{d}\mathbf{S}$ over the surface $S$ defined by $x^{2}+y^{2}+z^{2}=4$ for $z\geq0$ and bounded by $x=0$, $y=0$, $z=0$ and $\mathbf{F}=x\mathbf{i}+2z\mathbf{j}+y\mathbf{k}$.

The pictures below give an idea of what the vector field and surface both look like from a few different angles. This problem is asking to integrate this vector field over the surface to find the flux across the surface.

Flow of a vector field across a surface Flow of a vector field across a surface Flow of a vector field across a surface


For this problem, since the surface that we are integrating over is part of a sphere, it is convenient to change to spherical polar co-ordinates given by$$x=r\mathrm{sin}\theta\mathrm{cos}\phi$$ $$y=r\mathrm{sin}\theta\mathrm{sin}\phi$$ $$z=r\mathrm{cos}\theta$$.The integration itself is quite straightforward although some of the integrands look a bit of a pain at first glance, but some techniques from A-Level Further Maths courses should clear things up. I used some reduction formulae to deal with some of the integrals that I ended up with which really simplified things (which is always good). You can download and view my full solution here – Integrating Vector Fields.

Some differential equations are relatively easy to solve – even though at A-Level this may seem to be the norm, in reality it is very rare for a differential equation to have an explicit solution. The types of differential equations at A-Level can be solved by direct integration or by using a range of nifty techniques like making a substitution to change the differential equation into one that can be solved directly, separation of the variables can be used to solve some equations of the form $f(y)y’=g(x)$; sometimes the integrating factor $\mathrm{e}^{\int P(x) \mathrm{d}x}$ can be used to solve equations of the form $y^{\prime}+P(x)y=Q(x)$ and second-order equations of the form $ay^{\prime \prime}+by’+cy=f(x)$ where $a,b$ and $c$ are constants, can be solved by finding the complementary function and particular integral and combining the two.

However, these are very special cases. Many differential equations that come up in real-life situations are very non-linear in nature. Non-linear differential equations are generally very difficult to solve, but quite often, it is the behaviour of the solutions that interests us more than the actual solution. We may not be really interested in a particular solution but a family of solutions so that we can easily compare their behaviours particularly for equations that evolve over time. In some cases – even though it may be very difficult or even impossible to find the solution to a differential equation analytically, we can still analyse the behaviour through use of a phase-diagram.

A phase-diagram is a vector field that we can use to visually present the solutions to a differential equation. For example here is a second-order differential equation – (this is an example that I did that appears in the book by D. W. Jordan and P. Smith titled Nonlinear Ordinary Differential Equations – An Introduction for Scientists and Engineers Fourth Edition)
$$ \ddot{x} = x-x^{2}$$

This second order-differential equation can be separated into a system of first-order differential equations given by $$\dot{x}=y$$ $$\dot{y}=x-x^{2}$$

This is a very common technique for solving differential equations since first-order ODEs are generally easier to solve than second-order ODEs. Note that $$\frac{\mathrm{d}y}{\mathrm{d}x}= \frac{\dot{y}}{\dot{x}}=\frac{x-x^{2}}{y}$$

The phase-portrait of this looks as follows

Phase Portrait for a second order differential equation

A phase portrait drawn in SAGE Math for a second-order ODE

At a point on the phase-diagram there is attached a vector which are represented here by the arrows – the vector tells you which direction to move in. So given a starting point it is simply a case of moving from point to point in the direction that the vector tells you to move in. So different starting points will give different paths through the phase-portrait called trajectories, some of which are shown in the following diagram

Phase Portrait for a second order differential equation

Phase portrait of a second-order ODE with trajectories

Each trajectory is described by the equation $3y^{2}=-2x^{3}+3x^{2}+C$ for different values of $C$. The trajectories above are for $C=-8$, $-4$, $-0.5$, $0$, $1$, $2$, $4$ and $6$. Here the horizontal $x$-axis represents the position of a particle and the vertical $y$-axis represents the velocity of the particle – remember that $y=\dot{x}$.

So now if we start in the top left hand corner of the diagram on the curve corresponding to when $C=8$ (the right-most blue line) the arrows tell us to move to the right and down – so our displacement increases as we move to the right but our velocity decreases as we move down. Eventually we cross the $x$-axis and here our velocity is instantaneously zero; then we start moving backwards since our velocity becomes negative as we continue moving downwards, and our displacement decreases as we move to the left. So I have not even calculated the actual solution to the ODE but I still know how the solution behaves.


A Tautochrone is a curve described parametrically by the equations $x=a(\theta+\mathrm{sin}\theta)$ and $y=a(1-\mathrm{cos}\theta)$. Here is a graph of the curve for $\theta \in (-2\pi, 2\pi)$

A Tautochrone drawn using SAGE Math

The graph of a tautochrone drawn in SAGE math

This curve is not commonly encountered in mathematics, indeed, I only came across this curve for the first time about two months ago – but it has some very interesting properties. The curve itself is quite simple; however, if a particle were to be released from rest on one of the slopes of the tautochrone, then the time taken for the particle to reach the bottom of the slope is independent of its starting position assuming that the only force acting on it is the gravitational force. In other words, if you were to simultaneously set a ball rolling down the slope from the top of the slope and another ball from halfway down the slope then they would both arrive at the bottom of the slope at exactly the same time.

This property can be proved using energy considerations and some basic trigonometric identities to form the differential equation


where $g$ is the gravitational force – and then using integration by substitution we find that the time taken for the particle to reach the bottom is $T=\pi\sqrt{\frac{a}{g}}$ which is independent of starting position – you can download my full, detailed proof of this property here – Basic Properties of the Tautochrone. The derivation of the differential equation and then the integration that follows can all be done using techniques from the A-Level maths and further maths courses.

Since we are assuming that the only force acting on the particle is the gravitational force we can assume that all gravitational potential energy lost (remember that the particle will move down the slope and therefore lose gravitational potential energy) will be converted to kinetic energy. This is the starting point of the whole derivation of the above differential equation and although the resulting differential equation is non-linear, we are fortunate that it is nice enough to be able to solve – non-linear differential equations are notoriously difficult to solve and often impossible to solve analytically.