Euler's method refers to numerical methods that give a solution in the form of a table of approximate values ​​of the desired function y (x)... It is relatively crude and is used primarily for rough calculations. However, the ideas underlying the Euler method are the starting point for a number of other methods.

Consider the first-order differential equation

with the initial condition

x= x 0 , y(x 0 )= y 0 (3.2)

It is required to find a solution to the equation on the segment [ a, b].

We split the segment [ a, b] into n equal parts and obtain the sequence NS 0 , NS 1 , NS 2 ,…, NS n, where x i = x 0 + ih (i=0,1,…, n), a h=(b- a)/ n- step of integration.

In Euler's method, the approximate values y (x i +1 )  y i +1 are calculated sequentially by the formulas:

y i + 1 = at i + hf (x i , y i ) (i = 0,1,2 ...) (3.3)

In this case, the sought-for integral curve y = y (x) passing through the point M 0 (NS 0 , at 0 ), replaced by a broken line M 0 M 1 M 2 … with peaks M i (x i , y i ) (i=0,1,2,…); every link M i M i +1 this broken line, called Euler's broken line, has a direction coinciding with the direction of the integral curve of equation (1), which passes through the point M i(see figure 2):

Figure 2. View of Euler's broken line

Modified Euler's method more accurate. First, the auxiliary values ​​of the required function are calculated at k + 1/2 in points NS k + 1/2, then the value of the right-hand side of equation (3.1) is found at the midpoint y k + 1/2 = f ( xk + 1/2 , y k + 1/2 ) and define at to + :

Then:
(3.4)

Formulas (3.4) are recurrent formulas of the Euler method.

To estimate the error at the point NS To carry out calculations at To step by step h, then with a step 2 h and take 1/3 of the difference between these values:

,

where y (x)- exact solution of the differential equation.

Euler's method is easily extended to systems of differential equations and to higher-order differential equations. The latter must be preliminarily reduced to a system of first-order differential equations.

3.2. Runge-Kutta method

Runge-Kutta methods have the following properties:

These methods are one-step: to find at k + 1 need information about the previous point (x To  y To ) 

The methods agree with the Taylor series up to the terms of the order h p where the degree R is different for different methods and is called a serial number or method order

They do not require the calculation of the derivatives of f (x y) but require the calculation of the function itself

Runge-Kutta Algorithm third order:

(3.5)

Runge-Kutta Algorithm fourth order:

(3.6)

Algorithms of the third and fourth orders require three and four calculations of the function at each step, respectively, but they are very accurate.

3.3. Adams method

The Adams method refers to multi-step schemes for solving differential equations, characterized by the fact that the solution at the current node depends on the data not in one previous or subsequent grid node, as is the case in one-step methods, but depends on the data in several neighboring nodes.

The idea behind the Adams methods is to use the values ​​calculated already in the previous steps to improve the accuracy

Y k -1 , Y k -2 , Y k -3 …

If values ​​are used in k previous nodes, then we speak of the k-step method of integrating the equation. One way to build multi-step methods is as follows. An interpolation polynomial of degree (k-1) -L k -1 (x) , which is used when integrating the differential equation by the expression:

In this case, the integral is expressed through the quadrature formula:

where λ l - quadrature coefficients.

The family of formulas obtained in this way is called explicitk -step Adams scheme... As you can see, at k=1 Euler's formula is obtained as a special case.

For example, for a formula of the 4th order we have:

(3.7)

y ( p ) k +1 - “forecast”, calculated using the values ​​at the previous points, f ( p ) k +1 - the approximate value of the function, calculated at the point of obtaining the forecast, y ( c ) k +1 - "correction" of the forecast value, y k +1 Is the required value according to Adams.

The advantage of this method for solving differential equations is that at each point only one value of the function is calculated F (x, y). The disadvantages include the impossibility of starting the multistep method from a single starting point, since for calculations based on k-step formula requires the value of the function value in k nodes. Therefore it is necessary (k-1) solution at first nodes x 1 , x 2 , ..., x k-1 to obtain with the help of any one-step method, for example, the Runge-Kutta method of the 4th order.

Another problem is the impossibility of changing the step in the solution process, which is easily implemented in one-step methods.

4. Brief description of the C ++ program and presentation of the results of its execution

Many problems of science and technology are reduced to solving ordinary differential equations (ODE). ODEs are equations that contain one or more derivatives of the desired function. In general

the desired function, n is the order of the equation. The general solution of an n-th order ODE contains n arbitrary constants

