Numerical solution of ordinary differential equations. Numerical solution of differential equations (1) Order of the numerical differential equation method

To solve differential equations, it is necessary to know the value of the dependent variable and its derivatives for certain values ​​of the independent variable. If additional conditions are specified for one value of the unknown, i.e. independent variable., then such a problem is called the Cauchy problem. If the initial conditions are specified for two or more values ​​of the independent variable, then the problem is called a boundary value problem. When solving differential equations of various types, the function whose values ​​need to be determined is calculated in the form of a table.

Classification of numerical methods for solving differentials. Lv. Types.

Cauchy problem – one-step: Euler methods, Runge-Kutta methods; – multi-step: Main method, Adams method. Boundary problem – a method of reducing a boundary problem to the Cauchy problem; – finite difference method.

When solving the Cauchy problem, dif. must be specified. ur. order n or system of dif. ur. first order of n equations and n additional conditions for its solution. Additional conditions must be specified for the same value of the independent variable. When solving a boundary problem, equations must be specified. nth order or a system of n equations and n additional conditions for two or more values ​​of the independent variable. When solving the Cauchy problem, the required function is determined discretely in the form of a table with a certain specified step . When determining each successive value, you can use information about one previous point. In this case, the methods are called one-step, or you can use information about several previous points - multi-step methods.

Ordinary differential equations. Cauchy problem. One-step methods. Euler's method.

Given: g(x,y)y+h(x,y)=0, y=-h(x,y)/g(x,y)= f(x,y), x 0 , y( x 0)=y 0 . It is known: f(x,y), x 0 , y 0 . Determine the discrete solution: x i , y i , i=0,1,…,n. Euler's method is based on the expansion of a function into a Taylor series in the vicinity of the point x 0 . The neighborhood is described by step h. y(x 0 +h)y(x 0)+hy(x 0)+…+ (1). Euler's method takes into account only two terms of the Taylor series. Let us introduce some notation. Euler's formula will take the form: y i+1 =y i +y i, y i =hy(x i)=hf(x i,y i), y i+1 =y i +hf(x i,y i) (2), i= 0,1,2…, x i+1 =x i +h

Formula (2) is the formula of the simple Euler method.

Geometric interpretation of Euler's formula

To obtain a numerical solution, the tangent line passing through the equation is used. tangent: y=y(x 0)+y(x 0)(x-x 0), x=x 1,

y 1 =y(x 0)+f(x 0 ,y 0)  (x-x 0), because

x-x 0 =h, then y 1 =y 0 +hf(x 0 ,y 0), f(x 0 ,y 0)=tg £.

Modified Euler method

Given: y=f(x,y), y(x 0)=y 0 . It is known: f(x,y), x 0 , y 0 . Determine: the dependence of y on x in the form of a tabular discrete function: x i, y i, i=0.1,…,n.

Geometric interpretation

1) calculate the tangent of the angle of inclination at the starting point

tg £=y(x n ,y n)=f(x n ,y n)

2) Calculate the value  y n+1 on

end of step according to Euler's formula

 y n+1 =y n +f(x n ,y n) 3) Calculate the tangent of the angle of inclination

tangent at n+1 point: tg £=y(x n+1 ,  y n+1)=f(x n+1 ,  y n+1) 4) Calculate the arithmetic mean of the angles

tilt: tg £=½. 5) Using the tangent of the slope angle, we recalculate the value of the function at n+1 point: y n+1 =y n +htg £= y n +½h=y n +½h – formula of the modified Euler method. It can be shown that the resulting f-la corresponds to the expansion of the f-i in a Taylor series, including terms (up to h 2). The modified Eilnra method, unlike the simple one, is a method of second order accuracy, because the error is proportional to h 2.

We consider only the solution to the Cauchy problem. A system of differential equations or one equation must be converted to the form

Where ,
n-dimensional vectors; y– unknown vector function; x– independent argument,
. In particular, if n= 1, then the system turns into one differential equation. The initial conditions are set as follows:
, Where
.

