Department of Physical Chemistry SFU (RSU)
NUMERICAL METHODS AND PROGRAMMING
Materials for the lecture course
Lecturer – Art. Rev. Shcherbakov I.N.

SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Formulation of the problem

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.

Ordinary differential equation(ODE) of the nth order is the following equation, which contains one or more derivatives of the desired function y(x):

Here y(n) denotes the derivative of order n of some function y(x), x is the independent variable.

In some cases, a differential equation can be transformed into a form in which the highest derivative is expressed explicitly. This form of notation 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):

It is this form of recording that is accepted as standard when considering numerical methods for solving ODEs.

Linear differential equation is an equation that is linear with respect to the function y(x) and all its derivatives.

For example, below are linear ODEs of the first and second orders

Solving 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 by integrating the differential equation.

General solution of the ODE The nth order contains n arbitrary constants C 1 , C 2 , …, C n

This obviously follows from the fact that the indefinite integral is equal to the antiderivative of the integrand plus the constant of integration

Since n integrations are necessary to solve nth-order differential equations, n integration constants appear in the general solution.

Private solution The ODE is obtained from the general one if the constants of integration are given certain values ​​by defining some additional conditions, the number of which allows us to calculate all the uncertain constants of integration.

Exact (analytical) solution (general or particular) of a differential equation implies obtaining the desired solution (function y(x)) in the form of an expression from elementary functions. This is not always possible even for first-order equations.

Numerical solution DE (quotient) consists in calculating the function y(x) and its derivatives at some given points lying on a certain segment. That is, in fact, the solution to a nth-order differential equation of the form is obtained in the form of the following table of numbers (the column of values ​​of the highest derivative is calculated by substituting the values ​​into the equation):

For example, for a first-order differential equation, the solution table will have two columns - x and y.

The set of abscissa values ​​in which the value of a function is determined is called mesh, on which the function y(x) is defined. The coordinates themselves are called grid nodes. Most often, for convenience, they are used uniform grids, in which the difference between neighboring nodes is constant and is called grid spacing or integration step differential equation

Or , i= 1, …, N

For determining private solution it is necessary to set additional conditions that will allow the integration constants to be calculated. Moreover, there should be exactly n such conditions. For first order equations - one, for second - 2, etc. Depending on the way they are specified when solving differential equations, there are three types of problems:

· Cauchy problem (initial problem): Need to find something like this private solution differential equation that satisfies certain initial conditions specified at one point:

that is, a certain value of the independent variable (x 0), and the value of the function and all its derivatives up to order (n-1) at this point are given. This point (x 0) is called primary. For example, if a 1st order DE is being solved, then the initial conditions are expressed as a pair of numbers (x 0 , y 0)

This kind of problem occurs when solving ODE, which describe, for example, the kinetics of chemical reactions. In this case, the concentrations of substances at the initial moment of time are known ( t = 0) and it is necessary to find the concentrations of substances after a certain period of time ( t) . As an example, we can also cite the problem of heat transfer or mass transfer (diffusion), the equation of motion of a material point under the influence of forces, etc.

· Boundary value problem . In this case, the values ​​of the function and (or) its derivatives are known at more than one point, for example, at the initial and final moments of time, and it is necessary to find a particular solution to the differential equation between these points. The additional conditions themselves in this case are called regional (borderline) conditions. Naturally, the boundary value problem can be solved for ODEs of at least 2nd order. Below is an example of a second-order ODE with boundary conditions (function values ​​at two different points are given):

· Sturm-Liouville problem (eigenvalue problem). Problems of this type are similar to boundary value problems. When solving them, it is necessary to find at what values ​​of any parameter the solution DU satisfies boundary conditions (eigenvalues) and functions that are a solution to the DE for each parameter value (eigenfunctions). For example, many problems in quantum mechanics are eigenvalue problems.

Numerical methods for solving the Cauchy problem of first-order ODE

Let's consider some numerical methods for solving Cauchy problems(initial problem) ordinary differential equations of the first order. Let us write this equation in general form, resolved with respect to the derivative (the right side of the equation does not depend on the first derivative):

(6.2)

It is necessary to find the values ​​of the function y at given points of the grid if the initial values ​​are known, where there is the value of the function y(x) at the initial point x 0.

Let's transform the equation by multiplying by d x

And we integrate the left and right sides between the i-th and i+ 1st grid nodes.

(6.3)

We have obtained an expression for constructing a solution at the i+1 integration node through the values ​​of x and y at the i-th grid node. The difficulty, however, lies in the fact that the integral on the right side is an integral of an implicitly given function, which is generally impossible to find in analytical form. Numerical methods for solving ODEs in various ways approximate (approximate) the value of this integral to construct formulas for the numerical integration of ODEs.

Of the many methods developed for solving first-order ODEs, we consider methods , and . They are quite simple and give an initial idea of ​​approaches to solving this problem within the framework of a numerical solution.