c 1, c 2, ..., c n, i.e. the general solution has the form y x, c 1, c 2, ..., c n. To select a single solution, n additional conditions must be specified. Depending on the method of assignment

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 for n 1 we can only talk about the Cauchy problem. Examples of setting the Cauchy problem:

Examples of boundary value problems:

some special types of equations, so the use of approximate solution methods is a necessity.

Approximate Methods for Solving the Cauchy Problem for First Order ODEs

It is required to find the solution y (x) of the first order ODE

In order to find all the values ​​of y i, you need to somehow

calculate the integral on the right-hand side of (5.4). Applying various quadrature formulas, we will obtain methods for solving problem (5.2), (5.3) of different orders of accuracy.

Euler's method

If to calculate the integral in (5.4) we use the simplest formula for left rectangles of the first order

The explicit Euler method has a first order of approximation. Method implementation. Since x 0, y 0, f x 0, y 0

are known, applying (5.5) successively, we define all y i: y 1 y 0 hf x 0, y 0, y 2 y 1 hf x 1, y 1,….

tangent, obtained by substituting the value x 1 x 0 h into the right-hand side, should differ little from the ordinate y x 1 of the solution

y x of the Cauchy problem. Therefore, the point x 1, y 1 of the intersection of the tangent with the straight line x x 1 can be approximately taken

Introduction

When solving scientific and engineering problems, it is often necessary to describe a dynamic system mathematically. This is best done in the form of differential equations ( DU) or a system of differential equations. Most often, such a problem arises when solving problems associated with 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, the differential equation can be transformed to a form in which the highest derivative is expressed in an explicit form. This form of notation is called an equation resolved with respect to the highest derivative (while the highest derivative is absent on the right side of the equation):

The solution of 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 the integration of the 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 the finite increments of the dependent (y) and independent (x) variables between the nodes of a uniform grid:

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

The accuracy of the Euler method can be increased by using a more accurate integration formula 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 for y i + 1, which can be solved, for example, numerically, using an iterative method (in such form it can be considered as an iterative formula of the simple iteration method).

The composition of the term paper: The term paper consists of three parts. In the first part, a brief description of the methods. In the second part, the statement and solution of the problem. In the third part - software implementation in the 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

An equation containing one or more derivatives is called a differential equation. Differential equations fall into two categories depending on the number of independent variables.

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

Any function that, after substituting it into the equation, turns it into an identity is called a solution to the differential equation (1).

The main problem associated with a linear ODE is known as the Kashi 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 Kashi problem means: it is required to construct a table of function values ​​that satisfies 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 set at the left end of the segment.

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

Let equation (2) be given with the initial condition, that is, the Kashi problem is posed. Let us first solve the following problem. Find in the simplest way the approximate value of the solution at some point where is a sufficiently small step. Equation (2) together with the initial condition (3) define the direction of the tangent to the sought-for integral curve at the point with coordinates

The tangent equation has the form

Moving along this tangent, we get the approximate value of the solution at the 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 straight line is not tangent to the real integral curve, since the point is not available to us, however, if small enough, then the resulting approximations will be close to the exact values ​​of the solution.

Continuing this idea, we construct a system of equally spaced points

Getting a table of values ​​of the required function

according to the Euler method consists in the cyclic application of the formula

Figure 1. Graphical interpretation of the Euler method

Methods for the 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 next step increases systematically. The most used method for assessing the accuracy of step-by-step methods for the approximate numerical solution of ODE is the method of double passage of a given segment 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 slope of the straight line replacing the integral curve on the segment will be calculated not along the left edge (that is, at a point), but along the center of the segment. But since the value of the derivative between the points is not calculated, we move on to the double sections of 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 the Euler method, while to obtain a more accurate result, they do this: from the beginning, using formula (5), find the value

(8)

At the point and then it is found by formula (7) with a step

(9)

After finding further calculations for is produced according to the formula (7)

Let us know the input dynamic sequence X(input signal) and model (method of converting an input signal into an output signal). The problem of determining the output signal y(t) (see Fig.10.1).

A dynamic system model can be represented by a differential equation. Basic equation of dynamics:

y" = f(x(t), y(t), t) .

The initial conditions are known at time zero t 0 : y(t 0) , x(t 0). To determine the output signal, note that by definition of the derivative:

We know the position of the system at point "1", we need to determine the position of the system at point "2". The points are separated from each other by a distance Δ t(fig.10.2). That is, the calculation of the system's behavior is done in steps. From point "1" we jump (discretely) to point "2", the distance between points along the axis t called calculation step Δ t .

The last formula is called Euler's formula.

Obviously to know the state of the system in the future y(t + Δ t) , it is necessary to the present state of the system y(t) add the change Δ y elapsed in time Δ t .

Let's take another look at this important relationship, deriving it from geometric considerations (Figure 10.3).

Let A be the point at which the state of the system is known. This is the "real" state of the system.

Draw a tangent to the trajectory of the system at point A. The tangent is the derivative of the function f(x(t), y(t), t) by variable t... The derivative at a point is always easy to calculate, it is enough to substitute the known variables (at the moment "Present" they are known) into the formula y" = f(x(t), y(t), t) .

Note that, by definition, the derivative is related to the angle of inclination of the tangent: y" = tg ( α ) so the angle α easy to calculate ( α = arctg ( y" ) ) and draw a tangent.

Draw a tangent to the intersection with the line t + Δ t... Moment t + Δ t corresponds to the "future" state of the system. Draw a line parallel to the axis t from point A to intersection with the line t + Δ t... The lines form a right-angled triangle ABC, one leg of which is equal to Δ t(famous). Also known is the angle α ... Then the second leg in a right-angled triangle ABC is: a = Δ t Tg ( α ) ... Now it is easy to calculate the ordinate of point B. It consists of two line segments - y(t) and a... The ordinate symbolizes the position of the system at the point y(t + Δ t) ... That is y(t + Δ t) = y(t) + a or further y(t + Δ t) = y(t) + Δ t Tg ( α ) or, substituting further, we have: y(t + Δ t) = y(t) + Δ t · y" and finally y(t + Δ t) = y(t) + Δ t · f(x(t), y(t), t) ... Again we got Euler's formula (from geometric considerations).

This formula can give accurate results only for very small Δ t(say when Δ t-> 0). At Δ t≠ 0 the formula gives the discrepancy between the true value y and calculated, equal ε , therefore, it must contain an approximate equal sign, or it must be written like this:

y(t + Δ t) = y(t) + Δ t · f(x(t), y(t), t) + ε .

Indeed. Take another look at fig. 10.3. Let's mentally move the line t + Δ t to the left (in fact, we will approximate the value of Δ t to zero). As it is easy to see, the distance BB * = ε , - that is, an error! - will shrink. In the limit (at Δ t-> 0) error value ε will be equal to zero.

So, replacing the real curve with a straight (tangent) line on the segment Δ t, we introduce an error into the solution, getting as a result not at the point "2" (see Fig. 10.2), but next to it, at the point "3". It is obvious that this numerical method at each step has a calculation error ε .

It can be seen from the figure that the smaller the value of Δ t, the smaller the calculation error will be ε ... That is, to calculate the behavior of the system for any length of time (for example, from t 0 to t k) to reduce the error at each step, steps Δ t make it as small as possible. To reach the point t k section (t k – t 0) is divided into segments of length Δ t; so everything will work out N = (t k – t 0)/Δ t steps. As a result of the calculation, you will have to apply the Euler formula for each step, that is N once. But it should be borne in mind that errors ε i on each i th step (in the simplest case) add up, and the total error quickly accumulates (see Fig. 10.4). And this is a significant disadvantage of this method. Although using this method it is possible to obtain (in numerical form) the solution of any differential equation (including analytically undecidable). By decreasing the step, we get more accurate solutions, but we should not forget that an increase in the number of steps leads to computational costs and a decrease in performance. In addition, with a large number of iterations, another significant error is introduced into the calculation due to the limited accuracy of computers and round-off errors.

Objective 1. A differential equation is given y" = 2ty ... The initial position of the system is set: y(0) = 1. It is required to find y(t), that is, the behavior of the system on the time interval t from 0 to 1.

Analytical method for solving problem 1

y" = 2ty .

Using the method of separating variables, we find:

y" /y = 2t

We will integrate from 0 to t i, then, according to the rules of integration, we have:

The obtained analytical solution is characterized by the fact that it is absolutely exact, but if the equation turns out to be in any way complicated, then the solution will not be found at all. The analytical solution path is not universal.

Numerical method for solving problem 1

The numerical solution assumes that the calculation will be carried out according to the Euler formula in a number of successive steps. At each step, the solution has its own error (see Fig. 10.2), since at each step the curve is replaced by a straight line.

