# One-step process

One-step methods approximate the solution (blue) of an initial value problem by successively determining points , etc. from the given starting point${\ displaystyle A_ {0}}$${\ displaystyle A_ {1}}$${\ displaystyle A_ {2}}$

In numerical mathematics, one-step methods are, in addition to the multi-step methods, a large group of calculation methods for solving initial value problems . This task, in which a common differential equation is given together with a starting condition, plays a central role in all natural and engineering sciences and is also gaining more and more importance, for example, in economics and social sciences. Initial value problems are used to analyze, simulate, or predict dynamic processes.

The name-giving basic idea of ​​the one-step method is that, starting from the given starting point, it calculates approximation points step by step along the sought solution. In doing so, they only use the most recently determined approximation for the next step, in contrast to the multi-step method, which also include points that were further back in the calculation. The one-step methods can be roughly divided into two groups: the explicit methods, which calculate the new approximation directly from the old one, and the implicit methods, in which an equation has to be solved. The latter are also suitable for so-called stiff initial value problems.

The simplest and oldest one-step method, the explicit Euler method , was published by Leonhard Euler in 1768 . After a group of multi-step processes had been presented in 1883, Carl Runge , Karl Heun and Wilhelm Kutta developed significant improvements to Euler's process around 1900 . From these emerged the large group of the Runge-Kutta processes , which form the most important class of one-step processes. Further developments of the 20th century are, for example, the idea of extrapolation , but above all considerations for step size control, i.e. for the selection of suitable lengths for the individual steps of a method. These concepts form the basis to be able to solve difficult initial value problems, as they occur in modern applications, efficiently and with the required accuracy using computer programs.

## introduction

### Ordinary differential equations

The development of differential and integral calculus by the English physicist and mathematician Isaac Newton and independently of it by the German polymath Gottfried Wilhelm Leibniz in the last third of the 17th century was an essential impetus for the mathematization of science in the early modern period . These methods formed the starting point of the mathematical branch of analysis and are of central importance in all natural and engineering sciences. While Leibniz was led by the geometric problem of determining tangents to given curves to differential calculus, Newton started from the question of how changes in a physical quantity can be determined at a certain point in time.

