EntranceNumerical methods for solving differential equations. Modified Euler method
Differential equations are equations in which an unknown function appears under the derivative sign. The main task of the theory differential equations-- study of functions that are solutions to such equations.
Differential equations can be divided into ordinary differential equations, in which the unknown functions are functions of one variable, and partial differential equations, in which the unknown functions are functions of two and more variables.
The theory of partial differential equations is more complex and is covered in more complete or specialized mathematics courses.
Let's start studying differential equations with the simplest equation - a first-order equation.
Equation of the form
F(x,y,y") = 0,(1)
where x is an independent variable; y - the required function; y" - its derivative, is called a first-order differential equation.
If equation (1) can be resolved with respect to y", then it takes the form
and is called a first-order equation resolved with respect to the derivative.
In some cases, it is convenient to write equation (2) in the form f (x, y) dx - dy = 0, which is a special case of the more general equation
P(x,y)dx+Q(x,y)dy=O,(3)
where P(x,y) and Q(x,y) are known functions. The equation in symmetric form (3) is convenient because the variables x and y in it are equal, that is, each of them can be considered as a function of the other.
Let us give two basic definitions of the general and particular solutions of the equation.
A general solution to equation (2) in a certain region G of the Oxy plane is a function y = μ(x,C), depending on x and an arbitrary constant C, if it is a solution to equation (2) for any value of the constant C, and if for any initial conditions y x=x0 =y 0 such that (x 0 ;y 0)=G, there is a unique value of the constant C = C 0 such that the function y=q(x,C 0) satisfies the given initial conditions y=q(x 0 ,C).
A particular solution to equation (2) in the domain G is the function y=ts(x,C 0), which is obtained from the general solution y=ts(x,C) at a certain value of the constant C=C 0.
Geometrically, the general solution y = μ (x, C) is a family of integral curves on the Oxy plane, depending on one arbitrary constant C, and the particular solution y = μ (x, C 0) is one integral curve of this family passing through a given point (x 0; y 0).
Approximate solution of first order differential equations by Euler's method. The essence of this method is that the desired integral curve, which is a graph of a particular solution, is approximately replaced by a broken line. Let the differential equation be given
and initial conditions y |x=x0 =y 0 .
Let us find an approximate solution to the equation on the interval [x 0 ,b] that satisfies the given initial conditions.
Let's divide the segment [x 0 ,b] with points x 0<х 1 ,<х 2 <...<х n =b на n равных частей. Пусть х 1 --х 0 =х 2 -- x 1 = ... =x n -- x n-1 = ?x. Обозначим через y i приближенные значения искомого решения в точках х i (i=1, 2, ..., n). Проведем через точки разбиения х i - прямые, параллельные оси Оу, и последовательно проделаем следующие однотипные операции.
Let's substitute the values x 0 and y 0 into the right side of the equation y"=f(x,y) and calculate the slope y"=f(x 0,y 0) of the tangent to the integral curve at the point (x 0;y 0). To find the approximate value y 1 of the desired solution, we replace the integral curve on the segment [x 0 , x 1 ,] with a segment of its tangent at the point (x 0 ; y 0). In this case we get
y 1 - y 0 =f(x 0 ;y 0)(x 1 - x 0),
from where, since x 0, x 1, y 0 are known, we find
y1 = y0+f(x0;y0)(x1 - x0).
Substituting the values x 1 and y 1 into the right side of the equation y"=f(x,y), we calculate the slope y"=f(x 1,y 1) of the tangent to the integral curve at the point (x 1;y 1). Next, replacing the integral curve on the segment with a tangent segment, we find the approximate value of the solution y 2 at point x 2:
y 2 = y 1 +f(x 1 ;y 1)(x 2 - x 1)
In this equality, x 1, y 1, x 2 are known, and y 2 is expressed through them.
Similarly we find
y 3 = y 2 +f(x 2 ;y 2) ?x, …, y n = y n-1 +f(x n-1 ;y n-1) ?x
Thus, the desired integral curve in the form of a broken line was approximately constructed and approximate values y i of the desired solution at points x i were obtained. In this case, the values of i are calculated using the formula
y i = y i-1 +f(x i-1 ;y i-1) ?x (i=1,2, …, n).
The formula is the main calculation formula of the Euler method. Its accuracy is higher, the smaller the difference?x.
Euler's method refers to numerical methods that provide a solution in the form of a table of approximate values of the desired function y(x). It is relatively rough and is used mainly for approximate calculations. However, the ideas underlying Euler's method are the starting point for a number of other methods.
The degree of accuracy of Euler's method is, generally speaking, low. There are much more accurate methods for approximate solving differential equations.
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)
| x | Accurate solution | 0,1 | 0,05 | 0,025 | 0,00625 | 0,0015625 | 0,0007813 | 0,0001953 |
|---|---|---|---|---|---|---|---|---|
| 1 | 2 | 4 | 16 | 64 | 128 | 512 | ||
| 0 | 0,000000 | 0,000000 | 0,000000 | 0,000000 | 0,000000 | 0,000000 | 0,000000 | 0,000000 |
| 0,1 | 0,004837 | 0,000000 | 0,002500 | 0,003688 | 0,004554 | 0,004767 | 0,004802 | 0,004829 |
| 0,2 | 0,018731 | 0,010000 | 0,014506 | 0,016652 | 0,018217 | 0,018603 | 0,018667 | 0,018715 |
| 0,3 | 0,040818 | 0,029000 | 0,035092 | 0,037998 | 0,040121 | 0,040644 | 0,040731 | 0,040797 |
| 0,4 | 0,070320 | 0,056100 | 0,063420 | 0,066920 | 0,069479 | 0,070110 | 0,070215 | 0,070294 |
| 0,5 | 0,106531 | 0,090490 | 0,098737 | 0,102688 | 0,105580 | 0,106294 | 0,106412 | 0,106501 |
| 0,6 | 0,148812 | 0,131441 | 0,140360 | 0,144642 | 0,147779 | 0,148554 | 0,148683 | 0,148779 |
| 0,7 | 0,196585 | 0,178297 | 0,187675 | 0,192186 | 0,195496 | 0,196314 | 0,196449 | 0,196551 |
| 0,8 | 0,249329 | 0,230467 | 0,240127 | 0,244783 | 0,248202 | 0,249048 | 0,249188 | 0,249294 |
| 0,9 | 0,306570 | 0,287420 | 0,297214 | 0,301945 | 0,305423 | 0,306284 | 0,306427 | 0,306534 |
| 1 | 0,367879 | 0,348678 | 0,358486 | 0,363232 | 0,366727 | 0,367592 | 0,367736 | 0,367844 |
Relative errors of calculated function values for different h |
||||||||
| x | h | 0,1 | 0,05 | 0,025 | 0,00625 | 0,0015625 | 0,0007813 | 0,0001953 |
| N | 1 | 2 | 4 | 16 | 64 | 128 | 512 | |
| 0,1 | 100,00% | 48,32% | 23,76% | 5,87% | 1,46% | 0,73% | 0,18% | |
| 0,2 | 46,61% | 22,55% | 11,10% | 2,74% | 0,68% | 0,34% | 0,09% | |
| 0,3 | 28,95% | 14,03% | 6,91% | 1,71% | 0,43% | 0,21% | 0,05% | |
| 0,4 | 20,22% | 9,81% | 4,83% | 1,20% | 0,30% | 0,15% | 0,04% | |
| 0,5 | 15,06% | 7,32% | 3,61% | 0,89% | 0,22% | 0,11% | 0,03% | |
| 0,6 | 11,67% | 5,68% | 2,80% | 0,69% | 0,17% | 0,09% | 0,02% | |
| 0,7 | 9,30% | 4,53% | 2,24% | 0,55% | 0,14% | 0,07% | 0,02% | |
| 0,8 | 7,57% | 3,69% | 1,82% | 0,45% | 0,11% | 0,06% | 0,01% | |
| 0,9 | 6,25% | 3,05% | 1,51% | 0,37% | 0,09% | 0,05% | 0,01% | |
| 1 | 5,22% | 2,55% | 1,26% | 0,31% | 0,08% | 0,04% | 0,01% | |
Let us 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
| x | Accurate | Gün's method | Euler method | ||
|---|---|---|---|---|---|
| y | rel. error | y | rel. error | ||
| 0 | 0,000000 | 0,00000 | 0,00000 | ||
| 0,1 | 0,004837 | 0,00500 | 3,36% | 0,00000 | 100,00% |
| 0,2 | 0,018731 | 0,01903 | 1,57% | 0,01000 | 46,61% |
| 0,3 | 0,040818 | 0,04122 | 0,98% | 0,02900 | 28,95% |
| 0,4 | 0,070320 | 0,07080 | 0,69% | 0,05610 | 20,22% |
| 0,5 | 0,106531 | 0,10708 | 0,51% | 0,09049 | 15,06% |
| 0,6 | 0,148812 | 0,14940 | 0,40% | 0,13144 | 11,67% |
| 0,7 | 0,196585 | 0,19721 | 0,32% | 0,17830 | 9,30% |
| 0,8 | 0,249329 | 0,24998 | 0,26% | 0,23047 | 7,57% |
| 0,9 | 0,306570 | 0,30723 | 0,21% | 0,28742 | 6,25% |
| 1 | 0,367879 | 0,36854 | 0,18% | 0,34868 | 5,22% |
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
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:
| x | y | y" = f (x,y) | k=h f(x,y) | Δy |
| x0 | y0 | f(x0,y0) | k1(0) | k1(0) |
| x0 + ½ h | y0 + ½ k1(0) | f(x0 + ½ h, y0 + ½ k1(0)) | k2(0) | 2k2(0) |
| x0 + ½ h | y0 + ½ k2(0) | f(x0 + ½ h, y0 + ½ k2(0)) | k3(0) | 2k3(0) |
| x0 + h | y0 + k3(0) | f(x0 + h, y0 + k3(0)) | k4(0) | k4(0) |
| Δy0 = Σ / 6 | ||||
| x1 | y1 = y0 + Δy0 | f(x1,y1) | k1(1) | k1(1) |
| x1 + ½ h | y1 + ½ k1(1) | f(x1 + ½ h, y1 + ½ k1(1)) | k2(1) | 2k2(1) |
| x1 + ½ h | y1 + ½ k2(1) | f(x1 + ½ h, y1 + ½ k2(1)) | k3(1) | 2k3(1) |
| x1 + h | y1 + k3(1) | f(x1 + h, y1 + k3(1)) | k4(1) | k4(1) |
| Δy1 = Σ / 6 | ||||
| x2 | y2 = y1 + Δy1 | etc. | until you receive all the required | y values |
· 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.
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).
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)