In case of algorithmic implementation, the calculation is carried out in a cycle in which the t(counter t) and y :

The block diagram for the implementation of the method on a computer is shown in Fig. 10.5.

In the implementation of Stratum, the record will look like this (the presence of the "~" symbol when t ):

We will look for the meaning y of the previously considered example in numerical form in the interval from T= 0 to T= 1. Take the number of steps n= 10, then the increment step Δ t will be: Δ t= (1 - 0) / n= (1 - 0) / 10 = 0.1.

Note that the numerically calculated value ( y i+ 1) differs from the exact ( y exact. ), and the error (the difference between the columns y i+ 1 and y exact. ) in the process of calculation increases similarly to that shown in Fig. 10.4.

Now let's calculate the relative error σ for the calculated value y(1) obtained numerically in comparison with theoretical exact y theor. by the following formula:

σ = (1 - y calc. / y theory) 100%

and compare σ at different values ​​of Δ t .

If we change the value of the step Δ t, for example, to decrease the step, then the relative calculation error will also decrease. This is what happens when calculating the value y(1) with different step values ​​(see Table 10.2).

As you can see, with a decrease in the increment step Δ t the value of the relative error decreases, which means that the accuracy of the calculation increases.

Note that changing the step by a factor of 10 (from 1/10 to 1/100) leads to a change in the error value by about 10 times (from 14% to 2%). When the step is changed by 100 times, the error will also decrease by about 100 times. In other words, the step size and error for the Euler method are linearly related. If you want to reduce the error by 10 times - reduce the step by 10 times and increase the number of calculations by 10 times, respectively. This fact in mathematics is usually denoted by the symbol ε = O(Δ t) , and Euler's method is called the first-order precision method.

Since the error in the Euler method is large enough and accumulates from step to step, and the accuracy is proportional to the number of calculations, the Euler method is usually used for rough calculations, to assess the behavior of the system in principle. For accurate quantitative calculations, more accurate methods are used.

Notes (edit)

1. Each numerical method is accurate because the result differs from the theoretical one. The accuracy of the method depends on the step size. Different methods have different accuracy. The order of the dependence of the accuracy on the step size is denoted as O(h). Euler's method has the first order of accuracy; the dependence of the error on the step size is linear.
2. If, when decreasing the step, the limit y n strives for meaning y theor. , then the method is said to be convergent. Researchers are interested in the rate of convergence of the method.
3. The method must be stable. Stability is associated with some critical step size. With the manifestation of instability, there is a complete distortion of the qualitative picture of the calculation, "looseness" of the result.
4. When choosing a method, it is recommended to first achieve stability, and within the stability region, convergence of the result. Sustainability provides a quality picture. Convergence provides a quantitative result (see also Figure 10.10).

Set out in paragraphs. 1-4 will be illustrated by an example.

An example. Let be

Qualitatively, these equations describe the process of heat transfer between two bodies, the temperatures of which at some point in time will be denoted as A and B... Generally A and B- variables that change over time t... Finding the behavior of the system means finding how the temperatures will change. A(t) and B(t) .

It is intuitively clear that at the initial temperature difference A= 8 and B= 5 the temperatures of bodies should gradually equalize over time, since a hotter body will give energy to a colder one, and its temperature will decrease, and a colder body will receive energy from a hotter one, and its temperature will increase. The heat transfer process will end (that is, the changes will stop) when the temperatures of the two bodies become the same.

Let's do some calculations of the behavior A(t) and B(t) with different step sizes Δ t .

We will take different step sizes Δ t and find the corresponding values A and B in time by the following Euler formulas:

A new = A prev. + ( B prev. - A prev) Δ t ,
B new = B prev. + ( A prev. - B prev) Δ t .

Calculation for Δ t= 2 (Table 10.3).

The phenomenon of "looseness" is observed (see Fig. 10.6). Unstable solution. From physical considerations, it is obvious that two bodies cannot behave this way in the process of heat exchange.

Calculation for Δ t= 1 (Table 10.4).

The behavior of the solution of the system at the stability boundary is observed (see Fig. 10.7).

Calculation for Δ t= 0.5 (Table 10.5).

The solution is stable, corresponds to the correct qualitative picture (see Figure 10.8). The temperatures of bodies gradually approach, become the same over time. But the solution still has a large margin of error.

Calculation for Δ t= 0.1 (Table 10.6).

The solution is stable. The solution is more accurate (see Figure 10.9).

The role of changing the step size is illustrated in Fig. 10.10.