For example, when a body moves, its average speed is simply the distance traveled divided by the time it took. However, in order to mathematically formulate the current speed of the body at a certain point in time , a limit crossing is necessary: ​​one considers short periods of length , the distances covered and the associated average speeds . If you let the time span converge to zero and if the average speeds also approach a fixed value , then this value is called the (instantaneous) speed at the given point in time . Designates the position of the body at the time , then one writes and names the derivative of . ${\ displaystyle v (t)}$${\ displaystyle t}$${\ displaystyle \ Delta t}$${\ displaystyle \ Delta x}$${\ displaystyle {\ tfrac {\ Delta x} {\ Delta t}}}$${\ displaystyle \ Delta t}$${\ displaystyle v (t)}$${\ displaystyle t}$${\ displaystyle x (t)}$${\ displaystyle t}$${\ displaystyle v (t) = x '(t)}$${\ displaystyle v}$${\ displaystyle x}$

The decisive step in the direction of differential equation models is now the opposite question: In the example of the moving body, the speed is known at all times and its position is to be determined from this. It is clearly clear that the initial position of the body must also be known at a point in time in order to be able to solve this problem clearly. It is a function of searching that the initial condition with given values and fulfilled. ${\ displaystyle t}$${\ displaystyle v (t)}$${\ displaystyle x (t)}$${\ displaystyle t_ {0}}$${\ displaystyle x (t)}$${\ displaystyle x '(t) = v (t)}$${\ displaystyle x (t_ {0}) = x_ {0}}$${\ displaystyle t_ {0}}$${\ displaystyle x_ {0}}$

In the example of determining the position of a body from its speed, the derivation of the function sought is explicitly given. Usually, however, there is the important general case of ordinary differential equations for a desired quantity : Based on the laws of nature or the model assumptions, a functional relationship is known that indicates how the derivation of the function to be determined can be calculated from and from the (unknown) value . In addition, an initial condition must be given again , which can be obtained, for example, from a measurement of the required variable at a fixed time. In summary, we have the following general type of problem: Find the function that contains the equations ${\ displaystyle x}$${\ displaystyle y}$${\ displaystyle y '(t)}$${\ displaystyle t}$${\ displaystyle y (t)}$${\ displaystyle y}$

${\ displaystyle y '(t) = f (t, y (t)), \ quad y (t_ {0}) = y_ {0}}$

where is a given function. ${\ displaystyle f}$

The solution of the differential equation of the Lorenz attractor shown is a very complicated curve in three-dimensional space

A simple example is a size that grows exponentially . This means that the current change, i.e. the derivative, is proportional to itself. So it applies with a growth rate and, for example, an initial condition . In this case, the solution you are looking for can already be found using elementary differential calculus and specified using the exponential function : It applies . ${\ displaystyle y}$${\ displaystyle y '(t)}$ ${\ displaystyle y (t)}$${\ displaystyle y '(t) = \ lambda y (t)}$${\ displaystyle \ lambda}$${\ displaystyle y (0) = y_ {0}}$${\ displaystyle y}$${\ displaystyle y (t) = y_ {0} e ^ {\ lambda t}}$

The function sought in a differential equation can be vector-valued , that is, for each there can be a vector with components. One then speaks of a -dimensional differential equation system. In the case of a moving body, its position in -dimensional Euclidean space and its speed are at the time . The differential equation therefore specifies the speed of the trajectory sought with direction and magnitude at each point in time and space . The path itself should be calculated from this. ${\ displaystyle y}$${\ displaystyle t}$${\ displaystyle y (t) = (y_ {1} (t), \ dotsc, y_ {d} (t))}$${\ displaystyle d}$${\ displaystyle d}$${\ displaystyle y (t)}$${\ displaystyle d}$${\ displaystyle y '(t)}$${\ displaystyle t}$

### Basic idea of ​​the one-step process

In the simple differential equation of exponential growth considered as an example above, the solution function could be specified directly. This is generally no longer possible with more complex problems. Under certain additional conditions, one can then show for the function that a uniquely determined solution of the initial value problem exists; However, this can then no longer be calculated explicitly using analysis methods (such as separation of variables , an exponential approach or variation of the constants ). In this case, numerical methods can be used to determine approximations for the solution sought. ${\ displaystyle f}$

The methods for the numerical solution of initial value problems of ordinary differential equations can be roughly divided into two large groups: the one-step and the multi-step method. Both groups have in common that they gradually calculate approximations for the desired function values at points . The defining property of the one-step method is that only the “current” approximation is used to determine the following approximation . In contrast, in the case of multi-step procedures, previously calculated approximations are also included; a three-step method would therefore, for example, also use and to determine the new approximation . ${\ displaystyle y_ {0}, y_ {1}, y_ {2}, \ dotsc}$${\ displaystyle y (t_ {0}), y (t_ {1}), y (t_ {2}), \ dotsc}$${\ displaystyle t_ {0} ${\ displaystyle y_ {j + 1}}$${\ displaystyle y_ {j}}$${\ displaystyle y_ {j}}$${\ displaystyle y_ {j-1}}$${\ displaystyle y_ {j-2}}$${\ displaystyle y_ {j + 1}}$

Two steps of the explicit Euler method

The simplest and most fundamental one-step method is the explicit Euler method , which the Swiss mathematician and physicist Leonhard Euler presented in his textbook Institutiones Calculi Integralis in 1768 . The idea of ​​this method is to approximate the solution sought by a piece-wise linear function in which the slope of the straight line is given in every step from point to point . More precisely: The problem already gives a value for the function sought, namely . But the derivation at this point is also known, because it applies . The tangent to the graph of the solution function can thus be determined and used as an approximation. At this point there is a step size${\ displaystyle t_ {j}}$${\ displaystyle t_ {j + 1}}$${\ displaystyle f (t_ {j}, y_ {j})}$${\ displaystyle y (t_ {0}) = y_ {0}}$${\ displaystyle y '(t_ {0}) = f (t_ {0}, y_ {0})}$${\ displaystyle t_ {1}> t_ {0}}$${\ displaystyle h_ {0}: = t_ {1} -t_ {0}}$

${\ displaystyle y (t_ {1}) \ approx y_ {0} + h_ {0} f (t_ {0}, y_ {0}) =: y_ {1}}$.

This procedure can now be continued in the following steps. Overall, this results in the calculation rule for the explicit Euler method

${\ displaystyle y_ {j + 1} = y_ {j} + h_ {j} f (t_ {j}, y_ {j}), \ quad j = 0,1,2, \ dotsc}$

with the step sizes . ${\ displaystyle h_ {j} = t_ {j + 1} -t_ {j}}$

The explicit Euler method is the starting point for numerous generalizations in which the slope is replaced by slopes that approximate the behavior of the solution between the points and more precisely. The implicit Euler method, which is used as the slope, provides an additional idea for one-step methods . At first glance, this choice seems unsuitable, as it is unknown. However, the equation is now obtained as a process step ${\ displaystyle f (t_ {j}, y_ {j})}$${\ displaystyle t_ {j}}$${\ displaystyle t_ {j + 1}}$${\ displaystyle f (t_ {j + 1}, y_ {j + 1})}$${\ displaystyle y_ {j + 1}}$

${\ displaystyle y_ {j + 1} = y_ {j} + h_ {j} f (t_ {j + 1}, y_ {j + 1})}$

from which (if necessary by a numerical method) can be calculated. If, for example, the arithmetic mean of the slopes of the explicit and the implicit Euler method is selected as the slope, then the implicit trapezoid method is obtained . An explicit method can be obtained from this, for example, if the unknown on the right-hand side of the equation is approximated by the explicit Euler method, the so-called Heun method . All these procedures and all further generalizations have the basic idea of ​​the one-step procedure in common: the step ${\ displaystyle y_ {j + 1}}$${\ displaystyle y_ {j + 1}}$

${\ displaystyle y_ {j + 1} = y_ {j} + h_ {j} \ Phi}$

with a slope that can depend on , and and (in the case of implicit methods) on. ${\ displaystyle \ Phi}$${\ displaystyle t_ {j}}$${\ displaystyle y_ {j}}$${\ displaystyle h_ {j}}$${\ displaystyle y_ {j + 1}}$

## definition

With the considerations from the introductory section of this article, the concept of the one-step procedure can be defined as follows: We are looking for the solution to the initial value problem ${\ displaystyle y}$

${\ displaystyle y '(t) = f (t, y (t))}$, .${\ displaystyle \ quad y (t_ {0}) = y_ {0}}$

It is assumed that the solution

${\ displaystyle y \ colon I \ to \ mathbb {R} ^ {d}}$

exists on a given interval and is uniquely determined. are ${\ displaystyle I = [t_ {0}, T]}$

${\ displaystyle t_ {0}

Intermediate points in the interval and the associated step sizes, then that means through ${\ displaystyle I}$${\ displaystyle h_ {j} = t_ {j + 1} -t_ {j}}$

${\ displaystyle y_ {j + 1} = y_ {j} + h_ {j} \ Phi (t_ {j}, y_ {j}, y_ {j + 1}, h_ {j})}$, ${\ displaystyle \ quad j = 0, \ dotsc, n-1}$

given procedures one-step procedure with procedure function . If it doesn't depend on, then it's called an explicit one-step procedure . Otherwise an equation for has to be solved in each step and the procedure is called implicit . ${\ displaystyle \ Phi}$${\ displaystyle \ Phi}$${\ displaystyle y_ {j + 1}}$${\ displaystyle j}$${\ displaystyle y_ {j + 1}}$

## Consistency and Convergence

### Convergence order

Global error at different step sizes for three one-step methods. In a
double-logarithmic plot , the relationship appears approximately linear, the slopes correspond to the convergence orders 1, 2 and 4.${\ displaystyle h}$

For a practicable one-step method, the calculated good approximations for the values ​​of the exact solution should be at the point . Since the quantities are generally -dimensional vectors, the quality of this approximation is measured with a vector norm as , the error at the point . It is desirable that these errors quickly converge to zero for all if the step sizes are allowed to converge to zero. In order to also cover the case of non-constant step sizes, one defines more precisely as the maximum of the used step sizes and considers the behavior of the maximum error at all points in comparison to powers of . It is said that the one-step procedure for solving the given initial value problem has the order of convergence if the estimate ${\ displaystyle y_ {j}}$${\ displaystyle y (t_ {j})}$${\ displaystyle y}$${\ displaystyle t_ {j}}$${\ displaystyle d}$${\ displaystyle \ | y_ {j} -y (t_ {j}) \ |}$${\ displaystyle t_ {j}}$${\ displaystyle j}$${\ displaystyle h}$${\ displaystyle t_ {j}}$${\ displaystyle h}$ ${\ displaystyle p \ geq 1}$

${\ displaystyle \ max _ {j = 0, \ dotsc, n} \ | y_ {j} -y (t_ {j}) \ | \ leq Ch ^ {p}}$

holds for all sufficiently small with a constant that is independent of . ${\ displaystyle h}$${\ displaystyle h}$${\ displaystyle C> 0}$

The order of convergence is the most important parameter for comparing different one-step methods. A method with a higher order of convergence generally delivers a smaller total error for a given step size or, conversely, fewer steps are required to achieve a given accuracy. In the case of a method with, it is to be expected that if the step size is halved, the error is only roughly halved. In the case of a method of the order of convergence , however, one can assume that the error is reduced by approximately the factor . ${\ displaystyle p}$${\ displaystyle p = 1}$${\ displaystyle p = 4}$${\ displaystyle {\ bigl (} {\ tfrac {1} {2}} {\ bigr)} ^ {4} = {\ tfrac {1} {16}}}$

### Global and local error

The errors considered in the definition of the order of convergence are made up of two individual components in a manner that initially appears complicated: Of course, they depend on the one hand on the error that the procedure makes in a single step, in that the unknown slope of the function sought is replaced by the Procedural function approximates. On the other hand, it must also be taken into account that the starting point of a step generally does not already coincide with the exact starting point ; the error after this step also depends on all errors that have already been made in the previous steps. Due to the uniform definition of the one-step procedure, which only differs in the choice of the procedural function , it can be proven that (under certain technical conditions ) one can infer the convergence order directly from the error order in a single step, the so-called consistency order. ${\ displaystyle \ | y_ {j} -y (t_ {j}) \ |}$${\ displaystyle (t_ {j}, y_ {j})}$${\ displaystyle (t_ {j}, y (t_ {j}))}$${\ displaystyle \ Phi}$${\ displaystyle \ Phi}$

The concept of consistency is a general and central concept in modern numerical mathematics. While one examines how well the numerical approximations fit the exact solution during the convergence of a procedure, when it comes to consistency, put simply, one asks the "opposite" question: How well does the exact solution fulfill the procedural specification? In this general theory it holds that a method is convergent if and only if it is consistent and stable . In order to simplify the notation, it should be assumed in the following consideration that an explicit one-step procedure

${\ displaystyle y_ {j + 1} = y_ {j} + h \ Phi (t_ {j}, y_ {j}, h)}$

with a constant step size . The true solution defines the local clipping error (also called local procedural error) as ${\ displaystyle h}$${\ displaystyle t \ mapsto y (t)}$${\ displaystyle \ eta}$

${\ displaystyle \ eta (t, h) = y (t) + h \ Phi (t, y (t), h) -y (t + h)}$.

So one assumes that the exact solution is known, starts a process step at that point and forms the difference to the exact solution at the point . This defines: A one-step procedure has the order of consistency if the estimate ${\ displaystyle (t, y (t))}$${\ displaystyle t + h}$ ${\ displaystyle p \ geq 1}$

${\ displaystyle \ | \ eta (t, h) \ | \ leq Ch ^ {p + 1}}$

holds for all sufficiently small with a constant that is independent of . ${\ displaystyle h}$${\ displaystyle h}$${\ displaystyle C> 0}$

The noticeable difference between the definitions of the order of consistency and order of convergence is the power instead of . This can be clearly interpreted in such a way that a power of the step size is “lost” during the transition from the local to the global error. The following principle applies, which is central to the theory of the one-step procedure: ${\ displaystyle h ^ {p + 1}}$${\ displaystyle h ^ {p}}$

If the process function is Lipschitz continuous and the associated one-step process has the order of consistency , then it also has the order of convergence .${\ displaystyle \ Phi}$ ${\ displaystyle p}$${\ displaystyle p}$

The Lipschitz continuity of the process function as an additional requirement for stability is generally always fulfilled when the function from the differential equation is itself Lipschitz continuous. This requirement has to be assumed for most applications anyway, in order to guarantee the unambiguous solvability of the initial value problem. According to the sentence, it is sufficient to determine the consistency order of a one-step procedure. In principle, this can be achieved by a Taylor expansion of to powers of . In practice, the resulting formulas for higher orders become very complicated and confusing, so that additional concepts and notations are required. ${\ displaystyle f}$${\ displaystyle \ eta (t, h)}$${\ displaystyle h}$

## Stiffness and A-stability

The convergence order of a method is an asymptotic statement that describes the behavior of the approximations when the step size converges to zero. However, it does not say anything about whether the method actually calculates a usable approximation for a given fixed step size. Charles Francis Curtiss and Joseph O. Hirschfelder first described in 1952 that this can actually be a big problem with certain types of initial value problems . They had observed that the solutions of some systems of differential equations of chemical reaction kinetics could not be calculated using explicit numerical methods, and they named them such initial value problems "stiff". There are numerous mathematical criteria for determining how stiff a given problem is. Clearly, rigid initial value problems are mostly systems of differential equations in which some components become constant very quickly, while other components change only slowly. Such behavior typically occurs when modeling chemical reactions. The most useful definition of stiffness for practical application is: An initial value problem is stiff if one would have to choose the step size “too small” in order to obtain a usable solution when solving it with explicit one-step methods. Such problems can only be solved with implicit procedures.

To calculate an exponentially falling solution (blue), the explicit Euler method (red) is completely useless if the step size is too large; the implicit Euler method (green) determines the solution qualitatively correctly for any step size.

This effect can be illustrated more precisely by examining how the individual methods cope with exponential decay . To do this, consider the test equation according to the Swedish mathematician Germund Dahlquist

${\ displaystyle y '(t) = \ lambda y (t), \ quad y (0) = 1}$

with the exponentially decreasing solution . The graphic on the right shows - as an example of the explicit and the implicit Euler method - the typical behavior of these two groups of methods in this seemingly simple initial value problem: If one uses too large a step size in an explicit method, then strongly oscillating values ​​result, which are in the Build up the calculation and move further and further away from the exact solution. Implicit methods, on the other hand, typically calculate the solution qualitatively correctly for any step sizes, namely as an exponentially decreasing sequence of approximate values. ${\ displaystyle \ lambda <0}$${\ displaystyle y (t) = e ^ {\ lambda t}}$

The above test equation is also considered somewhat more generally for complex values ​​of . In this case, the solutions oscillations whose amplitude remains bounded if if applicable, ie the real component of less than or equal to 0. This allows a desirable property of one-step methods to be formulated, which are to be used for stiff initial value problems: the so-called A-stability. A method is called A-stable if it is applied to the test equation for any step size and calculated for all with a sequence of approximations which (like the true solution) remains constrained. The implicit Euler method and the implicit trapezoid method are the simplest examples of A-stable one-step methods. On the other hand, it can be shown that an explicit procedure can never be A-stable. ${\ displaystyle \ lambda}$${\ displaystyle \ operatorname {Re} (\ lambda) \ leq 0}$${\ displaystyle \ lambda}$${\ displaystyle h> 0}$${\ displaystyle \ lambda}$${\ displaystyle \ operatorname {Re} (\ lambda) \ leq 0}$${\ displaystyle y_ {0}, y_ {1}, y_ {2}, \ dotsc}$

## Special processes and process classes

Some one-step procedures in comparison

### Simple procedures of order 1 and 2

As the French mathematician Augustin-Louis Cauchy proved around 1820, the Euler method has the order of convergence 1. If the slopes of the explicit Euler method and the implicit Euler method , as they exist at the two end points of a step, can be averaged hope to get a better approximation over the whole interval. In fact, it can be proven that the implicit trapezoidal method thus obtained${\ displaystyle f (t_ {j}, y_ {j})}$${\ displaystyle f (t_ {j + 1}, y_ {j + 1})}$

${\ displaystyle y_ {j + 1} = y_ {j} + {\ frac {h} {2}} {\ Big (} f (t_ {j}, y_ {j}) + f (t_ {j + 1 }, y_ {j + 1}) {\ Big)}}$

the convergence order has 2. This method has very good stability properties, but is implicit, so that an equation for has to be solved in each step . If this quantity is approximated on the right-hand side of the equation using the explicit Euler method, Heun's explicit method is created${\ displaystyle y_ {j + 1}}$

${\ displaystyle y_ {j + 1} = y_ {j} + {\ frac {h} {2}} {\ Big (} f (t_ {j}, y_ {j}) + f {\ big (} t_ {j + 1}, y_ {j} + hf (t_ {j}, y_ {j}) {\ big)} {\ Big)}}$,

which also has the convergence order 2. Another simple explicit method of order 2, the improved Euler method , can be obtained by considering the following: A “mean” slope in the process step would be the slope of the solution in the middle of the step, i.e. at the point . However, since the solution is unknown, it is approximated by an explicit Euler step with half the step size. The procedural rule results ${\ displaystyle y}$${\ displaystyle t_ {j} + {\ tfrac {h} {2}}}$

${\ displaystyle y_ {j + 1} = y_ {j} + hf {\ big (} t_ {j} + {\ tfrac {h} {2}}, y_ {j} + {\ tfrac {h} {2 }} f (t_ {j}, y_ {j}) {\ big)}}$.

These one-step methods of order 2 were all published in 1895 by the German mathematician Carl Runge as improvements to the Euler method .

### Runge-Kutta method

The classic fourth-order Runge-Kutta method averages four auxiliary slopes (red) in each step

The mentioned ideas for simple one-step methods lead, if they are further generalized, to the important class of Runge-Kutta methods. For example, Heun's method can be presented more clearly as follows: First an auxiliary slope is calculated, namely the slope of the explicit Euler method. This determines another auxiliary slope, here . The method gradient actually used then results as a weighted mean of the auxiliary slopes, in the Heun method . This procedure can be generalized to more than two auxiliary slopes. A -step Runge-Kutta method first calculates auxiliary slopes by evaluating them at suitable points and then as a weighted mean. In the case of an explicit Runge-Kutta method, the auxiliary slopes are calculated directly one after the other; in the case of an implicit method, they result as solutions to a system of equations. A typical example is the explicit classic Runge-Kutta method of order 4, which is sometimes simply referred to as the Runge-Kutta method: First, the four auxiliary slopes are used ${\ displaystyle k_ {1} = f (t_ {j}, y_ {j})}$${\ displaystyle k_ {2} = f (t_ {j} + h, y_ {j} + hk_ {1})}$${\ displaystyle \ Phi}$${\ displaystyle {\ tfrac {1} {2}} k_ {1} + {\ tfrac {1} {2}} k_ {2}}$${\ displaystyle s}$${\ displaystyle k_ {1}, \ dotsc, k_ {s}}$${\ displaystyle f}$${\ displaystyle \ Phi}$${\ displaystyle k_ {1}, k_ {2}, k_ {3}, \ dotsc}$

{\ displaystyle {\ begin {aligned} k_ {1} & = f (t_ {j}, y_ {j}) \\ k_ {2} & = f (t_ {j} + {\ tfrac {h} {2 }}, y_ {j} + {\ tfrac {h} {2}} k_ {1}) \\ k_ {3} & = f (t_ {j} + {\ tfrac {h} {2}}, y_ {j} + {\ tfrac {h} {2}} k_ {2}) \\ k_ {4} & = f (t_ {j} + h, y_ {j} + hk_ {3}) \ end {aligned }}}

and then the weighted average as the process slope

${\ displaystyle {\ tfrac {1} {6}} k_ {1} + {\ tfrac {1} {3}} k_ {2} + {\ tfrac {1} {3}} k_ {3} + {\ tfrac {1} {6}} k_ {4}}$

used. The German mathematician Wilhelm Kutta published this well-known method in 1901, after Karl Heun had found a three-step, one-step method of order 3 a year earlier .

The construction of explicit processes of even higher order with the smallest possible number of stages is a mathematically very demanding problem. As John C. Butcher was able to show in 1965, there are, for example, only a minimum of six-step procedures for order 5; an explicit Runge-Kutta method of the 8th order requires at least 11 steps. In 1978 the Austrian mathematician Ernst Hairer found a procedure of order 10 with 17 steps. The coefficients for such a procedure must satisfy 1205 determining equations. With implicit Runge-Kutta methods, the situation is simpler and clearer: for every number of stages there is a method of order ; this is at the same time the maximum attainable order. ${\ displaystyle s}$${\ displaystyle p = 2s}$

### Extrapolation method

The idea of extrapolation is not limited to solving initial value problems with one-step methods, but can be applied analogously to all numerical methods that discretize the problem to be solved with a step size . A well-known example of an extrapolation method is the Romberg integration for the numerical calculation of integrals. In general, let therefore be a value that is to be determined numerically, in the case of this article the value of the solution function of an initial value problem at a given point. A numerical method, for example a one-step method, calculates an approximate value for this , which depends on the choice of the step size . Here, it is assumed that the method is convergent, ie that against if converges converges to zero. This convergence, however, is only a purely theoretical statement, since in real application of the method approximate values ​​can be calculated for a finite number of different step sizes, but of course you cannot let the step size "converge towards zero". The calculated approximations for different step sizes can, however, be understood as information about the (unknown) function : In the extrapolation method , an interpolation polynomial is used as an approximation, i.e. a polynomial with ${\ displaystyle h}$${\ displaystyle v}$${\ displaystyle {\ tilde {v}} (h)}$${\ displaystyle h> 0}$${\ displaystyle {\ tilde {v}} (h)}$${\ displaystyle v}$${\ displaystyle h}$${\ displaystyle {\ tilde {v}} (h_ {1}), {\ tilde {v}} (h_ {2}), \ dotsc, {\ tilde {v}} (h_ {m})}$${\ displaystyle h_ {1}> h_ {2}> \ ldots> h_ {m}}$${\ displaystyle {\ tilde {v}}}$${\ displaystyle {\ tilde {v}}}$${\ displaystyle P}$

${\ displaystyle P (h_ {k}) = {\ tilde {v}} (h_ {k})}$

for . The value of the polynomial at the point is then used as a calculable approximation for the non-calculable limit value of for towards zero. Roland Bulirsch and Josef Stoer published an early successful extrapolation algorithm for initial value problems in 1966. ${\ displaystyle k = 1,2, \ dotsc, m}$${\ displaystyle P (0)}$${\ displaystyle h = 0}$${\ displaystyle {\ tilde {v}} (h)}$${\ displaystyle h}$

Extrapolation on in a procedure of order${\ displaystyle h = 0}$${\ displaystyle p = 2}$

A concrete example in the case of a one-step method of order can make the general procedure of extrapolation understandable. With such a method, the calculated approximation for small step sizes can be easily converted into a polynomial of the form ${\ displaystyle p}$${\ displaystyle h}$

${\ displaystyle P (h) = a + bh ^ {p}}$

with initially unknown parameters and approximate. Now is calculated by the method for a pitch and half pitch of two approximations and are obtained from the interpolation conditions , and two linear equations for the unknowns and . The value extrapolated to ${\ displaystyle a}$${\ displaystyle b}$${\ displaystyle h_ {1}}$${\ displaystyle h_ {2} = {\ tfrac {1} {2}} h_ {1}}$${\ displaystyle y_ {h_ {1}}}$${\ displaystyle y_ {h_ {2}}}$${\ displaystyle P (h_ {1}) = y_ {h_ {1}}}$${\ displaystyle P (h_ {2}) = y_ {h_ {2}}}$${\ displaystyle a}$${\ displaystyle b}$${\ displaystyle h = 0}$

${\ displaystyle P (0) = a = y_ {h_ {2}} + {\ frac {y_ {h_ {2}} - y_ {h_ {1}}} {2 ^ {p} -1}}}$

then generally represents a much better approximation than the two values ​​initially calculated. It can be shown that the order of the one-step process obtained in this way is at least , i.e. at least 1 greater than in the original process. ${\ displaystyle p + 1}$

### Procedure with step size control

One advantage of the one-step method is that any step size can be used in each step independently of the other steps . In practice, the question obviously arises as to how to choose . In real applications there will always be an error tolerance with which the solution of an initial value problem is to be calculated; For example, it would be pointless to determine a numerical approximation that is significantly “more precise” than the data with measurement errors for starting values ​​and parameters of the given problem. The goal will therefore be to choose the step sizes so that, on the one hand, the specified error tolerances are adhered to, but on the other hand, to use as few steps as possible in order to keep the computing effort small. In general, this can only be achieved if the step sizes are adapted to the course of the solution: small steps where the solution changes significantly, large steps where it is almost constant. ${\ displaystyle j}$${\ displaystyle h_ {j}}$ ${\ displaystyle h_ {j}}$

In the case of well-conditioned initial value problems, it can be shown that the global procedural error is approximately equal to the sum of the local truncation errors in the individual steps. Therefore, as large a step size as possible should be selected for which is below a selected tolerance threshold. The problem here is that it cannot be calculated directly, since it depends on the unknown exact solution of the initial value problem at the point . The basic idea of ​​step size control is therefore to approximate with a method that is more precise than the basic method on which it is based. ${\ displaystyle \ eta _ {j}: = \ | \ eta (t_ {j}, h_ {j}) \ |}$${\ displaystyle h_ {j}}$${\ displaystyle \ eta _ {j}}$${\ displaystyle \ eta _ {j}}$${\ displaystyle y (t_ {j})}$${\ displaystyle t_ {j}}$${\ displaystyle y (t_ {j})}$

Two basic ideas for step size control are the step size halving and the embedded procedures . When halving the step size, the result for two steps with half the step size is calculated as a comparison value in addition to the actual process step. A more precise approximation for is then determined from both values ​​by extrapolation , and the local error is thus estimated. If this is too large, this step is discarded and repeated with a smaller step size. If it is significantly smaller than the specified tolerance, the step size can be increased in the next step. The additional computational effort for this step size halving method is relatively large; therefore, modern implementations mostly use so-called embedded methods for step size control. The basic idea is to calculate two approximations for each step using two one-step methods that have different orders of convergence, and thus to estimate the local error. In order to optimize the computing effort, the two procedures should have as many computing steps as possible in common: They should be “embedded in one another”. For example, embedded Runge-Kutta methods use the same auxiliary slopes and only differ in how they average them. Well-known embedded methods include the Runge-Kutta-Fehlberg methods ( Erwin Fehlberg , 1969) and the Dormand-Prince methods (JR Dormand and PJ Prince, 1980). ${\ displaystyle y (t_ {j})}$${\ displaystyle \ eta _ {j}}$${\ displaystyle y (t_ {j})}$

## Practical example: Solving initial value problems with numerical software

Numerous software implementations have been developed for the mathematical concepts presented in an overview in this article, which enable the user to numerically solve practical problems in a simple manner. As a concrete example, a solution of the Lotka-Volterra equations is to be calculated with the widely used numerical software Matlab . The Lotka-Volterra equations are a simple model from biology that describes the interactions between predator and prey populations . For this purpose the differential equation system is given

{\ displaystyle {\ begin {aligned} y_ {1} '(t) & = ay_ {1} (t) -by_ {1} (t) y_ {2} (t) \\ y_ {2}' (t ) & = cy_ {1} (t) y_ {2} (t) -dy_ {2} (t) \ end {aligned}}}

with the parameters and the initial condition , . Here and correspond to the temporal development of the prey or predator population. The solution should be calculated on the time interval . ${\ displaystyle a = 1, b = 2, c = 1, d = 1}$${\ displaystyle y_ {1} (0) = 3}$${\ displaystyle y_ {2} (0) = 1}$${\ displaystyle y_ {1}}$${\ displaystyle y_ {2}}$${\ displaystyle [0.20]}$

For the calculation using Matlab, the function on the right-hand side of the differential equation is first defined for the given parameter values : ${\ displaystyle f}$${\ displaystyle y '= f (t, y)}$

a = 1; b = 2; c = 1; d = 1;
f = @(t,y) [a*y(1) - b*y(1)*y(2); c*y(1)*y(2) - d*y(2)];


In addition, the time interval and the initial values ​​are required:

t_int = [0, 20];
y0 = [3; 1];


Then the solution can be calculated:

[t, y] = ode45(f, t_int, y0);


The Matlab function ode45implements a one-step method that uses two embedded explicit Runge-Kutta methods with orders of convergence 4 and 5 for step size control.

The solution can now be drawn, as a blue curve and as a red one; the calculated points are marked by small circles: ${\ displaystyle y_ {1}}$${\ displaystyle y_ {2}}$

figure(1)
plot(t, y(:,1), 'b-o', t, y(:,2), 'r-o')


The result is shown below in the picture on the left. The right picture shows the step sizes used by the process and was created with

figure(2)
plot(t(1:end-1), diff(t))

Solutions to the Lotka-Volterra equations
Step sizes used

This example can also be executed with the free numerical software GNU Octave without changes . With the method implemented there, however, a slightly different step size sequence results.

## literature

• John C. Butcher : Numerical Methods for Ordinary Differential Equations . John Wiley & Sons, Chichester 2008, ISBN 978-0-470-72335-7 .
• Wolfgang Dahmen , Arnold Reusken: Numerics for engineers and natural scientists . 2nd Edition. Springer, Berlin / Heidelberg 2008, ISBN 978-3-540-76492-2 , chap. 11: Ordinary differential equations .
• Peter Deuflhard , Folkmar Bornemann : Numerical Mathematics 2 - Ordinary Differential Equations . 3. Edition. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1 .
• David F. Griffiths, Desmond J. Higham: Numerical Methods for Ordinary Differential Equations - Initial Value Problems . Springer, London 2010, ISBN 978-0-85729-147-9 .
• Robert Plato: Numerical Mathematics compact . 4th edition. Vieweg + Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2 , chap. 7: One-step procedure for initial value problems .
• Hans-Jürgen Reinhardt: Numerics of ordinary differential equations . 2nd Edition. Walter de Gruyter, Berlin / Boston 2012, ISBN 978-3-11-028045-6 .
• Hans Rudolf Schwarz, Norbert Köckler: Numerical Mathematics . 8th edition. Vieweg + Teubner, Wiesbaden 2011, ISBN 978-3-8348-1551-4 , chap. 8: Initial value problems .
• Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerics of ordinary differential equations . 2nd Edition. Springer Spectrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8 .

## Individual evidence

1. ^ Thomas Sonar : 3000 Years of Analysis . Springer, Berlin / Heidelberg 2011, ISBN 978-3-642-17203-8 , pp. 378-388 and 401-426 .
2. Jean-Luc Chabert et al. a .: A History of Algorithms . Springer, Berlin / Heidelberg 1999, ISBN 978-3-540-63369-3 , p. 374-378 .
3. Wolfgang Dahmen, Arnold Reusken: Numerics for engineers and natural scientists . 2nd Edition. Springer, Berlin / Heidelberg 2008, ISBN 978-3-540-76492-2 , p. 386 f .
4. Wolfgang Dahmen, Arnold Reusken: Numerics for engineers and natural scientists . 2nd Edition. Springer, Berlin / Heidelberg 2008, ISBN 978-3-540-76492-2 , p. 386-392 .
5. ^ Hans Rudolf Schwarz, Norbert Köckler: Numerical Mathematics . 8th edition. Vieweg + Teubner, Wiesbaden 2011, ISBN 978-3-8348-1551-4 , pp. 350 f .
6. ^ Robert Plato: Numerical Mathematics compact . 4th edition. Vieweg + Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2 , pp. 157 .
7. ^ Robert Plato: Numerical Mathematics compact . 4th edition. Vieweg + Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2 , pp. 156 .
8. ^ Robert Plato: Numerical Mathematics compact . 4th edition. Vieweg + Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2 , pp. 157 .
9. Hans-Jürgen Reinhardt: Numerics of ordinary differential equations . 2nd Edition. Walter de Gruyter, Berlin / Boston 2012, ISBN 978-3-11-028045-6 , p. 42 f .
10. ^ John C. Butcher: Numerical Methods for Ordinary Differential Equations . John Wiley & Sons, Chichester 2008, ISBN 978-0-470-72335-7 , pp. 95-100 .
11. ^ JC Butcher: Numerical methods for ordinary differential equations in the 20th century . In: Journal of Computational and Applied Mathematics . tape 125 , no. 1–2 , December 15, 2000, pp. 21st f . ( online ).
12. ^ Peter Deuflhard, Folkmar Bornemann: Numerical Mathematics 2 - Ordinary Differential Equations . 3. Edition. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1 , p. 228 f .
13. ^ Peter Deuflhard, Folkmar Bornemann: Numerical Mathematics 2 - Ordinary Differential Equations . 3. Edition. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1 , p. 229-231 .
14. Wolfgang Dahmen, Arnold Reusken: Numerics for engineers and natural scientists . 2nd Edition. Springer, Berlin / Heidelberg 2008, ISBN 978-3-540-76492-2 , p. 443 f .
15. Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerics of ordinary differential equations . 2nd Edition. Springer Spectrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8 , p. 258 f .
16. Jean-Luc Chabert et al. a .: A History of Algorithms . Springer, Berlin / Heidelberg 1999, ISBN 978-3-540-63369-3 , p. 378 f .
17. Jean-Luc Chabert et al. a .: A History of Algorithms . Springer, Berlin / Heidelberg 1999, ISBN 978-3-540-63369-3 , p. 381-388 .
18. Wolfgang Dahmen, Arnold Reusken: Numerics for engineers and natural scientists . 2nd Edition. Springer, Berlin / Heidelberg 2008, ISBN 978-3-540-76492-2 , p. 406 f .
19. ^ JC Butcher: Numerical methods for ordinary differential equations in the 20th century . In: Journal of Computational and Applied Mathematics . tape 125 , no. 1–2 , December 15, 2000, pp. 4-6 ( online ).
20. ^ Peter Deuflhard, Folkmar Bornemann: Numerical Mathematics 2 - Ordinary Differential Equations . 3. Edition. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1 , p. 160-162 .
21. Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerics of ordinary differential equations . 2nd Edition. Springer Spectrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8 , p. 219-221 .
22. Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerics of ordinary differential equations . 2nd Edition. Springer Spectrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8 , p. 79 ff .
23. ^ JC Butcher: Numerical methods for ordinary differential equations in the 20th century . In: Journal of Computational and Applied Mathematics . tape 125 , no. 1–2 , December 15, 2000, pp. 26 ( online ).
24. ^ Robert Plato: Numerical Mathematics compact . 4th edition. Vieweg + Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2 , pp. 171-173 .
25. Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerics of ordinary differential equations . 2nd Edition. Springer Spectrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8 , p. 57-59 .
26. ^ Peter Deuflhard, Folkmar Bornemann: Numerical Mathematics 2 - Ordinary Differential Equations . 3. Edition. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1 , p. 199-204 .
27. ^ Robert Plato: Numerical Mathematics compact . 4th edition. Vieweg + Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2 , chap. 7: One-step method for initial value problems , p. 173-177 .
28. Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerics of ordinary differential equations . 2nd Edition. Springer Spectrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8 , p. 64-70 .
29. ode45: Solve nonstiff differential equations - medium order method. MathWorks, accessed November 23, 2017 .
 This article was added to the list of excellent articles on December 11, 2017 in this version .