If
in the vicinity of a point
is continuous and has continuous partial derivatives with respect to y, then the existence and uniqueness theorem guarantees that there is only one continuous vector function
, defined in some neighborhood of a point , satisfying equation (7) and the condition
.

Let us pay attention to the fact that the neighborhood of the point , where the solution is determined, can be very small. When approaching the boundary of this neighborhood, the solution can go to infinity, oscillate with an infinitely increasing frequency, in general, behave so badly that it cannot be continued beyond the boundary of the neighborhood. Accordingly, such a solution cannot be tracked by numerical methods on a larger segment, if one is specified in the problem statement.

Solving the Cauchy problem on [ a; b] is a function. In numerical methods, the function is replaced by a table (Table 1).

Table 1

Here
,
. The distance between adjacent table nodes is usually taken to be constant:
,
.

There are tables with variable steps. The table step is determined by the requirements of the engineering problem and not connected with the accuracy of finding a solution.

If y is a vector, then the table of solution values ​​will take the form of a table. 2.

Table 2

In the MATHCAD system, a matrix is ​​used instead of a table, and it is transposed with respect to the specified table.

Solve the Cauchy problem with accuracy ε means to get the values ​​in the specified table (numbers or vectors),
, such that
, Where
- exact solution. It is possible that the solution to the segment specified in the problem does not continue. Then you need to answer that the problem cannot be solved on the entire segment, and you need to get a solution on the segment where it exists, making this segment as large as possible.

It should be remembered that the exact solution
we don’t know (otherwise why use the numerical method?). Grade
must be justified on some other basis. As a rule, it is not possible to obtain a 100% guarantee that the assessment is being carried out. Therefore, algorithms are used to estimate the value
, which prove effective in most engineering problems.

The general principle for solving the Cauchy problem is as follows. Line segment [ a; b] is divided into a number of segments by integration nodes. Number of nodes k does not have to match the number of nodes m final table of decision values ​​(Tables 1, 2). Usually, k > m. For simplicity, we will assume that the distance between nodes is constant,
;h called the integration step. Then, according to certain algorithms, knowing the values at i < s, calculate the value . The smaller the step h, the lower the value will differ from the value of the exact solution
. Step h in this partition is already determined not by the requirements of the engineering problem, but by the required accuracy of solving the Cauchy problem. In addition, it must be selected so that at one step the table. 1, 2 fit an integer number of steps h. In this case the values y, obtained as a result of calculations with steps h at points
, are used accordingly in the table. 1 or 2.

The simplest algorithm for solving the Cauchy problem for equation (7) is the Euler method. The calculation formula is:

(8)

Let's see how the accuracy of the solution found is assessed. Let's pretend that
is the exact solution of the Cauchy problem, and also that
, although this is almost always not the case. Then where is the constant C depends on function
in the vicinity of a point
. Thus, at one step of integration (finding a solution) we get an error of the order . Because steps have to be taken
, then it is natural to expect that the total error at the last point
everything will be fine
, i.e. order h. Therefore, Euler's method is called the first order method, i.e. the error has the order of the first power of the step h. In fact, at one step of integration the following estimate can be justified. Let
– exact solution of the Cauchy problem with the initial condition
. It's clear that
does not coincide with the required exact solution
the original Cauchy problem of equation (7). However, at small h and "good" function
these two exact solutions will differ little. The Taylor remainder formula ensures that
, this gives the integration step error. The final error consists not only of errors at each integration step, but also of deviations of the desired exact solution
from exact solutions
,
, and these deviations can become very large. However, the final estimate of the error in the Euler method for a “good” function
still looks like
,
.

When applying Euler's method, the calculation proceeds as follows. According to specified accuracy ε determine the approximate step
. Determining the number of steps
and again approximately select the step
. Then again we adjust it downward so that at each step the table. 1 or 2 fit an integer number of integration steps. We get a step h. According to formula (8), knowing And , we find. By found value And
we find so on.

