One-step method

In numerical mathematics, one-step methods and multi-step methods are a large group of calculation methods for solving initial value problems. This problem, in which an ordinary differential equation is given together with an initial condition, plays a central role in all natural and engineering sciences and is also becoming increasingly important in the economic and social sciences, for example. Initial value problems are used to analyze, simulate or predict dynamic processes. The basic idea behind one-step methods is that they calculate approximation points step by step along the desired solution, starting from the given starting point. They only use the most recently determined approximation for the next step, in contrast to multi-step methods, which also include points 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, which require an equation 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 methods was presented in 1883, Carl Runge, Karl Heun and Wilhelm Kutta developed significant improvements to Euler's method around 1900.…

In numerical mathematics, one-step methods and multi-step methods are a large group of calculation methods for solving initial value problems. This problem, in which an ordinary differential equation is given together with an initial condition, plays a central role in all natural and engineering sciences and is also becoming increasingly important in the economic and social sciences, for example. Initial value problems are used to analyze, simulate or predict dynamic processes. The basic idea behind one-step methods is that they calculate approximation points step by step along the desired solution, starting from the given starting point. They only use the most recently determined approximation for the next step, in contrast to multi-step methods, which also include points 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, which require an equation 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 methods was presented in 1883, Carl Runge, Karl Heun and Wilhelm Kutta developed significant improvements to Euler's method around 1900. These gave rise to the large group of Runge-Kutta methods, which form the most important class of one-step methods. Further developments in the 20th century include the idea of extrapolation and, above all, considerations on step width control, i.e. the selection of suitable lengths for the individual steps of a method. These concepts form the basis for solving difficult initial value problems, as they occur in modern applications, efficiently and with the required accuracy using computer programs. For example, when a body moves, its average speed is simply the distance traveled divided by the time required to travel it. However, in order to mathematically formulate the instantaneous velocity v(t) of the body at a certain point in time t , a limit transition is necessary: Consider short time spans of length \Delta t , the distances traveled \Delta x and the corresponding average velocities \tfrac{\Delta x}{\Delta t} .If the time period Δ 𝑡 is now allowed to converge towards zero and if the average velocities also approach a fixed value, then this value is called the (instantaneous) velocity v(t) at the given time t . If \Delta t denotes the position of the body at time 𝑡 , then write v(t) = x'(t) and call v the derivative of x . The decisive step in the direction of differential equation models is now the reverse question: In the example of the moving body, let the velocity v(t) be known at every point in time 𝑡 and its position x(t) be determined from this. It is clear that the initial position of the body at a point in time 𝑡 0 must also be known in order to be able to solve this problem unambiguously. We are therefore looking for a function x(t) with x'(t) = v(t) that fulfills the initial condition x'(t) = v(t) with given values t_0 and x_0 . In the example of determining the position 𝑥 of a body from its velocity, the derivative of the function being searched for is explicitly given. In most cases, however, the important general case of ordinary differential equations exists for a sought-after variable y : Due to the laws of nature or the model assumptions, a functional relation is known that specifies how the deriativey y'(t) of the function to be determined can be calculated from t and from the (unknown) value y(t) . In addition, an initial condition must again be given, which can be obtained, for example, from a measurement of the required variable at a fixed point in time. To summarize, the following general type of task exists: Find the function y that satisfies the equations y'(t) = f(t, y(t)), \quad y(t_0) = y_0 is fulfilled, where f is a given function. A simple example is a variable y that grows exponentially. This means that the instantaneous change, i.e. the derivative y'(t) , is proportional to y(t) itself. Therefore, y'(t) = \lambda y(t) with a growth rate \lambda and, for example, an initial condition y(0) = y_0 . In this case, the required solution 𝑦 can already be found using elementary differential calculus and specified using the exponential function: y(t) = y_0 e^{\lambda t} . The required function y in a differential equation can be vector-valued, i.e. for each t , y(t) = (y_1(t), \dotsc, y_d(t)) can be a vector with d components. This is also referred to as an d -dimensional system of differential equations. In the case of a moving body, y(t) is its position in d -dimensional Euclidean space and y'(t) is its velocity at time t . The differential equation therefore specifies the velocity of the trajectory with direction and magnitude at each point in time and space. The trajectory itself is to be calculated from this. The simplest and most basic one-step method is the explicit Euler method, which was introduced by the Swiss mathematician and physicist Leonhard Euler in 1768 in his textbook Institutiones Calculi Integralis. The idea of this method is to approximate the solution sought by a piecewise linear function in which the gradient of the straight line piece is given by t_j in each step from the point t_{j+1} to the point f(t_j, y_j) . In more detail: The problem definition already gives a value of the function being searched for, namely y(t_0) = y_0 . However, the derivative at this point is also known, as y'(t_0) = f(t_0, y_0) applies. This allows the tangent to the graph of the solution function to be determined and used as an approximation. At the point t_1 > t_0 the following results with the step size h_0 := t_1 - t_0 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 following calculation rule for the explicit Euler method y_{j+1} = y_j + h_j f(t_j, y_j), \quad j = 0,1,2,\dotsc with the increments h_j = t_{j+1} - t_j . The explicit Euler method is the starting point for numerous generalizations in which the gradient f(t_j, y_j) is replaced by gradients that approximate the behaviour of the solution between the points t_j and t_{j+1} more precisely. An additional idea for one-step methods is provided by the implicit Euler method, which uses f(t_{j+1}, y_{j+1}) as the gradient. At first glance, this choice does not seem very suitable, as y_{j+1} is unknown. However, as a procedural step, we now obtain the equation y_{j+1} = y_j + h_j f(t_{j+1}, y_{j+1}) from which y_{j+1} can be calculated (using a numerical method if necessary). If, for example, the arithmetic mean of the slopes of the explicit and implicit Euler method is selected as the slope, the implicit trapezoidal method is obtained. In turn, an explicit method can be obtained from this if, for example, the unknown y_{j+1} on the right-hand side of the equation is approximated using the explicit Euler method, the so-called Heun method. All these methods and all other generalizations have the basic idea of one-step methods in common: the step y_{j+1} = y_j + h_j \Phi with a gradient \Phi that can depend on t_j , y_j and h_j as well as (for implicit methods) on y_{j+1} . y'(t) = f(t, y(t)) , \quad y(t_0) = y_0 . It is assumed that the solution y \colon I \to \R^d exists on a given interval I = [t_0, T] and is uniquely determined. Are t_0 Intermediate positions in the interval I and h_j = t_{j+1} - t_j the corresponding increments, then this is given by y_{j+1} = y_j + h_j \Phi(t_j, y_j, y_{j+1}, h_j) , \quad j=0,\dotsc,n-1 given method is a one-step method with method function \Phi . If \Phi does not depend on y_{j+1} , then it is called an explicit one-step method. Otherwise, an equation for j must be solved in each step j and the method is called implicit. \max_{j=0,\dotsc,n} \|y_j - y(t_j)\| \leq C h^p applies to all sufficiently small h with a constant C > 0 that is independent of h . The order of convergence is the most important parameter for comparing different one-step methods. A method with a higher order of convergence p generally delivers a smaller total error for a given step size or, conversely, fewer steps are required to achieve a given accuracy. For a method with p = 1 , it is to be expected that the error will only be approximately halved if the step size is halved. With a method of convergence order p=4 , on the other hand, it can be assumed that the error is reduced by a factor of approximately \bigl(\tfrac{1}{2}\bigr)^4 = \tfrac{1}{16} . y_{j+1} = y_j + h \Phi(t_j, y_j, h) with a constant step size h exists. With the true solution t \mapsto y(t) , the local truncation error (also called local process error) \eta is defined as \eta(t,h) = y(t) + h \Phi(t, y(t), h) - y(t+h) . Thus, one assumes that the exact solution is known, starts a method step at the point (t,y(t)) and forms the difference to the exact solution at the point t + h . This defines: A one-step method has the consistency order p \geq 1 if the estimate \|\eta(t,h)\| \leq C h^{p+1} applies to all sufficiently small h with a constant C > 0 that is independent of h . The striking difference between the definitions of the consistency order and the convergence order is the power h^{p+1} instead of h^p . This can be clearly interpreted as meaning that a power of the step size is "lost" during the transition from local to global error. The following theorem, which is central to the theory of one-step methods, applies: If the process function \Phi is Lipschitz-continuous and the associated one-step process has the consistency order p , then it also has the convergence order p . The Lipschitz continuity of the process function as an additional requirement for stability is generally always fulfilled if the function f from the differential equation itself is Lipschitz-continuous. This requirement must be assumed for most applications anyway in order to guarantee the unambiguous solvability of the initial value problem. According to the theorem, it is therefore sufficient to determine the consistency order of a one-step method. In principle, this can be achieved by Taylor expansion of \eta(t,h) to powers of h . In practice, the resulting formulas for higher orders become very complicated and confusing, so that additional concepts and notations are required. y'(t) = \lambda y(t), \quad y(0)=1 with the exponentially decreasing solution \lambda for y(t) = e^{\lambda t} . The adjacent diagram shows - as an example for the explicit and implicit Euler method - the typical behavior of these two groups of methods for this seemingly simple initial value problem: If too large a step size is used in an explicit method, this results in strongly oscillating values that build up over the course of the calculation and move further and further away from the exact solution. Implicit methods, on the other hand, typically calculate the solution for arbitrary step sizes qualitatively correctly, namely as an exponentially decreasing sequence of approximate values. More generally, the above test equation is also considered for complex values of \lambda . In this case, the solutions are oscillations whose amplitude remains limited precisely when \operatorname{Re}(\lambda) \leq 0 , i.e. the real part of \lambda is less than or equal to 0. This makes it possible to formulate a desirable property of one-step methods that are to be used for stiff initial value problems: the so-called A-stability. A method is called A-stable if it calculates a sequence of approximations h > 0 for any step size y_0, y_1, y_2,\dotsc applied to the test equation for all \lambda with \operatorname{Re}(\lambda) \leq 0 , which remains bounded (like the true solution). The implicit Euler method and the implicit trapezoidal method are the simplest examples of A-stable one-step methods. On the other hand, it can be shown that an explicit method can never be A-stable. y_{j+1} = y_j + \frac{h}{2} \Big( f(t_j, y_j) + f(t_{j+1},y_{j+1}) \Big) has a convergence order of 2. This method has very good stability properties, but is implicit, meaning that an equation for 𝑦 𝑗 + 1 must be solved in each step. If this variable is approximated on the right-hand side of the equation using the explicit Euler method, the result is the explicit method of Heun 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 convergence order 2. Another simple explicit method of order 2, the improved Euler method, is obtained by the following consideration: A "mean" slope in the method step would be the slope of the solution 𝑦 in the middle of the step, i.e. at the point y_{j+1} . However, as the solution is unknown, it is approximated by an explicit Euler step with half the step size. This results in the following procedure y_{j+1} = y_j + h f\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 as improvements of the Euler method in 1895¹ by the German mathematician Carl Runge. The aforementioned ideas for simple one-step methods lead to the important class of Runge-Kutta methods when generalized further. For example, Heun's method can be presented more clearly as follows: First, an auxiliary slope k_1 = f(t_j, y_j) is calculated, namely the slope of the explicit Euler method. This is used to determine a further auxiliary slope, here k_2 = f(t_j + h, y_j + h k_1) . The actual process gradient \Phi used is then calculated as a weighted average of the auxiliary gradients, i.e. \tfrac{1}{2} k_1 + \tfrac{1}{2} k_2 in Heun's method. This procedure can be generalized to more than two auxiliary slopes. An s - -stage Runge-Kutta method first calculates auxiliary slopes k_1, \dotsc, k_s by evaluating 𝑓 at suitable points and then \Phi as a weighted average. In an explicit Runge-Kutta method, the auxiliary slopes k_1, k_2, k_3, \dotsc are calculated directly one after the other; in an implicit method, they are obtained as solutions to a system of equations. A typical example is the explicit classical Runge-Kutta method of order 4, which is sometimes simply referred to as the Runge-Kutta method: First, the four auxiliary slopes \begin{align} 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{align} and then the weighted average is calculated as the process slope \tfrac{1}{6} k_1 + \tfrac{1}{3} k_2 + \tfrac{1}{3} k_3 + \tfrac{1}{6} k_4 is used. This well-known method was published by the German mathematician Wilhelm Kutta in 1901¹, after Karl Heun had found a three-step one-step method of order 3 a year earlier. The construction of explicit methods of even higher order with the smallest possible number of steps is a mathematically quite demanding problem. As John C. Butcher was able to show in 1965, there are, for example, only a minimum of six steps for order 5; an explicit Runge-Kutta method of order 8 requires at least 11 steps. In 1978, the Austrian mathematician Ernst Hairer found a method of order 10 with 17 levels. The coefficients for such a method must fulfill 1205 determinant equations. With implicit Runge-Kutta methods, the situation is simpler and clearer: for every number of steps s there is a method of order p = 2s ; this is also the maximum achievable order. P(h_k) = \tilde v(h_k) for k = 1,2,\dotsc,m . The value P(0) of the polynomial at the point h=0 is then used as a computable approximation for the non-computable limit value of \tilde v(h) for h towards zero. An early successful extrapolation algorithm for initial value problems was published by Roland Bulirsch and Josef Stoer in 1966¹. A concrete example in the case of a one-step method of order p can illustrate the general procedure of extrapolation. With such a method, the calculated approximation for small step sizes ℎ can be easily described by a polynomial of the form P(h) = a + b h^p with initially unknown parameters a and b . If you now calculate two approximations y_{h_1} and y_{h_2} using the method for a step size h_1 and for half the step size h_2 = \tfrac{1}{2} h_1 , two linear equations for the unknowns a and b are obtained from the interpolation conditions P(h_1) = y_{h_1} and P(h_2) = y_{h_2} . The value extrapolated to P(0) = a = y_{h_2} + \frac{y_{h_2} - y_{h_1}}{2^p-1} is then generally a significantly better approximation than the two values calculated initially. It can be shown that the order of the one-step method obtained in this way is at least p+1 , i.e. at least 1 greater than the original method. For well-conditioned initial value problems, it can be shown that the global process error is approximately equal to the sum of the local truncation errors \eta_j := \|\eta(t_j, h_j)\| in the individual steps. Therefore, the largest possible h_j should be selected as the step size, for which \eta_j is below a selected tolerance threshold. The problem here is that \eta_j cannot be calculated directly, as it depends on the unknown exact solution y(t_j) of the initial value problem at the point t_j . The basic idea of step size control is therefore to approximate y(t_j) with a method that is more accurate than the underlying basic method. Two basic ideas for step width control are step width halving and embedded processes. With step size halving, 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 y(t_j) is then determined from both values by extrapolation and the local error 𝜂 𝑗 is 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 width halving procedure is relatively high; this is why modern implementations usually use so-called embedded procedures for step width control. The basic idea is to calculate two approximations for y(t_j) in each step using two one-step methods that have different orders of convergence and thus estimate the local error. In order to optimize the computational effort, the two methods should have as many computational steps in common as possible: They should be "embedded in each other". Embedded Runge-Kutta methods, for example, use the same auxiliary slopes and differ only in how they average them. Well-known embedded methods include the Runge-Kutta-Fehlberg method ( , 1969) and the Dormand-Prince method (J. R. Dormand and P. J. Prince, 1980). y_1'(t) &= a y_1(t) - b y_1(t) y_2(t)\\ y_2'(t) &= c y_1(t) y_2(t) - d y_2(t) \end{align} with the parameters a = 1, b = 2, c = 1, d = 1 and the initial condition y_1(0) = 3 , y_2(0) = 1 . Here, y_1 and y_2 correspond to the temporal development of the prey and predator population respectively. The solution should be calculated on the time interval [0, 20] . For the calculation using Matlab, the function f is first defined for the given parameter values on the right-hand side of the differential equation 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)]; The time interval and the initial values are also required: t_int = [0, 20]; y0 = [3; 1]; The solution can then be calculated: [t, y] = ode45(f, t_int, y0); The Matlab function ode45 implements a one-step method that uses two embedded explicit Runge-Kutta methods with convergence orders 4 and 5 for step size control. The solution can now be plotted, y_1 as a blue curve and y_2 as a red curve; the calculated points are marked by small circles: figure(1) plot(t, y(:,1), 'b-o', t, y(:,2), 'r-o') The result is shown below in the left-hand image. The right-hand image shows the step sizes used by the method and was generated with figure(2) plot(t(1:end-1), diff(t)) This example can also be executed without changes using the free numerical software GNU Octave. However, the method implemented there results in a slightly different step size sequence.

Text derived from https://en.wikipedia.org/w/index.php?oldid=1366141094 under CC BY-SA 4.0; modifications have been made.