Euler method

Historically, the first and simplest way to numerically solve the Cauchy problem for first-order ODEs 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 the uniform grid:

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

If we now transform this equation and take into account the uniformity of the integration grid, we obtain an iterative formula by which we can calculate y i+1, if y i is known at point x i:

Comparing Euler's formula with the general expression obtained earlier, it is clear that to approximately calculate the integral in, Euler's method uses the simplest integration formula - the formula of rectangles along the left edge of the segment.

The graphical interpretation of Euler's method is also easy (see figure below). Indeed, based on the form of the equation being solved (), it follows that the value is the value of the derivative of the function y(x) at the point x=x i - , and, thus, is equal to the tangent of the tangent angle drawn to the graph of the function y(x) at the point x =x i .

From the right triangle in the figure you can find

This is where Euler's formula comes from. Thus, the essence of Euler's method is to replace the function y(x) on the integration segment with a straight line tangent to the graph at point x=x i. If the desired function differs greatly from linear on the integration segment, then the calculation error will be significant. The error of the Euler method is directly proportional to the integration step:

Error~h

The calculation process is structured as follows. For given initial conditions x 0 And y 0 can be calculated

Thus, a table of function values ​​y(x) is constructed with a certain step ( h) By x on the segment. Error in defining value y(x i) in this case, the smaller the step length chosen, the smaller it will be h(which is determined by the accuracy of the integration formula).

For large h, Euler's method is very inaccurate. It gives an increasingly accurate approximation as the integration step decreases. If the segment is too large, then each section is divided into N integration segments and the Euler formula is applied to each of them with a step, that is, the integration step h is taken less than the step of the grid on which the solution is determined.

Example:

Using Euler's method, construct an approximate solution for the following Cauchy problem:

On a grid with a step of 0.1 in the interval (6.5)

Solution:

This equation has already been written in standard form, resolved with respect to the derivative of the desired function.

Therefore, for the equation being solved we have

Let us take the integration step equal to the grid step h = 0.1. In this case, only one value will be calculated for each grid node (N=1). For the first four grid nodes the calculations will be as follows:

The full results (accurate to the fifth decimal place) are given in the third column - h =0.1 (N =1). For comparison, the second column of the table shows the values ​​calculated from the analytical solution of this equation .

The second part of the table shows the relative error of the solutions obtained. It can be seen that at h =0.1 the error is very large, reaching 100% for the first node x =0.1.

Table 1 Solution of the equation by the Euler method (for columns, the integration step and the number of integration segments N between grid nodes are indicated)

Let's reduce the integration step by half, h = 0.05, in this case, for each grid node, the calculation will be carried out in two steps (N = 2). So, for the first node x =0,1 we get:

(6.6)

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 a form, it can be considered as an iterative formula of the simple iteration method). However, you can do it differently and approximately calculate the value of a function at a node i+1 using the usual formula:

,

which can then be used in the calculation according to (6.6).

This gives the method Guna or Euler's method with recalculation. For each integration node the following chain of calculations is performed

(6.7)

Thanks to a more accurate integration formula, the error of the Hün method is proportional to the square of the integration step.

Error~ h 2

The approach used in Gün's method is used to construct so-called methods forecast and correction, which will be discussed later.

Example:

Let's carry out calculations for equation () using Hün's method.

With integration step h =0.1 at the first grid node x 1 we obtain:

Which is much more accurate than the values ​​obtained by the Euler method with the same integration step. Table 2 below shows the comparative results of calculations for h = 0.1 of the Euler and Gün methods.

Table 2 Solution of the equation by Euler and Gün methods

Let us note a significant increase in the accuracy of calculations of the Hün method compared to the Euler method. Thus, for node x =0.1, the relative deviation of the function value determined by the Huyn method turns out to be 30 (!) times less. The same accuracy of calculations using Euler's formula is achieved when the number of integration segments N is approximately 30. Consequently, when using the Hün method with the same accuracy of calculations, it will take approximately 15 times less computer time than when using the Euler method.

Checking the stability of the solution

A solution to an ODE at some point x i is called stable if the value of the function found at this point y i changes little as the integration step decreases. To check stability, therefore, it is necessary to carry out two calculations of the value ( y i) – with integration step h and with a reduced (for example, two) step size

As a stability criterion, you can use the smallness of the relative change in the obtained solution when the integration step is reduced (ε is a predetermined small value)

This check can be carried out for all solutions over the entire range of values x. If the condition is not met, then the step is again divided in half and a new solution is found, etc. until a stable solution is obtained.

Runge-Kutta methods

Further improvement in the accuracy of solving a first-order ODE is possible by increasing the accuracy of the approximate calculation of the integral in the expression.

We have already seen the advantage of moving from integrating using the rectangle formula () to using the trapezoid formula () when approximating this integral.

Using the well-proven Simpson formula, you can obtain an even more accurate formula for solving the Cauchy problem for first-order ODE - the Runge-Kutta method widely used in computing practice.