The resulting result may not, and generally will not, have the desired accuracy. Therefore, we reduce the step by half and again apply the Euler method. We compare the results of the first application of the method and the second in identical points . If all discrepancies are less than the specified accuracy, then the last calculation result can be considered the answer to the problem. If not, then we reduce the step by half again and apply Euler’s method again. Now we compare the results of the last and penultimate application of the method, etc.

Euler's method is used relatively rarely due to the fact that to achieve a given accuracy ε a large number of steps are required, in the order of
. However, if
has discontinuities or discontinuous derivatives, then higher order methods will produce the same error as Euler's method. That is, the same amount of calculations will be required as in the Euler method.

Of the higher order methods, the fourth order Runge–Kutta method is most often used. In it, calculations are carried out according to the formulas

This method, in the presence of continuous fourth derivatives of the function
gives an error on one step of the order , i.e. in the notation introduced above,
. In general, on the integration interval, provided that the exact solution is determined on this interval, the integration error will be of the order of .

The selection of the integration step occurs in the same way as described in Euler’s method, except that the initial approximate value of the step is selected from the relation
, i.e.
.

Most programs used to solve differential equations use automatic step selection. The gist of it is this. Let the value already be calculated . The value is calculated
in increments h, chosen during calculation . Then two integration steps are performed with step , i.e. extra node is added
in the middle between the nodes And
. Two values ​​are calculated
And
in nodes
And
. The value is calculated
, Where p– method order. If δ is less than the accuracy specified by the user, then it is assumed
. If not, then choose a new step h equal and repeat the accuracy check. If during the first check δ is much less than the specified accuracy, then an attempt is made to increase the step. For this purpose it is calculated
at the node
in increments h from node
and is calculated
in steps of 2 h from node . The value is calculated
. If is less than the specified accuracy, then step 2 h considered acceptable. In this case, a new step is assigned
,
,
. If more accuracy, then the step is left the same.

It should be taken into account that programs with automatic selection of the integration step achieve the specified accuracy only when performing one step. This occurs due to the accuracy of the approximation of the solution passing through the point
, i.e. approximation of the solution
. Such programs do not take into account how much the solution
differs from the desired solution
. Therefore, there is no guarantee that the specified accuracy will be achieved throughout the entire integration interval.

The described Euler and Runge–Kutta methods belong to the group of one-step methods. This means that to calculate
at the point
it is enough to know the meaning at the node . It is natural to expect that if more information about a decision is used, several previous values ​​of the decision will be taken into account
,
etc., then the new value
it will be possible to find more accurately. This strategy is used in multi-step methods. To describe them, we introduce the notation
.

Representatives of multi-step methods are the Adams–Bashforth methods:


Method k-th order gives a local order error
or global – order .

These methods belong to the group of extrapolation methods, i.e. the new meaning is clearly expressed through the previous ones. Another type is interpolation methods. In them, at each step, you have to solve a nonlinear equation for a new value . Let's take the Adams–Moulton methods as an example:


To use these methods, you need to know several values ​​at the beginning of the count
(their number depends on the order of the method). These values ​​must be obtained by other methods, for example the Runge–Kutta method with a small step (to increase accuracy). Interpolation methods in many cases turn out to be more stable and allow larger steps to be taken than extrapolation methods.

In order not to solve a nonlinear equation at each step in interpolation methods, Adams predictor-correction methods are used. The bottom line is that the extrapolation method is first applied at the step and the resulting value
is substituted into the right side of the interpolation method. For example, in the second order method

Numerical solution of differential equations

Many problems in science and technology come down to solving ordinary differential equations (ODEs). ODEs are those equations that contain one or more derivatives of the desired function. In general, the ODE can be written as follows:

Where x is an independent variable, is the i-th derivative of the desired function. n is the order of the equation. The general solution of an nth order ODE contains n arbitrary constants, i.e. the general solution has the form .

