An Essay on the Numerical Solution of Differential Equations - for Students

Don Estep and Claes Johnson


This essay is excerpted from the Introduction to our book, Computational Differential Equations, K. Eriksson, D. Estep, P. Hansbo, C. Johnson, Cambridge University Press, 1996, ISBN 91-44-49311-8. Follow the link to download a postscript file containing the table of contents, a chapter, and the index.

Knowing thus the Algorithm of this calculus, which I call Differential Calculus, all differential equations can be solved by a common method. (Leibniz)

When, several years ago, I saw for the first time an instrument which, when carried, automatically records the number of steps taken by a pedestrian, it occurred to me at once that the entire arithmetic could be subjected to a similar kind of machinery so that not only addition and subtraction, but also multiplication and division could be accomplished by a suitably arranged machine easily, promptly and with sure results.... For it is unworthy of excellent men to lose hours like slaves in the labour of calculations, which could safely be left to anyone else if the machine was used.... And now that we may give final praise to the machine, we may say that it will be desirable to all who are engaged in computations which, as is well known, are the mangers of financial affairs, the administrators of others estates, merchants, surveyors, navigators, astronomers, and those connected with any of the crafts that use mathematics. (Leibniz)

Building on tradition going back to the ancient Greek philosophers, Leibniz and Newton invented calculus in the late 17th century and thereby laid the foundation for the revolutionary development of science and technology that is shaping the world today. Calculus is a method for modeling physical systems mathematically by relating the state of a system to its neighboring states in space-time using differential and integral equations. Because calculus is inherently computational, this revolution began to accelerate tremendously in the 1940s when the electronic computer was created. Today, we are seeing what is essentially a ``marriage'' of calculus and computation in the creation of the field of computational mathematical modeling.

Actually, Leibniz himself sought to realize a unification of calculus and computation, but failed because the mechanical calculator he invented was not sufficiently powerful. The next serious effort was made in the 1830s by Babbage, who designed a steam powered mechanical computer he called the Analytical Engine. Again, technical difficulties and low speed stopped his ambitious plans.

The possibility of realizing Leibniz' and Babbage's visions of a universal computing machine came with the invention of the electronic valve in the 1930s, which enabled the construction of high speed digital computers. The development took a leap during the World War II spurred by the computing demands of the military. Until this time, large scale computations were performed by rooms of people using mechanical adding machines. The war provided an immediate pressure to speed up the process of scientific development by using mathematical modeling to hone a physical problem down to a manageable level, and mathematicians and physicists became interested in the invention of an electronic computing device. The logical design of programmable electronic computers was developed by the mathematician von Neumann, among others. By the late forties, von Neumann was using the first ENIAC (Electronic Numerical Integrator And Calculator) computer to address questions in fluid dynamics and aerodynamics.

The subsequent development of computer power that has resulted in desktop computers of far greater power than the ENIAC, has been paralleled by the rapid introduction of computational mathematical modeling into all areas of science and engineering. Questions routinely addressed computationally using a computer include: What is the weather going to do in three days? Will this airplane fly? Can this bridge carry a load of ten trucks? What happens during a car collision? How do we direct a rocket to pass close by Saturn? How can we create an image of the interior of the human body using very weak X-rays? What is the shape of a tennis racket that has the largest "sweet spot"? What is a design of a bicycle frame that combines low weight with rigidity? How can we create a sharp picture from a blurred picture? What will the deficit be in Sweden in the year 2000? How much would the mean temperature of the earth increase if the amount of carbon dioxide in the atmosphere increased by 20 percent?

The physical situations behind these kinds of questions are modeled by expressing the laws of mechanics and physics (or economics) in terms of equations that relate derivatives and integrals. Common variables in these models are time, position, velocity, acceleration, mass, density, momentum, energy, stress and force, and the basic laws express conservation of mass, momentum and energy, and balance of forces. Information about the physical process being modeled is gained by solving for some of the variables in the equation, i.e. by computing the solution of the differential/integral equation in terms of the others, which are assumed to be known data. Calculus is the basic study of differential/integral equations and their solutions.

