Solution of ordinary differential equations. Types of differential equations, solution methods
MATLAB 2 tutorial
1. Solution of ordinary differential equations (ODEs)
1.1. ODE solvers in MATLAB
1.2. First order ODE solution
1.3. ODE solution n-th order
1.4. Solution of ODE systems
2. Dynamic systems (DS)
2.1. Types of DC
2.2. Phase space DS
2.3. Kinematic interpretation of the control system
2.4. DS evolution
2.5. Pendulum equations
2.6. An example of solving the Van der Pol ODE system
3. Qualitative analysis of linear DS
3.1. Singular points of linear DE
3.2. Phase portrait of linear DCs
4. Qualitative analysis of nonlinear DS
4.1. Logistic display
4.2. Singular points of nonlinear DE, bifurcation
4.3. Phase portrait of nonlinear DS
5. Van der Pol equation
6. Lorentz attractor
7. Enon mapping
Lab #1
PURPOSE OF THE WORK
To form students' ideas about the use of remote control in various fields; to instill the ability to solve the Cauchy problem for remote control y¢ = f(x, y) on the segment [ a, b] for a given initial condition at 0 = f(x 0).
Solving Ordinary Differential Equations (ODEs)
ODE solvers in MATLAB
The analysis of the behavior of many systems and devices in dynamics, as well as the solution of many problems in the theory of oscillations, is usually based on the solution of systems of ODEs. They are usually represented as a system of differential equations (DE) of the first order in the Cauchy form:
Equations (1.2) cannot usually be reduced analytically to form (1.1). However, the numerical solution of special difficulties does not cause enough to determine f(y, t) solve (1.2) numerically with respect to the derivative for given y and t.
ODE solvers
To solve ODE systems in MATLAB, various numerical methods are implemented. Their implementations are named solvers ODE.
In this section, the generic name solver means one of the possible numerical methods for solving ODEs: ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, bvp4c, or pdepe.
The solvers implement the following methods for solving DE systems:
Ode45 one-step explicit 4th and 5th order Runge-Kutta methods modified by Dormand and Prinz. This is the classic method recommended for an initial trial of the solution. In many cases, it gives good results if the system of equations to be solved is not rigid.
Ode23 one-step explicit Runge-Kutta methods of the 2nd and 4th orders in the modification of Bogacki and Champin. With a moderate rigidity of the ODE system and low accuracy requirements, this method can give a gain in solution speed.
Ode113 multi-step Adams-Bashworth-Multon method of variable order predictor-corrector class. It is an adaptive method that can provide a high solution accuracy.
Ode15s multi-step variable order method (from 1 to 5, default 5) using numerical "backward differentiation" formulas. This is an adaptive method and should be used if the ode45 solver does not provide a solution and the control system is stiff.
Ode23s is a one-step method using a modified 2nd order Rosenbrock formula. It can provide high computational speed with low accuracy of solving a rigid DE system.
Ode23t implicit trapezoid method with interpolation. This method gives good results when solving problems describing oscillatory systems with an almost harmonic output signal. For moderately hard systems, remote control can give a high accuracy of the solution.
Ode23tb Runge Kutta's implicit method at the beginning of the solution and a method using 2nd order "backward differentiation" formulas thereafter. Despite the relatively low accuracy, this method may be more efficient than ode15s.
Bvp4c is used for the problem of boundary values of systems of control systems of the form y′ = f(t, y), F(y(a), y(b), p) = 0 (full form system of Cauchy equations). The problems he solves are called two-point boundary value problems, since the solution is sought by setting the boundary conditions both at the beginning and at the end of the solution interval.
All solvers can solve systems of explicit equations y′ = F(t, y), and to solve stiff systems of equations, it is recommended to use only special solvers ode15s, ode23s, ode23t, ode23tb.
Using ODE Solvers
tspan vector defining the integration interval [ t 0 t final]. To obtain solutions at specific points in time t 0 , t 1 , …, t final(ordered in decreasing or increasing order) should be used tspan = [t 0 t 1 … t final];
y 0 vector of initial conditions;
Options argument created by the odeset function (another odeget or bvpget function (bvp4c only) allows you to display the options set by default or using the odeset/bvpset function);
p 1, p 2,… arbitrary parameters passed to the function F;
T, Y decision matrix Y, where each row corresponds to the time returned in the column vector T.
Let's move on to describing the syntax of functions for solving remote control systems (the name solver means any of the functions presented above).
[T,Y]=solver(@ F,tspan,y 0) integrates a remote control system of the form y′ = F(t, y) on the interval tspan with initial conditions y 0 . @F descriptor of an ODE function (you can also specify a function in the form " F"). Each row in the solution array Y corresponds to the time value returned in the column vector T.
[T,Y]=solver(@ F,tspan,y 0 ,options) gives a solution similar to the one above, but with options determined by the values of the options argument created by the odeset function. Commonly used parameters include the relative error tolerance RelTol (default 1e3) and the absolute error tolerance vector AbsTol (all components default 1e6).
[T,Y]=solver(@ F,tspan,y 0 ,options, p 1 ,p 2 ...) gives a solution similar to the one above, passing additional parameters p 1 , p 2 , … in m-file F whenever it is called. Use options= if no options are given.
First order ODE solution
WORK PROCEDURE
· title page;
initial data of the variant;
· the solution of the problem;
the results of solving the problem.
Example
Find a solution to the differential equation on the segment for which at(1,7) = 5,3.
Create a user function in the Command Window
[email protected](x,y);
In the function syntax @(x,y) x independent variable, y dependent variable, x-cos( y/pi) the right side of the DE.
The solution process is carried out by calling the solver in the Command Window with the following statement:
Ode23(g,,);
The construction of a graph with a grid is carried out by the following operators:
The result is shown in fig. 1.1
Rice. 1.2.1. Visualization numerical solution
EXERCISE
1. Find first-order DE solutions , satisfying the initial conditions y(x 0 ) = y 0 on the interval [ a,b].
2. Construct graphs of the function.
Task options.
option number | y(x 0 )=y 0 | [a,b] | |
y 0 (1,8)=2,6 | |||
y 0 (0,6)=0,8 | |||
y 0 (2,1)=2,5 | |||
y 0 (0,5)=0,6 | |||
y 0 (1,4)=2,2 | |||
y 0 (1,7)=5,3 | |||
y 0 (1,4)=2,5 | |||
y 0 (1,6)=4,6 | |||
y 0 (1,8)=2,6 | |||
y 0 (1,7)=5,3 | |||
y 0 (0,4)=0,8 | |||
y 0 (1,2)=1,4 |
Lab #2
Solution of ODE systems
PURPOSE OF THE WORK
To form students' ideas about the use of remote control systems in various fields; to instill the ability to solve the Cauchy problem for remote control systems.
WORK PROCEDURE
1. Explore theoretical part. Complete the tasks corresponding to the number of your option and show them to the teacher.
2. Prepare a lab report, which should contain:
· title page;
initial data of the variant;
· the solution of the problem;
the results of solving the problem.
Example
Solve the system
using the ode23() solver.
Solution:
1. Create in the editor an m-file of the function for calculating the right parts of the DE.
Let the name in the file editor be sisdu.m, then the function might look like this:
function z=sisdu(t,y)
z1=-3*y(2)+cos(t)-exp(t);
z2=4*y(2)-cos(t)+2*exp(t);
>> t0=0;tf=5;y0=[-3/17,4/17];
>> =ode23("sisdu",,y0);
>>plot(t,y)
Rice. 1.3.1. Visualization of the numerical solution obtained using the ode23 function.
1. What does it mean to solve the Cauchy problem for a control system?
2. What methods exist for solving remote control systems?
EXERCISE
1. Find a remote control solution
satisfying the initial conditions on the interval ;
2. Build graphs of functions.
For example, the decision function of the 8th option is given:
function z=ssisdu(t,y)
% option 8
z1=-a*y(1)+a*y(2);
z2=a*y(1)-(a-m)*y(2)+2*m*y(3);
z3=a*y(2)-(a-m)*y(3)+3*m*y(4);
z4=a*y(3)-3*m*y(4);
>> =ode23("ssisdu",,);
>>plot(t,100*y)
Rice. 1.3.2. Visualization of the numerical solution obtained using the ode23 function.
Task options.
option number | Tasks | |
a | m | |
0,1 | 1,2 | |
0,2 | 1,5 | |
0,3 | 1,7 | |
0,4 | 1,9 | |
0,5 | ||
0,6 | 1,9 | |
0,7 | 2,3 | |
0,8 | 2,7 | |
0,9 | ||
0,1 | 1,5 | |
0,2 | 1,1 | |
0,3 |
Lab #3
1.4Solution of ODE n-th order
PURPOSE OF THE WORK
To form students' ideas about the use of higher-order control systems in various fields; to instill the ability to solve the Cauchy problem for higher-order control systems with the help of application programs; develop the skills to check the results obtained.
WORK PROCEDURE
1. Study the theoretical part. Complete the tasks corresponding to the number of your option and show them to the teacher.
2. Prepare a lab report, which should contain:
· title page;
initial data of the variant;
· the solution of the problem;
the results of solving the problem.
Example 1
Solve second order DE under given initial conditions .
Solution:
First, we bring the DE to the system:
1. Create an m-file of the function for calculating the right parts of the DE.
Let the file name be sisdu_3.m, then the function might look like this:
function z=sisdu_3(x,y)
z2=6*x*exp(x)+2*y(2)+y(1);
2. Do the following:
>> x0=0;xf=10;y0=;
>> =ode23("sisdu_3",,y0);
>>plot(x,y(:,1))
Rice. 1.4.1. Visualization of the numerical solution obtained using the ode23 function.
EXAMPLE QUESTIONS ON WORK DEFENSE
1. What does it mean to solve the Cauchy problem for higher-order DEs?
2. How to bring DU m th order to the remote control system?
EXERCISE
1. Find a DE solution that satisfies the initial conditions on the interval .
2. Build graphs of functions.
Task options.
option number | Tasks | |
Equations | Initial conditions | |
|
||
|
||
|
||
|
||
| ||
|
Lab #4 – 5
Dynamic Systems (DS)
PURPOSE OF THE WORK
Acquaintance of students with the basic concepts of DS, their classification, phase space of DS, kinematic interpretation of the DE system, evolution of DS. The equation of motion of the pendulum. Dynamics of the Van der Pol oscillator.
2. Dynamic system (DS) a mathematical object corresponding to real systems (physical, chemical, biological, etc.), whose evolution is uniquely determined by the initial state. The DS is determined by a system of equations (differential, difference, integral, etc.) that allow the existence of a unique solution for each initial condition on an infinite time interval.
The state of the DS is described by a set of variables chosen for reasons of naturalness of their interpretation, simplicity of description, symmetry, etc. The set of DS states forms a phase space, each state corresponds to a point in it, and evolution is represented by (phase) trajectories. To determine the proximity of states, the concept of distance is introduced in the phase space of DS. The set of states at a fixed point in time is characterized by a phase volume.
The description of the DS in the sense of specifying the law of evolution also allows for a great variety: it is carried out with the help of differential equations, discrete mappings, with the help of graph theory, the theory of Markov chains, etc. The choice of one of the description methods sets a specific type of mathematical model of the corresponding DS.
Mathematical model of DS is considered given if dynamic variables (coordinates) of the system are introduced that uniquely determine its state, and the law of evolution of the state in time is specified.
Depending on the degree of approximation, different mathematical models can be assigned to the same system. The study of real systems goes along the path of studying the corresponding mathematical models, the improvement and development of which is determined by the analysis of experimental and theoretical results when they are compared. In this regard, by a dynamic system we will understand its mathematical model. Investigating the same DS (for example, the movement of a pendulum), depending on the degree of consideration of various factors, we will obtain different mathematical models.
Make full rotations
Suppose we want to find, with an accuracy of 0.001, smallest value the initial speed required to make a pendulum, starting from its original position, perform a full rotation once. It will be useful to show solutions that correspond to several different initial speeds on the same graph.
First, we will consider integer speed values in the range from 5 to 10.
>> for a = 5:10
Ode45(g, , );
plot(xa(:, 1), xa(:, 2))
>> hold off
Initial velocities 5, 6, and 7 are not large enough (Fig. 4) to increase the angle by more than π, but initial velocities 8, 9, and 10 are sufficient to make the pendulum complete a revolution. Let's see what happens between 7 and 8.
Second step
>> for a = 7.0:0.2:8.0
plot(xa(:, 1), xa(:, 2))
>> hold off
You can see (Figure 6) that the answer is somewhere between 7.2 and 7.4. Let's make one more clarification.
Third step
>> for a = 7.2:0.05:7.4
Ode45(g, , );
plot(xa(:, 1), xa(:, 2))
It should be concluded that the lowest required speed with an accuracy of 0.01 is somewhere between 7.25 and 7.3 (Fig. 7 and 8).
Fourth step
for a = 7.25:0.01:7.3
Ode45(g, , );
plot(xa(:, 1), xa(:, 2))
For a more accurate analysis, let's increase the area of the graph where the oscillation mode changes.
It can be seen that the lowest required speed is somewhere between 7.29 and 7.3.
It is necessary to continue finding a more accurate value of the oscillation mode change rate.
Dynamics of the Van der Pol oscillator at w 2 = 2 and c = 1
limit cycle- stable mode of periodic oscillations in nonlinear systems after the completion of transients
The function file for these parameters has the following form (see Andreev 2013, p. 117):
function dydt = vdp1(t,y)
dydt = zeros(2,1); dydt(1) = y(2);
dydt(2) = 1*(1-y(1).^2).*y(2)-2*y(1);
With the initial position on the limit cycle (x0 = 2, v0 = 0), the file function call has the following form:
Ode23(@vdp1,,);
The following statement makes it possible to obtain the dependence of x and v on time
plot(t,y(:,1),t,y(:,2)), grid on
The result is shown in fig. one.
The following operator makes it possible to obtain a phase portrait of the system (Fig. 2):
plot(y(:,1),y(:,2)), grid on
If we take the initial data outside the loop (x 0 = -0.5, v 0 = 5), the file function call looks like this:
Ode23(@vdp1,,[-0.5;5]);
The result is shown in fig. 3.
The phase portrait is shown in Fig. four.
If we take the initial data inside the loop (x 0 \u003d -0.05, v 0 \u003d -0.05), the call to the file function looks like this:
Ode23(@vdp1,,[-0.05; -0.05]);
The result is shown in fig. 5.
The phase portrait is shown in Fig. 6.
Lab #6 – 7
PURPOSE OF THE WORK
Consider a linear homogeneous system with constant coefficients:
System (1) has a unique zero equilibrium position if the system matrix determinant:
notice, that a + d= tr A(matrix trace) and adbc= det A.
The classification of rest points in the case when det A ≠ 0 is given in the table:
Stability of rest points
The eigenvalues of the matrix of system (1) uniquely determine the nature of the stability of the equilibrium positions:
Phase portraits
An infinite number of resting points
If det A= 0, then system (1) has an infinite set of equilibrium positions. In this case, three cases are possible:
In the second case, any rest point is Lyapunov stable. In the first case, only if l 2< 0.
The direction on the phase curve indicates the direction in which the phase point moves along the curve as t increases.
Rules for determining the type of rest point
It is possible to determine the type of the rest point and the nature of its stability without finding the eigenvalues of the matrix of system (1), but knowing only its trace tr A and determinant det A.
Bifurcation diagram
Algorithm for constructing the LDS phase portrait (1)
1. Determine the equilibrium positions by solving the system of equations:
2. Find eigenvalues system matrices by solving the characteristic equation:
3. Determine the type of rest point and draw a conclusion about stability.
4. Find the equations of the main horizontal and vertical isoclines and plot them on the phase plane.
5. using a saddle or node, find those phase trajectories that lie on straight lines passing through the origin.
6. Draw phase trajectories.
7. Determine the direction of movement along the phase trajectories, indicating it with arrows on the phase portrait.
Main isoclines
Note that the rest point on the phase plane is the intersection of the principal isoclines. The vertical isocline on the phase plane will be marked with vertical strokes, and the horizontal with horizontal ones.
Phase trajectories
If the equilibrium position is saddle or node, then there are phase trajectories that lie on straight lines passing through the origin.
The equations of such lines can be sought in the form y = kx. Substituting y = kx into the equation:
for determining k we get:
(4) |
(The equations of straight lines containing phase trajectories can also be sought in the form x = ky . Then, to find the coefficients, one should solve the equation
We give a description of the phase trajectories depending on the number and multiplicity of the roots of equation (4).
Phase trajectories
* If the equations of lines are sought in the form x = ky, then it will be straight x= k 1 y and y = 0.
If the equilibrium position is center, then the phase trajectories are ellipses.
If the equilibrium position is focus, then the phase trajectories are spirals.
In the case when LDS has line of rest points, then you can find the equations of all phase trajectories by solving the equation:
Its first integral ax + by = C defines a family of phase lines.
Direction of travel
If the equilibrium position is node or focus, then the direction of motion along the phase trajectories is uniquely determined by its stability (to the origin) or instability (from the origin).
True, in the case focus you also need to set the direction of twisting (untwisting), spiraling clockwise or counterclockwise. This can be done, for example, like this. Determine the sign of the derivative y’(t) at the points of the axis x.
If the equilibrium position is center, then the direction of movement along the phase trajectories (clockwise or counterclockwise) can be determined in the same way as the direction of "twisting (untwisting)" of the trajectory in the case of focus is set.
Therefore, if the equilibrium position saddle, then it is enough to set the direction of motion along some trajectory. And then you can unambiguously establish the direction of movement along all other trajectories.
Direction of travel (saddle)
To set the direction of movement along phase trajectories in the case saddles, you can use one of the following methods:
1 way
Determine which of the two separatrices corresponds to a negative eigenvalue. Movement along it occurs to a point of rest.
2 way
Determine how the abscissa of the moving point changes along any of the separatrices. For example, for y = k 1 x we have:
If a x(t) → 0 as t→ +∞, then the motion along the separatrix y = k 1 x going to rest point.
If a x(t) → ±∞ as t→+∞, then the movement occurs from the rest point.
3 way
If the x-axis is not a separatrix, determine how the ordinate of the moving point changes along the phase trajectory when crossing the axis x.
When if x> 0, then the ordinate of the point increases and, therefore, the motion along the phase trajectories intersecting the positive part of the axis x, occurs from the bottom up. If the ordinate decreases, then the movement will occur from top to bottom.
If we determine the direction of motion along the phase trajectory that intersects the axis y, it is better to analyze the change in the abscissa of the moving point.
4 way
Example 1
1. The system has a single zero equilibrium position, since
det A = 6 ≠ 0.
2. Having constructed the corresponding characteristic equation l 2 6 = 0, we find its roots l 1,2 = ± Ö6. Real roots and different sign. Therefore, the equilibrium position saddle.
3. The separatrices of the saddle are sought in the form y = kx.
4. Vertical isocline: x + y = 0.
Horizontal isocline: x 2y = 0.
Example 2
A = 10 ≠ 0.
2. Having constructed the corresponding characteristic equation l 2 7l + 10 = 0,
find its roots l 1 = 2, l 2 = 5. Therefore, the equilibrium position unstable knot.
3. Direct: y = kx.
4. Vertical isocline: 2 x + y = 0.
Horizontal isocline: x + 3y = 0.
Example 3
1. The system has a single zero equilibrium position, since det A = 18 ≠ 0.
2. Having constructed the corresponding characteristic equation l 2 + 3l + 18 = 0,
find its discriminant D = 63. Since D< 0, то корни уравнения комплексные, причем Re l 1,2 = 3/2. Следовательно, положение равновесия steady focus.
3. Vertical isocline: x + 4y = 0.
Horizontal isocline: 2x y = 0.
Phase trajectories are spirals, the movement along which occurs to the origin. The directions of "twisting trajectories" can be determined as follows.
Example 4
1. The system has a single zero equilibrium position, since det A = 3 ≠ 0.
2. Having constructed the corresponding characteristic equation l 2 +3 = 0, we find its roots l 1,2 = ± iÖ3. Therefore, the equilibrium position center.
3. Vertical isocline: x 4y = 0.
Horizontal isocline: x y = 0.
Phase trajectories of the ellipse system.
The direction of movement along them can be set, for example, like this.
Example 5 (degenerate node)
1. The system has a single zero equilibrium position, since det A= 4 ¹ 0.
2. Having constructed the corresponding characteristic equation l 2 + 4l + 4 = 0, we find its roots l 1 = l 2 = 2. Therefore, the equilibrium position unstable degenerate knot.
4. Vertical isocline: 2x + y = 0.
Horizontal isocline: x + 3y = 0.
Example 6
Since the determinant of the system matrix det A= 0, then the system has infinitely many equilibrium positions. They all lie on the line y = 2x.
Having constructed the corresponding characteristic equation l 2 + 5l = 0,
find its roots l 1 =0, l 2 =5. Therefore, all equilibrium positions are stable in the sense of Lyapunov.
Let us construct the equations for the remaining phase trajectories:
Thus, the phase trajectories lie on straight lines
Example 7
Assuming in (1) the introduction n=1, we obtain a first-order ODE that is not resolved with respect to the derivative
If (1) is uniquely solvable with respect to or, then the equations
let's call ODEs of the first order resolved with respect to derivatives.
Given that, equations (2) and (3) can be written as:
which we will call first-order ODEs in differentials.
The general form of the first order ODE in differentials will be
where the functions M(x,y) and N(x,y) are defined in some open simply connected domain, and the points at which the functions M(x,y) and N(x,y) vanish simultaneously are called singular points of the equation (6).
First-order ODEs resolved with respect to the derivative
Consider the equation
Cauchy ordinary differential equation
where the function is defined and continuous in the domain, and the domain in is considered to be a connected open set, i.e.: a) any two points of this set G can be connected by a polyline that belongs entirely to G; b) any point M of the set G is contained in G together with some neighborhood of the point M.
Definition 1. A solution to ODE (2) on the interval PrxG is a function that satisfies the conditions:
Through we will denote a connected set on the real axis, which is one of the intervals:
Example. Differential equation
has solutions on the interval of continuity of the function the entire set of antiderivatives
It follows from Example 1 that ODE (2) can have an infinite number of solutions, and this situation is general. To select a specific solution, it is necessary to set additional conditions that distinguish this solution from the entire set of solutions. These conditions are the initial conditions:
The numbers are called initial data, and the problem of finding a solution to ODE (2) that satisfies the initial conditions (9) is called the Cauchy problem for ODE (2) with initial conditions (9).
Along with the Cauchy problem for ODE (2) with initial conditions (9), we consider the integral equation
Definition 2. A function defined on an interval is called a solution to equation (10) if the following conditions are met:
- 1) - continuous);
Theorem 1. A function is a solution to the Cauchy problem for ODE (2) with initial conditions (9) if and only if it is a solution to integral equation (10).
Proof. Let be a solution of the Cauchy problem for ODE (2) with initial conditions (9). Then condition 3) of Definition 1 is satisfied, and Conditions 1) and 2) of Definition 2 follow from conditions 1) and 2) of Definition 1. Integrating the identity in 3) of Definition 1 in the range from to, we obtain condition 3) of Definition 2). Therefore, is a solution to equation (10).
Let now be the solution of the integral equation (10) on the interval. By virtue of the identity in 3) of Definition 2, the function is differentiable for and This shows that the function satisfies the initial conditions (9) and condition 1) of Definition 1. Condition 2) of Definition 1 coincides with condition 2) of Definition 2. Differentiating with respect to x, the identity in 3 ) of Definition 2, we obtain an identity in 3) of Definition 1. It follows that the function is a solution of ODE (2) with initial conditions (9). The theorem has been proven.
Definition 3. We will say that a solution to the Cauchy problem for ODE (2) with initial conditions (9) exists if such an interval exists and there exists such a solution defined on this interval and satisfying the condition
Theorem 2. If the function in equation (2) is continuous in a domain, then a solution (at least one) of the Cauchy problem for ODE (2) with initial conditions (9) exists.
Definition 4. The graph of the solution of equation (2) in the XOY plane is called the integral curve of ODE (2).
Consider the space and assign to each point from the region G a sufficiently small segment of a straight line with a slope passing through the point. The resulting family of segments in the domain G will be called the field of directions defined by ODE (2). It follows from the definition of the solution and the integral curve of ODE (2) that a curve lying in the domain G is an integral curve of this equation if and only if it is smooth and the tangent at each of its points coincides with the direction of the field at this point. From here we obtain an approximate method for constructing integral curves of ODE (2) in the region G. For the convenience of this construction, a set of points in the region G with the same field slope is found.
Definition 5. An isocline of ODE (2) in the domain G is a curve, at all points of which the direction of the field defined by ODE (2) is the same.
It follows from this definition that the set of isoclines of the ODE (2) is given by the equation where C takes admissible real values. Having constructed a grid of isoclines, we can approximately construct the integral curves of equation (2) in the domain G. We also note that the isoclines and are called the isoclines of zero and infinity, respectively, that is, at the points of the first, the direction of the field is parallel to the OX axis, and at the points of the second, parallel to the OY axis.
Example. Approximately construct an integral curve of the equation passing through the origin. The isoclines of this equation will be circles. Assuming, we obtain circles with the direction of the field. Using this grid of isoclines, we construct an approximately integral curve. Note that the solutions given equation cannot be found in the form of integrals, so the isocline method is the most appropriate.
As can be seen from the formulation of Peano's theorem, it solves the local problem of the existence of a solution to ODE (2) passing through a point. What will happen outside the interval? To solve this problem, we introduce the concept of continuation of solutions to ODE (2).
Definition 6. We will say that the solution of ODE (2), defined on the interval PrxG, is extendable to the right (left), if there is a solution to this equation, defined on the interval PrxG, (RxG,), whose restriction to coincides with The solution to ODE (2) is called in this case, by continuing the solution to the right (left).
Theorem 3. If the solution to ODE (2) is defined on the interval, then it is extendable to the right (left).
Proof. Let us show the extension to the right. We pose the Cauchy problem for ODE (2) with initial data. According to the Peano theorem, there is an interval PrxG in which there is a solution that satisfies the initial condition. Set
It is easy to check what is the solution of ODE (2) on the interval, which means that the solution, according to Definition 6, will be extendable to the right. The extension to the left is proved similarly. The theorem has been proven.
Definition 7. A solution to ODE (2) is called complete if it is extendable neither to the right nor to the left.
It follows from Definition 7 and Theorem 3 that the domain of the complete solution is always the open interval PrxG, which is called the maximum interval for the existence of a solution to ODE (2).
Let us turn again to Peano's theorem and note that it asserts the existence of at least one, and not necessarily unique, solution of the Cauchy problem for ODE (2) with initial conditions (9). This implies the need to introduce the concepts of a point and a domain of uniqueness of ODE (2).
Definition 8. A point is called a point of uniqueness of ODE (2) if there exists a -neighborhood of this point such that one and only one integral curve of ODE (2) passes through the point inside.
The region consisting entirely of points of uniqueness of ODE (2) will be called the region of uniqueness of equation (2). This implies that any two solutions of ODE (2) from the domain D that coincide at some point coincide everywhere in the domain D.
Definition 9. A point is called a point of non-uniqueness of ODE (2) if more than one integral curve of ODE (2) passes through it in any neighborhood of this point.
Definition 10. A solution to ODE (2) is called partial if each point corresponding to this solution of the integral curve is a point of uniqueness to ODE (2).
The whole set of particular solutions of ODE (2) in the domain D is called the general solution of equation (2) in this domain.
Definition 11. The function where is an arbitrary constant is called the general solution of ODE (2) in the domain of uniqueness D if:
1) for any point, the equation is uniquely solvable with respect to C, that is,
2) the function is a solution of the Cauchy problem for ODE (2) with initial data
In the general case, integrating ODE (2) in the domain, we obtain the general solution in an implicit form:
Definition 12. The general solution of ODE (2) written in the form (11) or (12) is called common integral, and the function is the integral of ODE (2) in the domain D.
The main properties of the ODE integral (2) are:
- 1) the integral retains a constant value along any integral curve of equation (2) located in the region D;
- 2) identity holds for all
Example. For the equation, the function is an integral, and the function, where C is an arbitrary constant, is a general solution in the region. Check for yourself the feasibility of properties 1) and 2) of the integral in this case.
Theorem 4 (Cauchy-Picard). If the function is continuous together with its partial derivative in the domain, then there is a unique solution of the Cauchy ODE problem (2) with initial conditions, i.e., there is a unique integral curve of equation (2) that belongs entirely to the domain D and passes through the point
Definition 13. A solution to ODE (2) is called special if each point corresponding to this solution of the integral curve is a non-unique point of ODE (2). The integral curve corresponding to a special solution is called the special integral curve of ODE (2).
It follows from the definitions of the domain of uniqueness and the singular solution that singular solutions can only be the boundary of the domain D. It follows from the Cauchy-Picard theorem that at each point of a singular integral curve, at least one of the conditions of this theorem is violated.
Example. The equation has a general solution in the regions and the plane. Moreover, the function is also a solution to this equation, but it cannot be obtained from the general solution for any value of an arbitrary constant C. Since the derivative does not exist at the points of the integral curve, it is a special integral curve.
Construct a sketch of the location of the integral curves on the XOY plane and indicate the points of uniqueness and non-uniqueness.
Consider a one-parameter (C is a parameter) family
smooth curves completely filling the region and assume that the function is differentiable with respect to the variables x and y in this region.
Let us pose the following problem: to compose a first-order ODE in the domain D, for which each curve of the given family is an integral curve. Obviously, to solve the problem, it is sufficient to exclude the parameter C from the system of equations
Example. Let families of smooth curves be given and where C is a parameter. Let us construct the corresponding ODEs of the first order, excluding the parameter С from the systems
We obtain, respectively, the ODE and
Integration of the simplest first-order ODEs resolved with respect to the derivative
Definition 14. If the general solutions of equations (2)-(6) can be found in the form of a finite combination of integration operations, then we will say that the solution is found in quadratures.
Note that in some cases the left side of the equation
is the total differential of some function, that is
Then the general integral of ODE (6) will be the relation
where C is an arbitrary constant.
If, on multiplying both parts of ODE (6) by some function continuous in the domain G of continuity of functions, and the left side of the resulting equation
becomes a total differential of some function, then the relation
where C is an arbitrary constant, is the general integral of the ODE (6).
Separable Variable Equations
General view of the first order ODE with separable variables
where the functions and - are continuous in the interval and the functions and - are continuous in the interval
Let's consider the area. In this region, apart from the singular points at which the functions and vanish simultaneously, Eq. (19) generally has two types of solutions:
1) , if and, (20)
if and (21)
2) consider the area in which. Equation (19) in the region D is equivalent to the equation
the left side of which is the total differential of the function
Then the general integral of ODE (19) in the domain D will be the relation
where C is an arbitrary constant.
Thus, the entire set of solutions to ODE (19) consists of solutions (20), (21), and (24). Note that among the solutions (20) and (21) there may be special ones, and the integral line will be a special integral line of ODE (19), if one of the integrals
where small enough, is convergent. Similarly, the integral line will be a singular integral line of ODE (19) if and one of the integrals
is convergent.
Another method for finding singular solutions for a first-order ODE that is allowed with respect to the derivative is related to the violation of the conditions of the Cauchy-Picard theorem at the points of the solutions under study.
Let us indicate one more way of recognizing special solutions for ODE (19). If solutions (20) and (21) are not obtained from (24) for any particular numerical values of C, then they are special solutions to ODE (19).
The solution of the Cauchy problem for ODE (19) with initial data has the form:
Example. The entire set of solutions to the equation consists of a straight line and a set of curves, where C is an arbitrary constant. When deciding whether the integral line will be a special integral line, we turn to the Cauchy-Picard theorem. Since the function and its derivative are continuous in the domain, the line is not a special integral line.
Homogeneous ODEs of the first order
The equation
is called homogeneous if the functions and are homogeneous functions in the variables x and y of the same order, i.e.
By substitution, equation (6) is reduced to an equation with separable variables
Functions where the roots of the equation
are solutions of Eq. (6), and some of them may contain singular ones. Axle axles can be special:
Note that the point (0,0) is singular for ODE (6).
The equation
is always reduced to a homogeneous equation or to an equation with separable variables, and:
a) if then ODE (31) is homogeneous;
b) if or and, then after a linear substitution we obtain an equation with separable variables;
c) if or and then the system of equations
has a unique solution. The change of ODE (31) leads to a homogeneous equation
If equation (6) is not homogeneous, but after changing where it turns into the equation
where the functions and are homogeneous, then ODE (6) is called in this case a generalized homogeneous equation.
Example. The equation is a generalized homogeneous equation, since after the replacement it turns into the equation
which for will be a homogeneous equation.
Linear ODEs of the first order
General view of a first-order linear ODE
where the functions are continuous on the interval and
We divide both sides of equation (32) by a function and obtain the equivalent equation
where are also continuous on
Multiplying both parts of the ODE (33) by the function, we obtain the ODE
integrating which we obtain the general integral
and general solution
where C is an arbitrary constant.
Solving ODE (33) with respect to the derivative and applying the Cauchy-Picard theorem, we obtain that equation (33) has no special solutions.
Example. Linear Equation
after multiplying both parts by the function, it is transformed into the equation
Hence we obtain the general solution of ODE (37)
where C is an arbitrary constant.
Some first-order ODEs become linear if x is considered the desired function and y is the argument, i.e.
where - continuous functions on the interval
Example. J.Bernoulli's equation is transformed by substitution into a linear one, the general solution of which has the form. From here and from the substitution we obtain the general solution of the desired equation where C is an arbitrary constant.
J. Riccati's equation
where are functions continuous in the interval, in the general case it is not solved in quadratures. However, if at least one of its particular solutions is known, then by substitution it is reduced to the J. Bernoulli equation.
ODEs of the first order in total differentials
Definition 15. ODE of the first order
where the functions are continuous in the domain and the left side is the total differential of some differentiable function, let's call it an ODE in total differentials. The general integral of this equation is given by the relation
Let us pose the following questions: 1. How can one determine, by the form of equation (6), whether it is an equation in total differentials? 2. How to find the function?
These questions are answered by the following
Theorem 5. Equation (6), where - are continuous in the domain if and only then there will be an ODE in total differentials if the equality holds for all
Proof. Let us prove the first part of the theorem. Let the left side of equation (6) be the total differential of the function in the region, that is, the identity
From identity (45), due to the arbitrariness of u, we obtain the identities
differentiating which, respectively, with respect to y and x, we obtain
Due to the continuity of functions in the domain, from the theorem on the equality of mixed derivatives we obtain equality (44).
Let us prove the second part of the theorem. Let equality (44) hold in the domain. It is required to show that there exists such a function whose total differential in the domain is identically equal to the left side of equation (6), that is, identities (46) hold. The whole set of functions satisfying the first identity in (46) is given by the formula
where is an arbitrary differentiable function of y. Let us show what can be chosen so that the second identity in (46) also holds, i.e.,
Substituting (48) into (49), we obtain
It follows from this that it can be found if the right side in (50) does not depend on the variable x, that is,
Transforming the left side into (51), we obtain
(by condition (44) of theorem (5)) = 0.
Thus, the function, up to an arbitrary constant, has the form
The theorem has been proven.
Comment. The proof of the second part of Theorem 5 implies a practical method for solving ODEs in total differentials.
Example. The equation is an ODE in total differentials, since in the region Here By formula (52) we find the function
where is an arbitrary constant. Thus, the general integral of the desired equation will be where is an arbitrary constant.
Theorem 6. Solution of the Cauchy problem for ODE(6) in total differentials with initial conditions where the domain D is the domain of continuity of functions and in which one of the formulas is given:
and this solution is unique.
Note that at the points where the functions and and are called singular for ODE (6), the uniqueness of the solution of the Cauchy problem is not guaranteed.
Integrating factor
Consider equation (6) with functions continuous in the domain, which is not an equation in total differentials.
Definition 16. If, when both parts of ODE (6) are multiplied by a function that is continuous together with its partial derivatives and nonzero in the domain, equation (6) turns into an ODE in total differentials, that is, for all, the equality
then we call the function the integrating factor of ODE (6), and equation (55) - the equation of the integrating factor.
Note that equation (55) is no easier to solve than equation (6), so consider the cases where the integrating factor is found quite easily:
1) let, then equation (55) takes the form
If the right side in (56) depends only on x, that is,
then from (56) by integration we find up to a multiplicative constant
where is an arbitrary constant (usually assumed).
2) let, then from (55) we obtain
If the right side in (58) is a function of y alone, that is
then integrating (58) we obtain
where is an arbitrary constant (for convenience, it is usually considered)
3) let, where is a known function, continuous together with its partial derivatives in the region, then from (55) we obtain
If the right side in (60) is a function then
where is an arbitrary constant.
Let us clarify the question of the existence of an integrating factor in the general case. fair
Theorem 7. If ODE (6) has a common integral in the domain of uniqueness
where the function is continuous in the domain together with its partial derivatives up to the second order inclusive, then equation (6) has an integrating factor.
Proof. Since (62) is the general integral of ODE (6), the identity holds for all
Let us assume that in the region then from (63) we obtain an integrating factor in the form
Indeed, multiplying both parts of ODE (6) by function (64), we obtain
i.e. function (64) is an integrating factor of ODE (6). The theorem has been proven.
Consequence. It follows from Theorem 7 that each integrating factor of ODE (6) corresponds, up to an additive constant, to the integral of this equation.
It can be seen from the equation of the integrating factor (55) that there are an infinite number of integrating factors, if there is at least one), different from the identical constant, for example, where is a constant number Does this set exhaust the entire set of integrating factors of equation (6)? How to find the whole set of integrating factors if at least one is known? These questions are answered by the following theorems.
Theorem 8. If is the integrating factor of ODE (6), and is the corresponding first integral of this equation, then the function
where is an arbitrary differentiable function, will also be an integrating factor of equation (6).
Proof. We multiply the left side of the ODE (6) by the function (65)
This and Definition 16 imply the assertion of the theorem. The theorem has been proven.
Theorem 9. If is the integrating factor of ODE (6), and is the corresponding integral of equation (6), then any integrating factor of this equation is found by the formula
where is an arbitrary differentiable function.
Proof. Let - some integrating factor of ODE (6), and - the integral corresponding to it. Then we have the equalities
from which we get
Let us divide term by term the equalities of the first row by the corresponding equalities of the second row, we get
Thus, the Jacobian of the functions and is identically equal to zero. Hence, according to the theorem from the analysis of the dependence of functions, there is a functional dependence between these functions
where is a twice continuously differentiable function. Let us divide term by term the second equality in (67) by the first, we obtain
Differentiating (68) and substituting into (69), we have
Q.E.D.
Consequence. If and are two different integrating factors of ODE (6), the ratio of which is not identically constant, then the expression
where is an arbitrary constant, is the general integral of the ODE (6).
Example. Functions are integrating factors of the equation, and their ratio is not identically constant. So the expression
where is an arbitrary constant, is the general integral of this equation.
Proof of the Cauchy-Picard theorem on the existence and uniqueness of a solution to the Cauchy problem for a first-order ODE solved with respect to the derivative
Consider the equation
where the function - is continuous in some simply connected open domain, and we set the Cauchy problem for it with the initial conditions
Theorem 10 (Cauchy-Picard). If the function is continuous in a rectangle
together with its partial derivative, then there is this:
Differential equations of the second order and higher orders.
Linear DE of the second order with constant coefficients.
Solution examples.
We pass to the consideration of differential equations of the second order and differential equations of higher orders. If you have a vague idea of what a differential equation is (or don’t understand what it is at all), then I recommend starting with the lesson First order differential equations. Solution examples. Many solution principles and basic concepts of first-order diffurs are automatically extended to higher-order differential equations, so it is very important to first understand the first order equations.
Many readers may have a prejudice that DE of the 2nd, 3rd, and other orders is something very difficult and inaccessible for mastering. This is not true . Learning to solve higher-order diffuses is hardly more difficult than “ordinary” 1st-order DEs. And in some places it is even easier, since the material of the school curriculum is actively used in the decisions.
Most Popular second order differential equations. Into a second order differential equation necessarily includes the second derivative and not included
It should be noted that some of the babies (and even all at once) may be missing from the equation, it is important that the father was at home. The most primitive second-order differential equation looks like this:
Third order differential equations in practical tasks occur much less frequently, according to my subjective observations in State Duma they would get about 3-4% of the votes.
Into a third order differential equation necessarily includes the third derivative and not included derivatives of higher orders:
The simplest differential equation of the third order looks like this: - dad is at home, all the children are out for a walk.
Similarly, differential equations of the 4th, 5th and higher orders can be defined. In practical problems, such DE slips extremely rarely, however, I will try to give relevant examples.
Higher order differential equations that are proposed in practical problems can be divided into two main groups.
1) The first group - the so-called lower-order equations. Fly in!
2) The second group - higher-order linear equations with constant coefficients. Which we will begin to consider right now.
Second Order Linear Differential Equations
with constant coefficients
In theory and practice, two types of such equations are distinguished - homogeneous equation and inhomogeneous equation.
Homogeneous DE of the second order with constant coefficients has the following form:
, where and are constants (numbers), and on the right side - strictly zero.
As you can see, there are no special difficulties with homogeneous equations, the main thing is that solve the quadratic equation correctly.
Sometimes there are non-standard homogeneous equations, for example, an equation in the form , where at the second derivative there is some constant , different from unity (and, of course, different from zero). The solution algorithm does not change at all, one should calmly compose the characteristic equation and find its roots. If the characteristic equation will have two different real roots, for example: , then the general solution can be written in the usual way: .
In some cases, due to a typo in the condition, “bad” roots can turn out, something like . What to do, the answer will have to be written like this:
With "bad" conjugate complex roots like no problem either, general solution:
That is, a general solution exists in any case. Because any quadratic equation has two roots.
In the final paragraph, as I promised, we will briefly consider:
Higher Order Linear Homogeneous Equations
Everything is very, very similar.
The linear homogeneous equation of the third order has the following form:
, where are constants.
For this equation, you also need to compose a characteristic equation and find its roots. The characteristic equation, as many have guessed, looks like this:
, and it anyway It has exactly three root.
Let, for example, all roots be real and distinct: , then the general solution can be written as follows:
If one root is real, and the other two are conjugate complex, then we write the general solution as follows:
A special case is when all three roots are multiples (the same). Let's consider the simplest homogeneous DE of the 3rd order with a lonely father: . The characteristic equation has three coincident zero roots. We write the general solution as follows:
If the characteristic equation has, for example, three multiple roots, then the general solution, respectively, is:
Example 9
Solve a homogeneous differential equation of the third order
Solution: We compose and solve the characteristic equation:
, - one real root and two conjugate complex roots are obtained.
Answer: common decision
Similarly, we can consider a linear homogeneous fourth-order equation with constant coefficients: , where are constants.
Department of Physical Chemistry SFedU (RSU)
NUMERICAL METHODS AND PROGRAMMING
Materials for the lecture course
Lecturer - Art. teacher 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 any dynamic system. This is best done in the form of differential equations ( DU) or systems of differential equations. Most often, such a 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, the differential equation can be converted to a form in which the highest derivative is expressed explicitly. This form of writing is called an equation. allowed with respect to the highest derivative(at the same time, the highest derivative is absent on the right side of the equation):
It is this form of notation 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
By solving the ordinary differential equation a function y(x) is called such that for any x it satisfies this equation in a certain finite or infinite interval. The process of solving a differential equation is called integration of the differential equation.
General solution of ODE nth order contains n arbitrary constants C 1 , C 2 , …, C n
This obviously follows from the fact that indefinite integral equal to the antiderivative of the integrand plus the integration constant
Since it is necessary to carry out n integrations to solve the n-th order DE, then in common decision n integration constants appear.
Private decision The ODE is obtained from the general one, if the integration constants are given some values by defining some additional conditions, the number of which allows us to calculate all the indefinite integration constants.
Exact (analytical) solution (general or particular) of a differential equation implies obtaining the desired solution (of the function y(x)) as an expression from elementary functions. This is far from always possible even for first-order equations.
Numerical solution DE (private) is to calculate the function y(x) and its derivatives in some given points lying on a certain interval. That is, in fact, the solution of the DE of the nth order 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 be two columns - x and y.
The set of abscissa values in which the value of the function is determined is called grid, on which the function y(x) is defined. The coordinates themselves are called grid nodes. Most often, for convenience, uniform grids, in which the difference between neighboring nodes is constant and is called grid step or integration step differential equation
Or , i= 1, …, N
For determining private decision it is necessary to set additional conditions that will allow one to calculate the integration constants. Moreover, there must be exactly n such conditions. For equations of the first order - one, for the second - 2, etc. Depending on the way they are specified in solving differential equations, there are three types of problems:
· Cauchy problem (initial problem): It is necessary to find such private solution differential equation that satisfies certain initial conditions given at one point:
that is, given 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 that point. This point (x 0) is called primary. For example, if the DE of the 1st order is being solved, then the initial conditions are expressed as a pair of numbers (x 0 , y 0)
This kind of problem occurs when 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, one can also cite the problem of heat transfer or mass transfer (diffusion), the equation of motion material point under the action of forces, etc.
· Boundary 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 of the differential equation between these points. The additional conditions themselves in this case are called regional (boundary) conditions. Naturally, the boundary value problem can be solved for ODEs of at least order 2. Below is an example of a second-order ODE with boundary conditions (values of the function are given at two different points):
· Sturm-Liouville problem (problem for eigenvalues). Problems of this type are similar to a boundary value problem. When solving them, it is necessary to find for what values of any parameter the solution DU satisfies the boundary conditions (eigenvalues) and the functions that are the 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 ODEs
Consider some numerical methods for solving Cauchy problems(initial problem) ordinary differential equations of the first order. We write this equation in a 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 the given grid points if the initial values are known, where is the value of the function y(x) at the initial point x 0 .
We transform the equation by multiplying by d x
And let's integrate the left and right parts between the i -th and i + 1st grid nodes.
(6.3)
We have obtained an expression for constructing a solution at the i + 1 node of integration through the values of x and y at the i -th node of the grid. The difficulty, however, lies in the fact that the integral on the right-hand side is an integral of an implicitly given function, which cannot be found analytically in the general case. Numerical methods for solving ODEs in a different way approximate (approximate) the value of this integral to construct formulas for the numerical integration of the ODE.
From the set of methods developed for solving ODEs of the first order, we consider the methods , and . They are quite simple and give an initial idea of the approaches to solving a given problem within the framework of a numerical solution.
Euler method
Historically the first and most in a simple way The numerical solution of the Cauchy problem for ODEs of the first order 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 required value of the function at the point x i+1 .
If we now transform this equation and take into account the uniformity of the integration grid, we get an iterative formula by which we can calculate yi+1, if y i is known at the point x i :
Comparing the Euler formula with the general expression obtained earlier, it can be seen that for the approximate calculation of the integral in the Euler method uses the simplest integration formula - the formula of rectangles along the left edge of the segment.
The graphical interpretation of the Euler method is also not difficult (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 slope of the tangent 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
where does Euler's formula come from. Thus, the essence of the Euler method is to replace the function y(x) on the integration segment with a straight line tangent to the graph at the point x=x i . If the desired function is very different from the linear one on the integration interval, then the calculation error will be significant. The Euler method error is directly proportional to the integration step:
Error~ h
The calculation process is constructed as follows. For given initial conditions x0 and y 0 can be calculated
Thus, a table of values of the function y(x) is built with a certain step ( h) on x on the segment. Error in determining the value y(x i) in this case, it will be the smaller, the smaller the step length is chosen h(which is determined by the accuracy of the integration formula).
For large h, the Euler method is very inaccurate. It gives an increasingly accurate approximation as the integration step decreases. If the segment is too large, then each segment is divided into N integration segments and the Euler formula with a step is applied to each of them, that is, the integration step h is taken less than the grid step on which the solution is determined.
Example:
Using the Euler 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 the 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 (N=1 ) will be calculated for each grid node. For the first four grid nodes, the calculations will be as follows:
Full results (up to the fifth decimal place) are given in the third column - h = 0.1 (N = 1). In the second column of the table, for comparison, the values calculated from the analytical solution of this equation are given .
The second part of the table shows the relative error of the obtained solutions. It can be seen that for 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 | Exact 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 the calculated values of the function for various 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 for 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). However, you can do otherwise and approximately calculate the value of the function in the node i+1 using the usual formula:
,
which is then used in the calculation according to (6.6).
Thus, a method is obtained Gyuna 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 Gün method is already proportional to the square of the integration step.
Error~h2
The approach used in Gün's method is used to build the so-called methods forecast and correction, which will be discussed later.
Example:
We will carry out calculations for the equation () using the Gün method.
With the integration step h = 0.1 at the first node of the grid x 1 we get:
Which is much more accurate than the value obtained by the Euler method with the same integration step. Table 2 below shows the comparative results of calculations for h = 0.1 by the Euler and Gün methods.
Table 2 Solution of the equation by the Euler and Gün methods
x | Exact | Gün 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% |
We note a significant increase in the computational accuracy of the Gün method compared to the Euler method. So, for the node x =0.1, the relative deviation of the value of the function determined by the Gün method turns out to be 30 (!) times less. The same accuracy of calculations by the Euler formula is achieved when the number of integration segments N is approximately 30. Therefore, when using the Gün method with the same calculation accuracy, it will take about 15 times less computer time than when using the Euler method.
Checking the stability of the solution
The solution of 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 an integration step h and with a reduced (for example, two) step size
As a stability criterion, one can use the smallness of the relative change in the obtained solution with a decrease in the integration step (ε is a preassigned small value)
Such a check can also be carried out for all solutions on 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, and so on. until a stable solution is obtained.
Runge-Kutta Methods
Further improvement in the accuracy of solving the 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 integration using the rectangle formula () to using the trapezoid formula () when approximating this integral.
Using the well-established Simpson formula , one can obtain an even more accurate formula for solving the Cauchy problem for first-order ODEs - the Runge-Kutta method widely used in computational practice.
The advantage of the multi-step Adams methods in solving ODEs is that at each node only one value of the right 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 for calculations using a k-step formula, it is necessary to know the value of the function at k nodes. Therefore, it is necessary to obtain (k-1) solution at the first nodes x 1 , x 2 , …, x k-1 using some one-step method, for example, the method