To select a single solution, it is necessary to set n additional conditions. Depending on the method of specifying additional conditions, there are two different types of problems: the Cauchy problem and the boundary value problem. If additional conditions are specified at one point, then such a problem is called the Cauchy problem. Additional conditions in the Cauchy problem are called initial conditions. If additional conditions are specified at more than one point, i.e. for different values ​​of the independent variable, then such a problem is called a boundary value problem. The additional conditions themselves are called boundary or boundary conditions.

It is clear that when n=1 we can only talk about the Cauchy problem.

Examples of setting up the Cauchy problem:

Examples of boundary value problems:

It is possible to solve such problems analytically only for some special types of equations.

Numerical methods for solving the Cauchy problem for first-order ODEs

Formulation of the problem. Find a solution to the first order ODE

On the segment provided

When finding an approximate solution, we will assume that the calculations are carried out with a calculated step, the calculation nodes are the interval points [ x 0 , x n ].

The goal is to build a table

x i

x n

y i

y n

those. Approximate values ​​of y are sought at grid nodes.

Integrating the equation on the interval, we obtain

A completely natural (but not the only) way to obtain a numerical solution is to replace the integral in it with some quadrature formula of numerical integration. If we use the simplest formula for left rectangles of the first order

,

then we get explicit Euler formula:

Payment procedure:

Knowing, we find, then etc.

Geometric interpretation of Euler's method:

Taking advantage of what is at the point x 0 the solution is known y(x 0)= y 0 and the value of its derivative, we can write the equation of the tangent to the graph of the desired function at the point:. With a small enough step h the ordinate of this tangent, obtained by substituting into the right side of the value, should differ little from the ordinate y(x 1) solutions y(x) Cauchy problems. Therefore, the point of intersection of the tangent with the line x = x 1 can be approximately taken as the new starting point. Through this point we again draw a straight line, which approximately reflects the behavior of the tangent to at the point. Substituting here (i.e. the intersection with the line x = x 2), we obtain an approximate value y(x) at point x 2: etc. As a result for i-th point we obtain Euler’s formula.

The explicit Euler method has first order accuracy or approximation.

If you use the right rectangle formula: , then we come to the method

This method is called implicit Euler method, since calculating an unknown value from a known value requires solving an equation that is generally nonlinear.

The implicit Euler method has first order accuracy or approximation.

In this method, the calculation consists of two stages:

This scheme is also called the predictor-corrector method (predictive-correcting). In the first stage, the approximate value is predicted with low accuracy (h), and in the second stage this prediction is corrected so that the resulting value has second order accuracy.

Runge–Kutta methods: the idea of ​​constructing explicit Runge–Kutta methods p-th order is to obtain approximations to the values y(x i+1) according to a formula of the form

…………………………………………….

Here a n ,b nj , p n, – some fixed numbers (parameters).

When constructing the Runge–Kutta methods, the parameters of the function ( a n ,b nj , p n) are selected in such a way as to obtain the desired order of approximation.

Runge–Kutta scheme of fourth order of accuracy:

Example. Solve the Cauchy problem:

Consider three methods: explicit Euler method, modified Euler method, Runge–Kutta method.

Exact solution:

Calculation formulas using the explicit Euler method for this example:

Calculation formulas of the modified Euler method:

Calculation formulas for the Runge–Kutta method:

y1 – Euler’s method, y2 – modified Euler’s method, y3 – Runge Kutta’s method.

It can be seen that the most accurate is the Runge–Kutta method.

Numerical methods for solving systems of first-order ODEs

The methods considered can also be used to solve systems of first-order differential equations.

Let us show this for the case of a system of two first-order equations:

Explicit Euler method:

Modified Euler method:

Runge–Kutta scheme of fourth order of accuracy:

Cauchy problems for higher order equations are also reduced to solving systems of ODE equations. For example, consider Cauchy problem for a second order equation

Let's introduce a second unknown function. Then the Cauchy problem is replaced by the following:

Those. in terms of the previous problem: .

