# A numerical solution of problem for essentially loaded differential equations with an integro-multipoint condition

• Zhazira M. Kadirbayeva and Symbat S. Kabdrakhova
From the journal Open Mathematics

## Abstract

We study a linear boundary value problem for systems of essentially loaded differential equations with an integro-multipoint condition. We make use of the numerical implementation of the Dzhumabaev parametrization method to obtain the desired result, which is well supported by two numerical examples.

MSC 2010: 34B10; 45J05; 65L06

## 1 Introduction and problem statement

In the last few decades, loaded differential equations with boundary conditions have been studied by many researchers [1,2,3, 4,5,6, 7,8,9, 10,11]. This is because, loaded equations describe problems in optimal control, regulation of the layer of soil water and ground moisture, and underground fluid and gas dynamics. Monographs [12,13] contain various applications of loaded equations as a method for studying problems in mathematical biology, mathematical physics, the theory of mathematical modeling of nonlocal processes and phenomena, and the theory of elastic shells.

Various problems for loaded differential equations with integral conditions and methods for finding their solutions are considered in [14,15,16, 17,18,19, 20,21]. Boundary value problems with integro-multipoint boundary conditions have been studied by many researchers, for example, see [22,23].

In a recent article [24], the author discussed the numerical method for solving the boundary value problem for essentially loaded differential equations based on the Dzhumabaev parametrization method [25,26]. This is a constructive method originally developed to investigate and solve boundary value problems for ordinary differential equations. The Dzhumabaev parametrization method is based on dividing the interval [ 0 , T ] into N parts and introducing the additional parameters. In [25], coefficient criteria were established for the unique solvability of linear boundary value problems. An algorithm for finding their approximate solutions was developed. The Dzhumabaev parametrization method was later extended to boundary value problems, both linear and nonlinear, for various classes of equations. In [24], the term essentially loaded differential equation means that the right side of the differential equation depends on the value of the desired solution and its derivatives at given points, where the order of the derivatives is not less than the order of the differential part of the equation. In [24], by assuming the invertibility of the matrix compiled through the coefficients at the values of the derivative of the desired function at load points, we reduce the considering problem to a two-point boundary value problem for loaded differential equations.

Motivated by the research going on in this direction, in this article, we study the finding a numerical solution to the linear boundary value problem for systems of essentially loaded differential equations with an integro-multipoint condition of the form:

(1) d x d t = A 0 ( t ) x + i = 1 m A i ( t ) x ˙ ( θ i ) + f ( t ) , t ( 0 , T ) ,

(2) j = 0 m + 1 B j x ( θ j ) + θ 0 θ m + 1 C ( t ) x ( t ) d t = d , d R n , x R n .

Here, ( n × n ) -matrices A i ( t ) , ( i = 0 , m ¯ ) , C ( t ) , and n -vector-function f ( t ) are continuous on [ 0 , T ] , B j ( j = 0 , m + 1 ¯ ) are constant ( n × n ) -matrices, and 0 = θ 0 < θ 1 < < θ m < θ m + 1 = T ; x = max i = 1 , n ¯ x i .

Let C ( [ 0 , T ] , R n ) denote the space of continuous functions x : [ 0 , T ] R n with the norm x 1 = max t [ 0 , T ] x ( t ) .

A solution to problem (1), (2) is a continuously differentiable on ( 0 , T ) function x ( t ) C ( [ 0 , T ] , R n ) satisfying the system of essentially loaded differential equations (1) and the integro-multipoint condition (2).

The aim of this article is to propose a numerical implementation of the Dzhumabaev parametrization method for solving the boundary value problem for systems of essentially loaded differential equations with integro-multipoint condition (1), (2).

### Definition 1

Problem (1), (2) is called uniquely solvable, if for any function f ( t ) C ( [ 0 , T ] , R n ) and vector d R n , it has a unique solution.

The interval [ 0 , T ] is partitioned by the loading points: [ 0 , T ) = r = 1 m + 1 [ θ r 1 , θ r ) .