The advantage of Adams' multi-step methods for solving ODEs is that at each node only one value of the right-hand side of the ODE is calculated - the function F(x,y). The disadvantages include the impossibility of starting a multi-step method from a single starting point, since calculations using the k-step formula require knowledge of the value of the function at k nodes. Therefore, it is necessary to obtain a (k-1) solution at the first nodes x 1, x 2, ..., x k-1 using some one-step method, for example the method

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 steps 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 steps 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

Main issues discussed in the lecture:

1. Statement of the problem

2. Euler's method

3. Runge-Kutta methods

4. Multi-step methods

5. Solution of the boundary value problem for a linear differential equation of 2nd order

6. Numerical solution of partial differential equations

1. Statement of the problem

The simplest ordinary differential equation (ODE) is a first-order equation solved with respect to the derivative: y " = f (x, y) (1). The main problem associated with this equation is known as the Cauchy problem: find a solution to equation (1) in the form of a function y (x), satisfying the initial condition: y (x0) = y0 (2).
DE of nth order y (n) = f (x, y, y",:, y(n-1)), for which the Cauchy problem is to find a solution y = y(x) that satisfies the initial conditions:
y (x0) = y0 , y" (x0) = y"0 , :, y(n-1)(x0) = y(n-1)0 , where y0 , y"0 , :, y(n- 1)0 - given numbers, can be reduced to a first order DE system.

· Euler method

The Euler method is based on the idea of ​​graphically constructing a solution to a differential equation, but the same method also provides a numerical form of the desired function. Let equation (1) with initial condition (2) be given.
Obtaining a table of values ​​of the desired function y (x) using the Euler method involves cyclically applying the formula: , i = 0, 1, :, n. To geometrically construct Euler’s broken line (see figure), we select the pole A(-1,0) and plot the segment PL=f(x0, y0) on the ordinate axis (point P is the origin of coordinates). Obviously, the angular coefficient of the ray AL will be equal to f(x0, y0), therefore, in order to obtain the first link of the Euler broken line, it is enough to draw straight line MM1 from point M parallel to the ray AL until it intersects with the straight line x = x1 at some point M1(x1, y1). Taking the point M1(x1, y1) as the initial one, we plot the segment PN = f (x1, y1) on the Oy axis and draw a straight line through the point M1 M1M2 | | AN until the intersection at point M2(x2, y2) with the line x = x2, etc.

Disadvantages of the method: low accuracy, systematic accumulation of errors.

· Runge-Kutta methods

The main idea of ​​the method: instead of using partial derivatives of the function f (x, y) in working formulas, use only this function itself, but at each step calculate its values ​​at several points. To do this, we will look for a solution to equation (1) in the form:


Changing α, β, r, q, we will obtain various versions of the Runge-Kutta methods.
For q=1 we obtain Euler's formula.
With q=2 and r1=r2=½ we obtain that α, β= 1 and, therefore, we have the formula: , which is called the improved Euler-Cauchy method.
For q=2 and r1=0, r2=1 we obtain that α, β = ½ and, therefore, we have the formula: - the second improved Euler-Cauchy method.
For q=3 and q=4, there are also entire families of Runge-Kutta formulas. In practice, they are used most often, because do not increase errors.
Let's consider a scheme for solving a differential equation using the Runge-Kutta method of 4th order of accuracy. Calculations when using this method are carried out according to the formulas:

It is convenient to include them in the following table:

· Multi-step methods

The methods discussed above are the so-called methods of step-by-step integration of a differential equation. They are characterized by the fact that the value of the solution at the next step is sought using the solution obtained only at one previous step. These are the so-called one-step methods.
The main idea of ​​multi-step methods is to use several previous solution values ​​when calculating the solution value at the next step. Also, these methods are called m-step methods based on the number m used to calculate previous solution values.
In the general case, to determine the approximate solution yi+1, m-step difference schemes are written as follows (m 1):
Let's consider specific formulas that implement the simplest explicit and implicit Adams methods.

Explicit 2nd order Adams method (2-step explicit Adams method)

We have a0 = 0, m = 2.
Thus, these are the calculation formulas of the explicit Adams method of the 2nd order.
For i = 1, we have an unknown y1, which we will find using the Runge-Kutta method for q = 2 or q = 4.
For i = 2, 3, : all necessary values ​​are known.

Implicit 1st order Adams method

We have: a0 0, m = 1.
Thus, these are the calculation formulas of the implicit Adams method of the 1st order.
The main problem with implicit schemes is the following: yi+1 is included in both the right and left sides of the presented equality, so we have an equation for finding the value of yi+1. This equation is nonlinear and is written in a form suitable for an iterative solution, so we will use the simple iteration method to solve it:
If step h is chosen well, then the iterative process converges quickly.
This method is also not self-starting. So to calculate y1 you need to know y1(0). It can be found using Euler's method.