Example. Find a solution to the Cauchy problem:

On the segment.

Exact solution:

Really:

Let's solve the problem using the explicit Euler method, modified by the Euler and Runge-Kutta method with a step h=0.2.

Let's introduce the function.

Then we obtain the following Cauchy problem for a system of two first-order ODEs:

Explicit Euler method:

Modified Euler method:

Runge–Kutta method:

Euler circuit:

Modified Euler method:

Runge - Kutta scheme:

Max(y-y theory)=4*10 -5

Finite difference method for solving boundary value problems for ODE

Formulation of the problem: find a solution to a linear differential equation

satisfying the boundary conditions:. (2)

Theorem. Let . Then there is a unique solution to the problem.

This problem reduces, for example, to the problem of determining the deflections of a beam that is hinged at its ends.

Main stages of the finite difference method:

1) the area of ​​continuous change of argument () is replaced by a discrete set of points called nodes: .

2) The desired function of the continuous argument x is approximately replaced by the function of the discrete argument on a given grid, i.e. . The function is called a grid function.

3) The original differential equation is replaced by a difference equation with respect to the grid function. This replacement is called difference approximation.

Thus, solving a differential equation comes down to finding the values ​​of the grid function at grid nodes, which are found from solving algebraic equations.

Approximation of derivatives.

To approximate (replace) the first derivative, you can use the formulas:

- right difference derivative,

- left difference derivative,

Central difference derivative.

that is, there are many possible ways to approximate the derivative.

All these definitions follow from the concept of derivative as a limit: .

Based on the difference approximation of the first derivative, we can construct a difference approximation of the second derivative:

Similarly, we can obtain approximations of higher order derivatives.

Definition. The approximation error of the nth derivative is the difference: .

To determine the order of approximation, Taylor series expansion is used.

Let us consider the right-hand difference approximation of the first derivative:

Those. the right difference derivative has first by h order of approximation.

The same is true for the left difference derivative.

The central difference derivative has second order approximation.

The approximation of the second derivative according to formula (3) also has a second order of approximation.

In order to approximate a differential equation, it is necessary to replace all its derivatives with their approximations. Let's consider problem (1), (2) and replace the derivatives in (1):

As a result we get:

(4)

The order of approximation of the original problem is 2, because the second and first derivatives are replaced with order 2, and the rest - exactly.

So, instead of differential equations (1), (2), a system of linear equations is obtained for determination at grid nodes.

The diagram can be represented as:

i.e., we got a system of linear equations with a matrix:

This matrix is ​​tridiagonal, i.e. all elements that are not located on the main diagonal and the two diagonals adjacent to it are equal to zero.

By solving the resulting system of equations, we obtain a solution to the original problem.

Introduction

When solving scientific and engineering problems, it is often necessary to mathematically describe some dynamic system. This is best done in the form of differential equations ( DU) or systems of differential equations. Most often, this problem arises when solving problems related to modeling the kinetics of chemical reactions and various transfer phenomena (heat, mass, momentum) - heat transfer, mixing, drying, adsorption, when describing the movement of macro- and microparticles.

In some cases, a differential equation can be transformed into a form in which the highest derivative is expressed explicitly. This form of writing is called an equation resolved with respect to the highest derivative (in this case, the highest derivative is absent on the right side of the equation):

A solution to an ordinary differential equation is a function y(x) that, for any x, satisfies this equation in a certain finite or infinite interval. The process of solving a differential equation is called integrating a differential equation.

Historically, the first and simplest way to numerically solve the Cauchy problem for a first-order ODE is the Euler method. It is based on the approximation of the derivative by the ratio of finite increments of the dependent (y) and independent (x) variables between the nodes of a uniform grid:

where y i+1 is the desired value of the function at point x i+1.

The accuracy of Euler's method can be improved if a more accurate integration formula is used to approximate the integral - trapezoidal formula.