Let C ( [ 0 , T ] , θ , R n ( m + 1 ) ) denote the space of function systems x [ t ] = ( x 1 ( t ) , x 2 ( t ) , , x m + 1 ( t ) ) , where x r : [ θ r 1 , θ r ) R n are continuous on [ θ r 1 , θ r ) and have finite left-sided limits lim t θ r 0 x r ( t ) for all r = 1 , m + 1 ¯ , with the norm x [ ] 2 = max r = 1 , m + 1 ¯ sup t [ θ r 1 , θ r ) x r ( t ) .

The restriction of the function x ( t ) to the r th subinterval [ θ r 1 , θ r ) is denoted by x r ( t ) , i.e., x r ( t ) = x ( t ) for t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ , and the restriction of the function x ˙ ( t ) to the r th subinterval [ θ r 1 , θ r ) is denoted by x ˙ r ( t ) , i.e., x ˙ r ( t ) = x ˙ ( t ) for t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ . The problem (1), (2) is then transformed into the equivalent problem:

(3) d x r d t = A 0 ( t ) x r + i = 1 m A i ( t ) x ˙ i + 1 ( θ i ) + f ( t ) , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ ,

(4) j = 0 m B j x j + 1 ( θ j ) + B m + 1 lim t θ m + 1 0 x m + 1 ( t ) + r = 1 m + 1 θ r 1 θ r C ( t ) x r ( t ) d t = d ,

(5) lim t θ i 0 x i ( t ) = x i + 1 ( θ i ) , i = 1 , m ¯ ,

where (5) are conditions for matching the solution at the interior partition points. From conditions (5) and the assumption of the continuity of the coefficients A i ( t ) , ( i = 0 , m ¯ ) , it follows that the derivatives of the solution will also be continuous at the partition points.

A solution of problem (3)–(5) is a system of functions x [ t ] = ( x 1 ( t ) , x 2 ( t ) , , x m + 1 ( t ) ) C ( [ 0 , T ] , θ , R n ( m + 1 ) ) , where functions x r ( t ) , r = 1 , m + 1 ¯ , are continuously differentiable on [ θ r 1 , θ r ) , which satisfies system (3) and conditions (4) and (5).

Problem (1), (2), and (3)–(5) are equivalent. If a system of functions x ˜ [ t ] = ( x ˜ 1 ( t ) , x ˜ 2 ( t ) , , x ˜ m + 1 ( t ) ) C ( [ 0 , T ] , θ , R n ( m + 1 ) ) is a solution of problem (3)–(5), then the function x ˜ ( t ) defined by the equalities x ˜ ( t ) = x ˜ r ( t ) , x ˜ ˙ ( t ) = x ˜ ˙ r ( t ) , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ , and x ˜ ( θ m + 1 ) = lim t θ m + 1 0 x ˜ m + 1 ( t ) is a solution of the original problem (1), (2). Conversely, if x ( t ) is a solution of problem (1), (2), then the system of functions x [ t ] = ( x 1 ( t ) , x 2 ( t ) , , x m + 1 ( t ) ) , where x r ( t ) = x ( t ) , x ˙ r ( t ) = x ˙ ( t ) , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ , and lim t θ m + 1 0 x m + 1 ( t ) = x ( θ m + 1 ) is a solution of problem (3)–(5).