Sometimes it is possible to find an exact formula for the solution of a differential/integral equation. For example, the solution might be expressed in terms of the data as a combination of elementary functions or as a trigonometric or power series. This is the classical mathematical method of solving a differential equation, which is now partially automated in mathematical software for symbolic computation such as Maple or Mathematica. However, this approach only works on a relatively small class of differential equations. In more realistic models, solutions of differential equations cannot be found explicitly in terms of known functions, and the alternative is to determine an approximate solution for given data through numerical computations on a computer. The basic idea is to discretize a given differential/integral equation to obtain a system of equations with a finite number of unknowns, which may be solved using a computer to produce an approximate solution. The finite-dimensional problem is referred to as a discrete problem and the corresponding differential/integral equation as a continuous problem. A good numerical method has the property that the error decreases as the number of unknowns, and thus the computational work, increases. Discrete problems derived from physical models are usually computationally intensive, and hence the rapid increase of computer power has opened entirely new possibilities for this approach. Using a desktop computer, we can often obtain more information about physical situations by numerically solving differential equations than was obtained over all the previous centuries of study using analytical methods.

The progress in weather prediction is a good example for this discussion. Historically, weather forecasting was based on studying previous patterns to predict future behavior. A farmer's almanac gives predictions based on the past behavior, but involves so many variables related to the weather that determining meaningful correlations is an overwhelming task. By modeling the atmosphere with a set of differential equations, the number of variables is reduced to a handful that can be measured closely, albeit at many locations. This was envisioned by the English pioneer of numerical weather prediction Richardson in the 1920s, who proposed the formation of a department of 64,000 employees working in shifts to perform the necessary calculations using mechanical calculators more quickly than the weather changed. After this proposal, the attitude toward numerical weather prediction became pessimistic. Not until the development of the modern computer, could the massive computations required be performed sufficiently rapidly to be useful. The first meaningful numerical forecasts were made by von Neumann and Charney in the late forties using the ENIAC, but of course the reliability was very low due to the extremely coarse discretization of the earth's system they had to use. The most recent model for the global weather uses a discretization grid with roughly 50,000 points horizontally and 31 layers vertically giving a total of five million equations that are solved in a couple of hours on a super-computer.

There are three sources of errors affecting the reliability of a numerical weather forecast: (i) measurement errors in data (or lack of data) (ii) approximation errors in modeling and (iii) approximation errors in computation. The initial data at the start of the computer simulation are always measured with some error; the set of differential equations in the computer model only approximately describes the evolution of the atmosphere; and finally the numerical solution of the differential equations is only an approximation of the true solution. These sources add up to form the total prediction error. It is essential to be able to estimate the total error by estimating individually the contributions from the sources (i)-(iii) and improve the precision where possible. This is a basic issue in all applications in computational mathematical modeling.

Our experience tells that forecasts of the daily weather become very unreliable in predictions for more than say a week. This was discussed in the 1960s by the meteorolgist Lorenz, who coined the phrase "the butterfly effect" to describe situations in which a small cause can have a large effect after some time. Lorenz gave a simple example displaying this phenomenon in the form of the Lorenz system of ordinary differential equations with only three unknowns. We plot a typical solution below, showing the trajectory of a "particle" being ejected away from the origin to be attracted into a slowly diverging orbit to the left, then making a loop on the right, returning to a few orbits to the left, then back to the right etc. The trajectory is very sensitive to perturbations as to the number of loops to the left or right, and thus is difficult to compute accurately over a longer time interval, just as the evolution of the weather may be difficult to predict for more than a week.

Evolution of the Solution of
the Lorenz Equations

A solution of the Lorenz system computed with an error of .1 or less over the time interval (0,30).

If we summarize the Leibniz vision as a fusion of mathematical modeling, mathematical analysis and computation, then there are three fundamental issues to be addressed:

Our books try to answer these questions for a specific set of problems and provide a set of tools that can be used to tackle the large variety of problems met in applications.


Today we know that the program-controlled computer began during the last century with Babbage. But he was so far ahead his time that his machine was nearly completely forgotten. So in Germany when I started in 1934 nobody knew of his work. I was a student then in civil engineering in Berlin. Berlin is a nice town and there were many opportunities for a student to spend his time in an agreeable manner, for instance with nice girls. But instead of that we had to perform big and awful calculations. Also later as an engineer in the aircraft industry, I became aware of the tremendous number of monotonous calculations necessary for the design of static and aerodynamic structures. Therefore I decided to design and construct calculating machines suited to solving these problems automatically. (Konrad Zuse)