This formula turns out to be implicit with respect to y i+1 (this value is on both the left and right sides of the expression), that is, it is an equation with respect to y i+1, which can be solved, for example, numerically, using some iterative method (in such form, it can be considered as an iterative formula of the simple iteration method).

Composition of the course work: The course work consists of three parts. The first part contains a brief description of the methods. In the second part, the formulation and solution of the problem. In the third part - software implementation in computer language

The purpose of the course work: to study two methods for solving differential equations - the Euler-Cauchy method and the improved Euler method.

1. Theoretical part

Numerical differentiation

A differential equation is an equation containing one or more derivatives. Depending on the number of independent variables, differential equations are divided into two categories.

    Ordinary differential equations (ODE)

    Partial differential equations.

Ordinary differential equations are those equations that contain one or more derivatives of the desired function. They can be written as

independent variable

The highest order included in equation (1) is called the order of the differential equation.

The simplest (linear) ODE is equation (1) of order resolved with respect to the derivative

A solution to the differential equation (1) is any function that, after its substitution into the equation, turns it into an identity.

The main problem associated with linear ODE is known as the Kasha problem:

Find a solution to equation (2) in the form of a function satisfying the initial condition (3)

Geometrically, this means that it is required to find the integral curve passing through the point ) when equality (2) is satisfied.

Numerical from the point of view of the Kasha problem means: it is required to construct a table of function values ​​satisfying equation (2) and the initial condition (3) on a segment with a certain step. It is usually assumed that that is, the initial condition is specified at the left end of the segment.

The simplest numerical method for solving a differential equation is the Euler method. It is based on the idea of ​​graphically constructing a solution to a differential equation, but this method also provides a way to find the desired function in numerical form or in a table.

Let equation (2) with the initial condition be given, that is, the Kasha problem has been posed. Let's solve the following problem first. Find in the simplest way the approximate value of the solution at a certain point where is a fairly small step. Equation (2) together with the initial condition (3) specify the direction of the tangent of the desired integral curve at the point with coordinates

The tangent equation has the form

Moving along this tangent, we obtain an approximate value of the solution at point:

Having an approximate solution at a point, you can repeat the previously described procedure: construct a straight line passing through this point with an angular coefficient, and from it find the approximate value of the solution at the point

. Note that this line is not tangent to the real integral curve, since the point is not available to us, but if it is small enough, the resulting approximate values ​​will be close to the exact values ​​of the solution.

Continuing this idea, let's build a system of equally spaced points

Obtaining a table of values ​​of the required function

Euler's method consists of cyclically applying the formula

Figure 1. Graphical interpretation of Euler's method

Methods for numerical integration of differential equations, in which solutions are obtained from one node to another, are called step-by-step. Euler's method is the simplest representative of step-by-step methods. A feature of any step-by-step method is that starting from the second step, the initial value in formula (5) is itself approximate, that is, the error at each subsequent step systematically increases. The most used method for assessing the accuracy of step-by-step methods for approximate numerical solution of ODEs is the method of passing a given segment twice with a step and with a step

1.1 Improved Euler method

The main idea of ​​this method: the next value calculated by formula (5) will be more accurate if the value of the derivative, that is, the angular coefficient of the straight line replacing the integral curve on the segment, is calculated not along the left edge (that is, at point), but at the center of the segment. But since the value of the derivative between points is not calculated, we move on to the double sections with the center, in which the point is, and the equation of the straight line takes the form:

And formula (5) takes the form

Formula (7) is applied only for , therefore, values ​​​​cannot be obtained from it, therefore they are found using Euler’s method, and to obtain a more accurate result they do this: from the beginning, using formula (5) they find the value

(8)

At point and then found according to formula (7) with steps

(9)

Once found further calculations at produced by formula (7)

Lab 1

Numerical solution methods

ordinary differential equations (4 hours)

When solving many physical and geometric problems, one has to search for an unknown function based on a given relationship between the unknown function, its derivatives and independent variables. This ratio is called differential equation , and finding a function that satisfies the differential equation is called solving a differential equation.