Introducing the additional parameters λ r = lim t θ r 1 + 0 x r ( t ) , μ r = lim t θ r 1 + 0 x ˙ r ( t ) , r = 1 , m + 1 ¯ , and performing a replacement of the function u r ( t ) = x r ( t ) λ r on the each interval [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ , we obtain the boundary value problem with parameters λ r , r = 1 , m + 1 ¯ :

(6) d u r d t = A 0 ( t ) ( u r + λ r ) + i = 1 m A i ( t ) μ i + 1 + f ( t ) , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ ,

(7) u r ( θ r 1 ) = 0 , r = 1 , m + 1 ¯ ,

(8) j = 0 m B j λ j + 1 + B m + 1 λ m + 1 + B m + 1 lim t θ m + 1 0 u m + 1 ( t ) + r = 1 m + 1 θ r 1 θ r C ( t ) ( u r ( t ) + λ r ) d t = d ,

(9) λ i + lim t θ i 0 u i ( t ) = λ i + 1 , i = 1 , m ¯ ,

(10) lim t θ r 1 + 0 u ˙ r ( t ) = μ r , r = 1 , m + 1 ¯ .

A solution to problem (6)–(10) is a triple ( λ , μ , u [ t ] ) , with elements λ = ( λ 1 , λ 2 , , λ m + 1 ) R n ( m + 1 ) , μ = ( μ 1 , μ 2 , , μ m + 1 ) R n ( m + 1 ) , and u [ t ] = ( u 1 ( t ) , u 2 ( t ) , , u m + 1 ( t ) ) C ( [ 0 , T ] , θ , R n ( m + 1 ) ) , where u r ( t ) are continuously differentiable on [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ , and satisfying the system of ordinary differential equations (6) and conditions (7)–(10) at the λ r = λ r , μ r = μ r , j = 1 , m + 1 ¯ .

Problem (1), (2) are equivalent to problem (6)–(10). If the function x ( t ) is a solution to problem (1), (2), then the triple ( λ , μ , u [ t ] ) , where λ = ( x ( θ 0 ) , x ( θ 1 ) , , x ( θ m ) ) , μ = ( x ˙ ( θ 0 ) , x ˙ ( θ 1 ) , , x ˙ ( θ m ) ) , and u [ t ] = ( x ( t ) x ( θ 0 ) , x ( t ) x ( θ 1 ) , , x ( t ) x ( θ m ) ) is a solution to problem (6)–(10). Conversely, if the triple ( λ ˜ , μ ˜ , u ˜ [ t ] ) , with elements λ ˜ = ( λ ˜ 1 , λ ˜ 2 , , λ ˜ m + 1 ) R n ( m + 1 ) , μ ˜ = ( μ ˜ 1 , μ ˜ 2 , , μ ˜ m + 1 ) R n ( m + 1 ) , u ˜ [ t ] = ( u ˜ 1 ( t ) , u ˜ 2 ( t ) , , u ˜ m + 1 ( t ) ) C ( [ 0 , T ] , θ , R n ( m + 1 ) ) , is a solution to problem (6)–(10), then the function x ˜ ( t ) defined by the equalities x ˜ ( t ) = u ˜ r ( t ) + λ ˜ r , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ , and x ˜ ( T ) = λ ˜ m + 1 + lim t T 0 u ˜ m + 1 ( t ) will be the solution of the original problem (1), (2).

First, find the values of μ j , j = 1 , m + 1 ¯ . By passing on the right-hand side of (6) to the limit as t θ r 1 + 0 and r = 1 , m + 1 ¯ , and substituting the expression into (10), we obtain the system of equations for the unknown parameters μ j , j = 1 , m + 1 ¯ :

(11) μ r i = 1 m A i ( θ r 1 ) μ i + 1 = A 0 ( θ r 1 ) λ r + f ( θ r 1 ) , r = 1 , m + 1 ¯ ,

that is, we can rewrite equation (11) in the following form:

(12) G ( θ ) μ = H ( θ , λ ) , μ R n ( m + 1 ) .

Here,

G ( θ ) = I A 1 ( θ 0 ) A 2 ( θ 0 ) A m 1 ( θ 0 ) A m ( θ 0 ) O I A 1 ( θ 1 ) A 2 ( θ 1 ) A m 1 ( θ 1 ) A m ( θ 1 ) O A 1 ( θ 2 ) I A 2 ( θ 2 ) A m 1 ( θ 2 ) A m ( θ 2 ) O A 1 ( θ m 1 ) A 2 ( θ m 1 ) I A m 1 ( θ m 1 ) A m ( θ m 1 ) O A 1 ( θ m ) A 2 ( θ m ) A m 1 ( θ m ) I A m ( θ m ) ,

where I is the identity matrix of dimension n and O is the zero matrix of dimension n ,

H ( θ , λ ) = ( H 1 ( θ , λ ) , H 2 ( θ , λ ) , , H m + 1 ( θ , λ ) ) , H r ( θ , λ ) = A 0 ( θ r 1 ) λ r + f ( θ r 1 ) , r = 1 , m + 1 ¯ .

Assume that the matrix G ( θ ) is invertible. Denote by S ( θ ) the inverse matrix G ( θ ) , i.e., S ( θ ) = [ G ( θ ) ] 1 , where S ( θ ) = s p , k ( θ ) , p , k = 1 , m + 1 ¯ . Then from (12), we can uniquely determine μ :

μ = [ G ( θ ) ] 1 H ( θ , λ ) = S ( θ ) H ( θ , λ ) ,

i.e.,

(13) μ r = k = 1 m + 1 s r , k ( θ ) { A 0 ( θ k 1 ) λ k + f ( θ k 1 ) } , r = 1 , m + 1 ¯ .

In (6), substituting the right-hand side of (13) instead of μ i + 1 , i = 1 , m + 1 ¯ , we obtain

(14) d u r d t = A 0 ( t ) ( u r + λ r ) + j = 1 m + 1 D j ( t ) λ j + F ( t ) , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ ,

where

D j ( t ) = i = 1 m A i ( t ) s i + 1 , j ( θ ) A 0 ( θ j 1 ) , j = 1 , m + 1 ¯ , F ( t ) = i = 1 m k = 1 m + 1 A i ( t ) s i + 1 , k ( θ ) f ( θ k 1 ) + f ( t ) .

### Definition 2

The Cauchy problem (14) and (7) is called uniquely solvable, if for any λ R n ( m + 1 ) , f ( t ) C ( [ 0 , T ] , R n ) it has a unique solution.

Let Φ r ( t ) be a fundamental matrix of the differential equation d x d t = A 0 ( t ) x on [ θ r 1 , θ r ] , r = 1 , m + 1 ¯ .

Then the unique solution to the Cauchy problem for the system of ordinary differential equations (14) and (7) with the fixed values

(15) u r ( t ) = Φ r ( t ) θ r 1 t Φ r 1 ( τ ) A 0 ( τ ) d τ λ r + Φ r ( t ) θ r 1 t Φ r 1 ( τ ) j = 1 m + 1 D j ( τ ) d τ λ j + Φ r ( t ) θ r 1 t Φ r 1 ( τ ) F ( τ ) d τ , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ .

As a rule, construction of fundamental matrices for the systems of ordinary differential equations with variable coefficients fails. Therefore, later, we propose a numerical method for solving problem (1), (2). For this purpose, we consider the Cauchy problems for ordinary differential equations on subintervals

(16) d z d t = A ( t ) z + P ( t ) , z ( θ r 1 ) = 0 , t [ θ r 1 , θ r ] , r = 1 , m + 1 ¯ ,

where P ( t ) is either ( n × n ) matrix, or n vector, both continuous on [ θ r 1 , θ r ] , r = 1 , m + 1 ¯ . Consequently, the solution to problem (16) is a square matrix or a vector of dimension n . Denote by a r ( P , t ) a solution to the Cauchy problem (16). Obviously,

(17) a r ( P , t ) = Φ r ( t ) θ r 1 t Φ r 1 ( τ ) P ( τ ) d τ , t [ θ r 1 , θ r ] , r = 1 , m + 1 ¯ ,

where Φ r ( t ) is a fundamental matrix of the differential equation (16) on the r th subinterval.

Now, taking into account (17), we can rewrite (15) in the following form:

(18) u r ( t ) = a r ( A 0 , t ) λ r + j = 1 m + 1 a r ( D j , t ) λ j + a r ( F , t ) , t [ θ r 1 , θ r ) , r = 1 , m + 1 ¯ .

Substituting the right-hand side of (18) into conditions (8) and (9), at the corresponding limit values, we obtain the following system of linear algebraic equations with respect to parameters λ r , r = 1 , m + 1 ¯ :

(19) j = 0 m B j λ j + 1 + B m + 1 λ m + 1 + B m + 1 { a m + 1 ( A 0 , θ m + 1 ) λ m + 1 + j = 1 m + 1 a m + 1 ( D j , θ m + 1 ) λ j } + r = 1 m + 1 θ r 1 θ r C ( t ) λ r d t + r = 1 m + 1 θ r 1 θ r C ( t ) { a r ( A 0 , t ) λ r + j = 1 m + 1 a r ( D j , t ) λ j } d t = d B m + 1 a m + 1 ( F , θ m + 1 ) r = 1 m + 1 θ r 1 θ r C ( t ) a r ( F , t ) d t ,

(20) λ i + a i ( A 0 , θ i ) λ i + j = 1 m + 1 a i ( D j , θ i ) λ j λ i + 1 = a i ( F , θ i ) , i = 1 , m ¯ .

We denote the matrix corresponding to the left-hand side of the system of equations (19) and (20) by Q ( θ ) and write the system in the following form:

(21) Q ( θ ) λ = F ( θ ) , λ R n ( m + 1 ) ,

where F ( θ ) = ( d B m + 1 a m + 1 ( F , θ m + 1 ) r = 1 m + 1 θ r 1 θ r C ( t ) a r ( F , t ) d t , a 1 ( F , θ 1 ) , , a m ( F , θ m ) ) ,

Q ( θ ) = ( Q p , k ( θ ) ) , p , k = 1 , m + 1 ¯ , i.e. , Q 1 , 1 ( θ ) = B 0 + B m + 1 a m + 1 ( D 1 , θ m + 1 ) + θ 0 θ 1 C ( t ) d t + θ 0 θ 1 C ( t ) a 1 ( A 0 , t ) d t + r = 1 m + 1 θ r 1 θ r C ( t ) a r ( D 1 , t ) d t , Q 1 , k ( θ ) = B k 1 + B m + 1 a m + 1 ( D k , θ m + 1 ) + θ k 1 θ k C ( t ) d t + θ k 1 θ k C ( t ) a k ( A 0 , t ) d t + r = 1 m + 1 θ r 1 θ r C ( t ) a r ( D k , t ) d t , k = 2 , m ¯ , Q 1 , m + 1 ( θ ) = B m + B m + 1 + B m + 1 a m + 1 ( A 0 , θ m + 1 ) + B m + 1 a m + 1 ( D m + 1 , θ m + 1 ) + θ m θ m + 1 C ( t ) d t + θ m θ m + 1 C ( t ) a m + 1 ( A 0 , t ) d t + r = 1 m + 1 θ r 1 θ r C ( t ) a r ( D m + 1 , t ) d t , Q p , p 1 ( θ ) = I + a p 1 ( A 0 , θ p 1 ) + a p 1 ( D p 1 , θ p 1 ) , p = 2 , m + 1 ¯ , Q p , p ( θ ) = a p 1 ( D p , θ p 1 ) I , p = 2 , m + 1 ¯ , Q p , k ( θ ) = a p 1 ( D k , θ p 1 ) , k p , k p 1 , k = 1 , m + 1 ¯ , p = 2 , m + 1 ¯ .

### 1.2 Numerical implementation of the method

We offer the following numerical implementation of the Dzhumabaev parametrization method for solving linear boundary value problem for essentially loaded differential equations with integro-multipoint condition based on the Runge-Kutta method of fourth-order and Simpson’s method.

1. Suppose we have a partition: 0 = θ 0 < θ 1 < θ 2 < < θ m 1 < θ m < θ m + 1 = T . Divide each r th interval [ θ r 1 , θ r ] , r = 1 , m + 1 ¯ , into N r parts with step size h r = ( θ r θ r 1 ) N r . Assume that on each interval [ θ r 1 , θ r ] , r = 1 , m + 1 ¯ , the variable θ ^ takes its discrete values: θ ^ = θ r 1 , θ ^ = θ r 1 + h r , , θ ^ = θ r 1 + ( N r 1 ) h r , θ ^ = θ r , and denote by { θ r 1 , θ r } , r = 1 , m + 1 ¯ , the set of such points.

2. We find the values of ( n × n ) matrices a r ( A 0 , θ ^ ) , a r ( D i , θ ^ ) , i = 1 , m + 1 ¯ , and n vector a r ( F , θ ^ ) on { θ r 1 , θ r } , r = 1 , m + 1 ¯ .

3. Applying Simpson’s method on the set { θ r 1 , θ r } , r = 1 , m + 1 ¯ , we evaluate the definite integrals

W r ( C ) = θ r 1 θ r C ( t ) d t , W r ( F ) = θ r 1 θ r C ( t ) a r ( F , t ) d t , r = 1 , m + 1 ¯ , W r ( A 0 ) = θ r 1 θ r C ( t ) a r ( A 0 , t ) d t , W r ( D j ) = θ r 1 θ r C ( t ) a r ( D j , t ) d t , j = 1 , m + 1 ¯ .

4. Construct the system of linear algebraic equations in parameters

(22) Q h ˜ ( θ ) λ = F h ˜ ( θ ) , λ R n ( m + 1 ) ,

and find its solution λ h ˜ . As noted earlier, the elements of λ h ˜ = ( λ 1 h ˜ , λ 2 h ˜ , , λ m h ˜ , λ m + 1 h ˜ ) are the values of an approximate solution to problem (1), (2) at the left end-points of the subintervals: x h ˜ r ( θ r 1 ) = λ r h ˜ , r = 1 , m + 1 ¯ .

5. To define the values of an approximate solution at the remaining points of set { θ r 1 , θ r } , r = 1 , m + 1 ¯ , we solve the Cauchy problems:

d x d t = A 0 ( t ) x + j = 1 m + 1 D j ( t ) λ j h ˜ + F ( t ) , x ( θ r 1 ) = λ r h ˜ , t [ θ r 1 , θ r ] , r = 1 , m + 1 ¯ .

We found the solutions to Cauchy problems by using the fourth-order Runge-Kutta method. Thus, the algorithm allows us to find the numerical solution to the problem (1), (2), when the matrices G ( θ ) , Q ( θ ) are invertible.

### 1.3 Solvability of the problem

In this section, we establish conditions for the unique solvability of problem (1), (2). To prove the main assertion, we need the following lemma.

### Lemma 1

If G ( θ ) is invertible, then the following assertions hold:

1. If λ ˜ R n ( m + 1 ) and μ ˜ R n ( m + 1 ) are solutions of systems (21) and (12), respectively, and the function system u ˜ [ t ] C ( [ 0 , T ] , θ , R n ( m + 1 ) ) is a solution to the Cauchy problem (6), (7) with λ r = λ ˜ r , μ r = μ ˜ r , r = 1 , m + 1 ¯ , then the function x ˜ ( t ) , defined by the equalities:

x ˜ ( t ) = u ˜ r ( t ) + λ ˜ r , t [ θ r 1 , θ r ) , lim t θ r 1 + 0 x ˜ ˙ ( t ) = μ ˜ r , r = 1 , m + 1 ¯ , x ˜ ( T ) = λ ˜ m + 1 + lim t T 0 u ˜ m + 1 ( t ) ,

is a solution to problem (1), (2);

2. The vectors λ R n ( m + 1 ) and μ R n ( m + 1 ) , composed of the values of the solution x ( t ) and its derivative x ˙ ( t ) to problem (1), (2) at the partition points λ r = lim t θ r 1 + 0 x ( t ) , and μ r = lim t θ r 1 + 0 x ˙ ( t ) , r = 1 , m + 1 ¯ , satisfy the systems (21) and (12), respectively.

The evidence is similar in concept to the proof of Lemma 1 with slight variations [27].

We present the main theorem on the existence of a unique solution to problem (1), (2) in terms of matrices G ( θ ) and Q ( θ ) .

### Theorem 1

Let the matrices G ( θ ) : R n ( m + 1 ) R n ( m + 1 ) and Q ( θ ) : R n ( m + 1 ) R n ( m + 1 ) be invertible. Then boundary value problem (1), (2) has a unique solution x ( t ) for any f ( t ) C ( [ 0 , T ] , R n ) , d R n .

### Proof

Let f ( t ) C ( [ 0 , T ] , R n ) and d R n . We find an unique solution to the systems (12) and (21), using the invertibility of matrices G ( θ ) and Q ( θ ) : μ = [ G ( θ ) ] 1 P ( θ , λ ) and λ = [ Q ( θ ) ] 1 F ( θ ) . Solving the Cauchy problem (6), (7) with λ = λ and μ = μ , define the function system u [ t ] . The invertibility of the matrices G ( θ ) and Q ( θ ) lead to the existence of unique function system u [ t ] with the elements u r [ t ] , defined by the right-hand side of (15) at λ = λ . Then, according to Lemma 1, the function x ( t ) , defined by the equalities: x ( t ) = u r ( t ) + λ r , t [ θ r 1 , θ r ) , lim t θ r 1 + 0 x ˙ ( t ) = μ r , r = 1 , m + 1 ¯ , x ( T ) = λ m + 1 + lim t T 0 u m + 1 ( t ) , is a solution to problem (1), (2). Uniqueness of the solution is proved by contradiction. Theorem 1 is proved.□

To illustrate the fulfillment of the theorem’s conditions and the proposed approach of the numerical solving of the boundary value problem for systems of essentially loaded differential equations with integro-multipoint condition (1), (2) based on the Dzhumabaev parametrization method, let us consider the following examples.

### Example 1

We consider a linear boundary value problem for essentially loaded differential equations with an integro-multipoint condition:

(23) d x d t = A 0 ( t ) x + A 1 ( t ) x ˙ ( θ 1 ) + A 2 ( t ) x ˙ ( θ 2 ) + A 3 ( t ) x ˙ ( θ 3 ) + f ( t ) , t ( 0 , 1 ) ,

(24) j = 0 4 B j x ( θ j ) + θ 0 θ 4 C ( t ) x ( t ) d t = d , d R 3 , x R 3 .

Here,

θ 0 = 0 , θ 1 = 1 4 , θ 2 = 1 2 , θ 3 = 3 4 , θ 4 = T = 1 , f ( t ) = 36 t 108 t 3 132 t 2 48 t 4 396 48 t 2 90 t 3 12 t 4 252 t 696 360 t 4 228 t 3 108 t 2 408 t 54 ,

A 0 ( t ) = t 2 t 2 t 1 t 3 0 3 t 2 5 6 t , A 1 ( t ) = 1 t t 2 4 t 6 0 3 t 4 t 2 , A 2 ( t ) = t 0 6 t 2 4 5 0 t 5 7 t , A 3 ( t ) = t 2 3 t 0 t + 2 6 7 0 t 3 1 , B 0 = 1 4 5 6 3 2 8 0 1 , B 1 = 3 0 5 3 2 1 5 3 5 , B 2 = 1 6 8 3 0 7 6 4 2 , B 3 = 5 0 6 2 1 6 4 6 8 , B 4 = 2 6 7 6 0 5 5 11 3 , C ( t ) = t 8 0 4 2 1 0 3 t 5 , d = 368 869 454 .

We use the numerical implementation of algorithm. The accuracy of the solution depends on the accuracy of solving the Cauchy problem on subintervals and evaluating definite integrals. We provide the results of the numerical implementation of algorithm by partitioning the subintervals [ 0 , 0.25 ] , [ 0.25 , 0.5 ] , [ 0.5 , 0.75 ] , [ 0.75 , 1 ] with step h = 0.025 . Solving the system of equations (22), we obtain the numerical values of the parameters

λ 1 h ˜ = 12.00002405 6.00001267 0.0000129 , λ 2 h ˜ = 12.00002727 3.00001292 6.00000624 , λ 3 h ˜ = 24.00002692 0.00001382 18.00000584 , λ 4 h ˜ = 48.00001942 2.99998431 36.00000785 .

Exact solution of problem (23) and (24) is the following: x ( t ) = 96 t 2 24 t + 12 12 t 6 48 t 2 + 12 t .

The differences between the exact and approximate solutions to problem (23) and (24)  ε i = x ( i ) ( t k ) x ˜ ( i ) ( t k ) , i = 1 , 3 ¯ , are provided in Table 1. Table 2 provides the difference between numerical and exact solutions of problem (23) and (24), where we solve Cauchy problems by using the Runge-Kutta Fehlberg method.

Table 1

Error analysis in Example 1 (using the fourth-order Runge-Kutta method)

k t k ε 1 ε 2 ε 3 k t k ε 1 ε 2 ε 3
0 0 2.4 × 1 0 5 1.27 × 1 0 5 1.29 × 1 0 5 20 0.5 2.69 × 1 0 5 1.38 × 1 0 5 0.58 × 1 0 5
1 0.025 2.45 × 1 0 5 1.27 × 1 0 5 1.19 × 1 0 5 21 0.525 2.66 × 1 0 5 1.39 × 1 0 5 0.61 × 1 0 5
2 0.05 2.49 × 1 0 5 1.27 × 1 0 5 1.09 × 1 0 5 22 0.55 2.62 × 1 0 5 1.41 × 1 0 5 0.63 × 1 0 5
3 0.075 2.52 × 1 0 5 1.27 × 1 0 5 1.01 × 1 0 5 23 0.575 2.57 × 1 0 5 1.42 × 1 0 5 0.66 × 1 0 5
4 0.1 2.56 × 1 0 5 1.27 × 1 0 5 0.93 × 1 0 5 24 0.6 2.51 × 1 0 5 1.44 × 1 0 5 0.69 × 1 0 5
5 0.125 2.59 × 1 0 5 1.27 × 1 0 5 0.86 × 1 0 5 25 0.625 2.44 × 1 0 5 1.46 × 1 0 5 0.72 × 1 0 5
6 0.15 2.63 × 1 0 5 1.27 × 1 0 5 0.8 × 1 0 5 26 0.65 2.36 × 1 0 5 1.48 × 1 0 5 0.74 × 1 0 5
7 0.175 2.66 × 1 0 5 1.28 × 1 0 5 0.75 × 1 0 5 27 0.675 2.28 × 1 0 5 1.5 × 1 0 5 0.77 × 1 0 5
8 0.2 2.68 × 1 0 5 1.28 × 1 0 5 0.7 × 1 0 5 28 0.7 2.18 × 1 0 5 1.52 × 1 0 5 0.78 × 1 0 5
9 0.225 2.71 × 1 0 5 1.29 × 1 0 5 0.66 × 1 0 5 29 0.725 2.07 × 1 0 5 1.54 × 1 0 5 0.79 × 1 0 5
10 0.25 2.73 × 1 0 5 1.29 × 1 0 5 0.62 × 1 0 5 30 0.75 1.94 × 1 0 5 1.57 × 1 0 5 0.79 × 1 0 5
11 0.275 2.74 × 1 0 5 0.13 × 1 0 5 0.59 × 1 0 5 31 0.775 1.8 × 1 0 5 1.6 × 1 0 5 0.76 × 1 0 5
12 0.3 2.76 × 1 0 5 1.31 × 1 0 5 0.57 × 1 0 5 32 0.8 1.65 × 1 0 5 1.63 × 1 0 5 0.71 × 1 0 5
13 0.325 2.77 × 1 0 5 1.31 × 1 0 5 0.55 × 1 0 5 33 0.825 1.47 × 1 0 5 1.67 × 1 0 5 0.63 × 1 0 5
14 0.35 2.77 × 1 0 5 1.32 × 1 0 5 0.54 × 1 0 5 34 0.85 1.28 × 1 0 5 1.71 × 1 0 5 0.51 × 1 0 5
15 0.375 2.77 × 1 0 5 1.33 × 1 0 5 0.54 × 1 0 5 35 0.875 1.06 × 1 0 5 1.75 × 1 0 5 0.34 × 1 0 5
16 0.4 2.77 × 1 0 5 1.34 × 1 0 5 0.54 × 1 0 5 36 0.9 0.82 × 1 0 5 1.8 × 1 0 5 0.1 × 1 0 5
17 0.425 2.76 × 1 0 5 1.35 × 1 0 5 0.54 × 1 0 5 37 0.925 0.54 × 1 0 5 1.85 × 1 0 5 0.22 × 1 0 5
18 0.45 2.74 × 1 0 5 1.36 × 1 0 5 0.55 × 1 0 5 38 0.95 0.23 × 1 0 5 1.91 × 1 0 5 0.66 × 1 0 5
19 0.475 2.72 × 1 0 5 1.37 × 1 0 5 0.57 × 1 0 5 39 0.975 0.12 × 1 0 5 1.98 × 1 0 5 1.23 × 1 0 5
20 0.5 2.69 × 1 0 5 1.38 × 1 0 5 0.58 × 1 0 5 40 1 0.53 × 1 0 5 2.06 × 1 0 5 1.97 × 1 0 5
Table 2

Error analysis in Example 1 (using the Runge-Kutta Fehlberg method)

k t k ε 1 ε 2 ε 3 k