Ordinary differential equation called equality

, (1)

in which

is an independent variable that changes in a certain segment, and - unknown function y ( x ) and her first n derivatives. called order of the equation .

The task is to find a function y that satisfies equality (1). Moreover, without stipulating this separately, we will assume that the desired solution has one or another degree of smoothness necessary for the construction and “legal” application of one or another method.

There are two types of ordinary differential equations

Equations without initial conditions

Equations with initial conditions.

Equations without initial conditions are equations of the form (1).

Equation with initial conditions is an equation of the form (1), in which it is required to find such a function

, which for some satisfies the following conditions: ,

those. at the point

the function and its first derivatives take on predetermined values.

Cauchy problems

When studying methods for solving differential equations using approximate methods main task counts Cauchy problem.

Let's consider the most popular method for solving the Cauchy problem - the Runge-Kutta method. This method allows you to construct formulas for calculating an approximate solution of almost any order of accuracy.

Let us derive the formulas of the Runge-Kutta method of second order accuracy. To do this, we represent the solution as a piece of a Taylor series, discarding terms with an order higher than the second. Then the approximate value of the desired function at the point x 1 can be written as:

(2)

Second derivative y "( x 0 ) can be expressed through the derivative of the function f ( x , y ) , however, in the Runge-Kutta method, instead of the derivative, the difference is used

selecting parameter values ​​accordingly

Then (2) can be rewritten as:

y 1 = y 0 + h [ β f ( x 0 , y 0 ) + α f ( x 0 + γh , y 0 + δh )], (3)

Where α , β , γ And δ – some parameters.

Considering the right-hand side of (3) as a function of the argument h , let's break it down into degrees h :

y 1 = y 0 +( α + β ) h f ( x 0 , y 0 ) + αh 2 [ γ f x ( x 0 , y 0 ) + δ f y ( x 0 , y 0 )],

and select the parameters α , β , γ And δ so that this expansion is close to (2). It follows that

α + β =1, αγ =0,5, α δ =0,5 f ( x 0 , y 0 ).

Using these equations we express β , γ And δ via parameters α , we get

y 1 = y 0 + h [(1 - α ) f ( x 0 , y 0 ) + α f ( x 0 +, y 0 + f ( x 0 , y 0 )], (4)

0 < α ≤ 1.

Now, if instead of ( x 0 , y 0 ) in (4) substitute ( x 1 , y 1 ), we get a formula for calculating y 2 approximate value of the desired function at the point x 2 .

In the general case, the Runge-Kutta method is applied to an arbitrary partition of the segment [ x 0 , X ] on n parts, i.e. with variable pitch

x 0 , x 1 , …, x n ; h i = x i+1 – x i , x n = X. (5)

Options α are chosen equal to 1 or 0.5. Let us finally write down the calculation formulas of the second order Runge-Kutta method with variable steps for α =1:

y i+1 =y i +h i f(x i + , y i + f(x i , y i)), (6.1)

i = 0, 1,…, n -1.

And α =0,5:

y i+1 =y i + , (6.2)

i = 0, 1,…, n -1.

The most used formulas of the Runge-Kutta method are formulas of the fourth order of accuracy:

y i+1 =y i + (k 1 + 2k 2 + 2k 3 + k 4),

k 1 =f(x i , y i), k 2 = f(x i + , y i + k 1), (7)

k 3 = f(x i + , y i + k 2), k 4 = f(x i +h, y i +hk 3).

For the Runge-Kutta method, Runge's rule is applicable to estimate the error. Let y ( x ; h ) – approximate value of the solution at the point x , obtained by formulas (6.1), (6.2) or (7) with step h , A p the order of accuracy of the corresponding formula. Then the error R ( h ) values y ( x ; h ) can be estimated using an approximate value y ( x ; 2 h ) solutions at a point x , obtained in increments 2 h :

(8)

Where p =2 for formulas (6.1) and (6.2) and p =4 for (7).