A Collocation Method Based on the Bernoulli Operational Matrix for Solving High-Order Linear Complex Differential Equations in a Rectangular Domain

Abstract

This paper contributes a new matrix method for the solution of high-order linear complex differential equations with variable coefficients in rectangular domains under the considered initial conditions. On the basis of the presented approach, the matrix forms of the Bernoulli polynomials and their derivatives are constructed, and then by substituting the collocation points into the matrix forms, the fundamental matrix equation is formed. This matrix equation corresponds to a system of linear algebraic equations. By solving this system, the unknown Bernoulli coefficients are determined and thus the approximate solutions are obtained. Also, an error analysis based on the use of the Bernoulli polynomials is provided under several mild conditions. To illustrate the efficiency of our method, some numerical examples are given.

1. Introduction

Complex differential equations have a great popularity in science and engineering. In real world, physical events can be modeled by complex differential equations usually. For instance, the vibrations of a one-mass system with two DOFs are mostly described using differential equations with a complex dependent variable [1, 2]. The various applications of differential equations with complex dependent variables are introduced in [2]. Since a huge size of such equations cannot be solved explicitly, it is often necessary to resort to approximation and numerical techniques.

In recent years, the studies on complex differential equations, such as a geometric approach based on meromorphic function in arbitrary domains [3], a topological description of solutions of some complex differential equations with multivalued coefficients [4], the zero distribution [5], growth estimates [6] of linear complex differential equations, and also the rational together with the polynomial approximations of analytic functions in the complex plane [7, 8], were developed very rapidly and intensively.

Since the beginning of 1994, the Laguerre, Chebyshev, Taylor, Legendre, Hermite, and Bessel (matrix and collocation) methods have been used in the works in [9–19] to solve linear differential, integral, and integrodifferential-difference equations and their systems. Also, the Bernoulli matrix method has been used to find the approximate solutions of differential and integrodifferential equations [20–22].

In this paper, in the light of the above-mentioned methods and by means of the matrix relations between the Bernoulli polynomials and their derivatives, we develop a new method called the Bernoulli collocation method (BCM) for solving high-order linear complex differential equation 𝑓^(𝑚)(𝑧)+^𝑚−1_𝑘=0𝑃_𝑘(𝑧)𝑓^(𝑘)(𝑧)(𝑧),[][],=𝐺𝑚≥1,𝑧=𝑥+𝑖𝑦,𝑥∈𝑎,𝑏,𝑦∈𝑐,𝑑 (1) with the initial conditions 𝑓^(𝑟)(0)=𝛽_𝑟,𝑟=0,1,…,𝑚−1. (2) Evidently, we will let 𝑧 denote a variable point in the complex plane. Its real and imaginary parts will be denoted by 𝑥 and 𝑦, respectively.

Here, the coefficients 𝑃_𝑘(𝑧) and the known function 𝐺(𝑧) together with the unknown function 𝑓(𝑧) are holomorphic (or analytic) functions in the domain 𝐷={𝑧∈ℂ∶𝑧=𝑥+𝑖𝑦,  𝑎≤𝑥≤𝑏,  𝑐≤𝑦≤𝑑∶𝑎,𝑏,𝑐,𝑑∈ℝ} where the coefficients 𝛽_𝑟 are appropriate complex constants.

We assume that the solution of (1) under the conditions (2) is approximated in the form 𝑓(𝑧)≈𝑓_𝑁(𝑧)=^𝑁_𝑛=0𝑓_𝑛𝐵_𝑛(𝑧),𝑧∈𝐷, (3) which is the truncated Bernoulli series of the unknown function 𝑓(𝑧), where all of 𝑓_𝑛(𝑛=0,1,…,𝑁) are the Bernoulli coefficients to be determined. We also use the collocation points 𝑧_𝑝𝑝=𝑥_𝑝+𝑖𝑦_𝑝,𝑥_𝑝=𝑎+𝑏−𝑎𝑁𝑦𝑝,_𝑝=𝑐+𝑑−𝑐𝑁𝑝,𝑝=0,1,…,𝑁. (4)

In this paper, by generalizing the methods [20, 21] from real calculus to the complex calculus, we propose a new matrix method which is based on the Bernoulli operational matrix of differentiation and a uniform collocation scheme. It should be noted that since an ordinary complex differential equation equals to a system of partial differential equations (see Section 4) the methods that based on high-order Gauss quadrature rules [23, 24] could not be effective. Needing to more CPU time from one side and ill conditioning of the associated algebraic problem from another side are two disadvantages of such methods. Therefore, implementing an easy to use approach such as the methods that based on operational matrices is necessary for solving any practical problem.

The rest of this paper is organized as follows. In Section 2, we review some notations from complex calculus and also provide several properties of the Bernoulli polynomials. Section 3 is devoted to the proposed matrix method. Error analysis and accuracy of the approximated solution by the aid of the Bernoulli polynomials is given in Section 4. Several illustrative examples are provided in Section 5 for confirming the effectiveness of the presented method. Section 6 contains some conclusions and notations about the future works.

2. Review on Complex Calculus and the Bernoulli Polynomials

This Section is divided into two subsections. In the first subsection we review some notations from complex calculus specially the concept of differentiability in the complex plane under some remarks. Then we recall several properties of the Bernoulli polynomials and introduce the operational matrix of differentiation of the Bernoulli polynomials in the complex form.

2.1. Review on Complex Calculus

From the definition of derivative in the complex form, it is immediate that a constant function is differentiable everywhere, with derivative 0, and that the identity function (the function 𝑓(𝑧)=𝑧) is differentiable everywhere, with derivative 1. Just as in elementary calculus one can show from the last statement, by repeated applications of the product rule, that, for any positive integer 𝑛, the function 𝑓(𝑧)=𝑧^𝑛 is differentiable everywhere, with derivative 𝑛𝑧^𝑛−1. This, in conjunction with the sum and product rules, implies that every polynomial is everywhere differentiable: if 𝑓(𝑧)=𝑐_𝑛𝑧^𝑛+𝑐_𝑛−1𝑧^𝑛−1+⋯+𝑐₁𝑧+𝑐₀, where 𝑐₀,…,𝑐_𝑛 are complex constants, then 𝑓^′(𝑧)=𝑛𝑐_𝑛𝑧^𝑛−1+(𝑛−1)𝑐_𝑛−1𝑧^𝑛−2+⋯+𝑐₁.

Remark 1. The function 𝑓=𝑢+𝑖𝑣 is differentiable (in the complex sense) at 𝑧₀ if and only if 𝑢 and 𝑣 are differentiable (in the real sense) at 𝑧₀ and their first partial derivatives satisfy the relations 𝜕𝑢(𝑧₀)/𝜕𝑥=𝜕𝑣(𝑧₀)/𝜕𝑦,𝜕𝑢(𝑧₀)/𝜕𝑦=−𝜕𝑣(𝑧₀)/𝜕𝑥. In that case 𝑓^′𝑧₀=𝜕𝑢𝑧₀𝑧𝜕𝑥+𝑖𝜕𝑣₀=𝑧𝜕𝑥𝜕𝑣₀𝑧𝜕𝑦−𝑖𝜕𝑢₀.𝜕𝑦 (5)

Remark 2. Two partial differential equations 𝜕𝑢=𝜕𝑥𝜕𝑣,𝜕𝑦𝜕𝑢𝜕𝑦=−𝜕𝑣𝜕𝑥 (6) are called the Cauchy-Riemann equations for the pair of functions 𝑢,𝑣. As seen above (i.e., the Remark 1), the equations are satisfied by the real and imaginary parts of a complex-valued function at each point where that function is differentiable.

Remark 3 (sufficient condition for complex differentiability). Let the complex-valued function 𝑓=𝑢+𝑖𝑣 be defined in the open subset 𝐺 of the complex plane ℂ and assume that 𝑢 and 𝑣 have first partial derivatives in 𝐺. Then 𝑓 is differentiable at each point where those partial derivatives are continuous and satisfy the Cauchy-Riemann equations.

Definition 4. A complex-valued function that is defined in an open subset 𝐺 of the complex plane ℂ and differentiable at every point of 𝐺 is said to be holomorphic (or analytic) in 𝐺. The simplest examples are polynomials, which are holomorphic in ℂ, and rational functions, which are holomorphic in the regions where they are defined. Moreover, the elementary functions such as exponential function, the logarithm function, trigonometric and inverse trigonometric functions, and power functions all have complex versions that are holomorphic functions. It should be noted that if the real and imaginary parts of a complex-valued function have continuous first partial derivatives obeying the Cauchy-Riemann equations, then the function is holomorphic.

Remark 5 (complex partial differential operators). The partial differential operators 𝜕/𝜕𝑥 and 𝜕/𝜕𝑦 are applied to a complex-valued function 𝑓=𝑢+𝑖𝑣 in the natural way: 𝜕𝑓=𝜕𝑥𝜕𝑢𝜕𝑥+𝑖𝜕𝑣,𝜕𝑥𝜕𝑓=𝜕𝑦𝜕𝑢𝜕𝑦+𝑖𝜕𝑣.𝜕𝑦 (7) We define the complex partial differential operators 𝜕/𝜕𝑧 and 𝜕/𝜕𝑧 by 𝜕=1𝜕𝑧2𝜕𝜕𝜕𝑥−𝑖,𝜕𝜕𝑦𝜕𝑧=12𝜕𝜕𝜕𝑥+𝑖.𝜕𝑦 (8) Thus, 𝜕/𝜕𝑥=𝜕/𝜕𝑧+𝜕/𝜕𝑧and𝜕/𝜕𝑦=𝑖(𝜕/𝜕𝑧−𝜕/𝜕𝑧).
Intuitively one can think of a holomorphic function as a complex-valued function in an open subset of ℂ that depends only on 𝑧, that is, independent of 𝑧. We can make this notion precisely as follows. Suppose the function 𝑓=𝑢+𝑖𝑣 is defined and differentiable in an open set. One then has 𝜕𝑓=1𝜕𝑧2𝜕𝑢+𝜕𝑥𝜕𝑣+𝑖𝜕𝑦2𝜕𝑣−𝜕𝑥𝜕𝑢,𝜕𝑦𝜕𝑓𝜕𝑧=12𝜕𝑢−𝜕𝑥𝜕𝑣+𝑖𝜕𝑦2𝜕𝑣+𝜕𝑥𝜕𝑢.𝜕𝑦 (9)

The Cauchy-Riemann equations thus can be written 𝜕𝑓/𝜕𝑧=0. As this is the condition for 𝑓 to be holomorphic, it provides a precise meaning for the statement: a holomorphic function is one that is independent of 𝑧. If 𝑓 is holomorphic, then (not surprisingly) 𝑓^′=𝜕𝑓/𝜕𝑧, as the following calculation shows: 𝑓^′=𝜕𝑓=𝜕𝑥𝜕𝑓+𝜕𝑧𝜕𝑓𝜕𝑧=𝜕𝑓.𝜕𝑧 (10)

2.2. The Bernoulli Polynomials and Their Operational Matrix

The Bernoulli polynomials play an important role in different areas of mathematics, including number theory and the theory of finite differences. The classical Bernoulli polynomials 𝐵_𝑛(𝑥) are usually defined by means of the exponential generating functions (see [21]) 𝑤𝑒^𝑥𝑤𝑒^𝑤=−1^∞_𝑘=0𝐵_𝑘(𝑥)𝑤^𝑘,𝑤|).𝑘!(|<2𝜋 (11) The following familiar expansion (see [20]): ^𝑛𝑘=0𝑘𝐵𝑛+1_𝑘(𝑥)=(𝑛+1)𝑥^𝑛 (12) is the most primary property of the Bernoulli polynomials. The first few Bernoulli polynomials are 𝐵₀(𝑥)𝐵=1,₁(𝑥)1=𝑥−2,𝐵₂(𝑥)=𝑥²1−𝑥+6,𝐵₃(𝑥)=𝑥³−32𝑥²+12𝐵𝑥,₄(𝑥)=𝑥⁴−2𝑥³+𝑥²−1.30 (13) The Bernoulli polynomials satisfy the well-known relations (see [21]) 𝑑𝐵_𝑛(𝑥)𝑑𝑥=𝑛𝐵_𝑛−1(𝑥),(𝑛≥1₁₀),𝐵_𝑛(𝑥)𝑑𝑥=0,(𝑛≥1). (14)

The Bernoulli polynomials have another specific property that is satisfied in the following linear homogeneous recurrence relation [20]: 𝐵_𝑛(𝑥)=1𝑥−2𝐵_𝑛−1(𝑥)−1𝑛^𝑛−2_𝑘=0𝑛𝑘𝐵_𝑛−𝑘(0)𝐵_𝑘(𝑥). (15) Also, the Bernoulli polynomials satisfy the following interesting property [25]: ¹₀𝐵_𝑚(𝑥)𝐵_𝑛(𝑥)=()𝑑𝑥−1^𝑚+𝑛𝑚!𝑛!()!𝐵𝑚+𝑛^𝑚+𝑛(0). (16) Moreover, 𝐵_𝑛(𝑥) satisfy the differential equation [20] 𝐵_𝑛(0)𝑦𝑛!^(𝑛)(𝑥)+𝐵_𝑛−1(0)()!𝑦𝑛−1^(𝑛−1)(𝑥)𝐵+⋯+²(0)𝑦2!^′′(𝑥)+12−𝑥𝑦^′(𝑥)+𝑛𝑦(𝑥)=0. (17) According to the discussions in [20], the Bernoulli polynomials form a complete basis over the interval [0,1].

If we introduce the Bernoulli vector 𝐵(𝑥) in the form 𝐵(𝑥)=[𝐵₀(𝑥),𝐵₁(𝑥),…,𝐵_𝑁(𝑥)], then the derivative of 𝐵(𝑥), with the aid of the first property of (14), can be expressed in the matrix form by 𝐵^′₀(𝑥)𝐵^′₁(𝑥)𝐵^′₂(𝑥)⋮𝐵^′_𝑁−1(𝑥)𝐵^′_𝑁(𝑥)𝐵^′(𝑥)^𝑇=000⋯000100⋯000020⋯000⋮⋮⋮⋱⋮⋮⋮000…𝑁−100000…0𝑁0𝑀𝐵₀(𝑥)𝐵₁(𝑥)𝐵₂(𝑥)⋮𝐵_𝑁−1(𝑥)𝐵_𝑁(𝑥)𝐵^𝑇(𝑥), (18) where 𝑀 is the (𝑁+1)×(𝑁+1) operational matrix of differentiation. Note that if we replace the real variable 𝑥 by the complex variable 𝑧 in the above relation, again we reach to the same result, since (𝑧^𝑛)^′=𝑛𝑧^𝑛−1.

Accordingly, the 𝑘th derivative of 𝐵(𝑥) can be given by 𝐵^′(𝑥)^𝑇=𝑀𝐵(𝑥)^𝑇⟹𝐵⁽¹⁾(𝑥)=𝐵(𝑥)𝑀^𝑇,𝐵⁽²⁾(𝑥)=𝐵⁽¹⁾(𝑥)𝑀^𝑇=𝐵(𝑥)𝑀^𝑇2,𝐵⁽³⁾(𝑥)=𝐵⁽¹⁾(𝑥)𝑀^𝑇2=𝐵(𝑥)𝑀^𝑇3,⋮𝐵^(𝑘)(𝑥)=𝐵(𝑥)𝑀^𝑇𝑘, (19) where 𝑀 is defined in (18).

We recall that the Bernoulli expansions and Taylor series are not based on orthogonal functions; nevertheless, they possess the operational matrices of differentiations (and also integration). However, since the integration of the cross product of two Taylor series vectors is given in terms of a Hilbert matrix, which are known to be ill conditioned, the applications of Taylor series in the integration form of view are limited. But in differential form of view, (see for instance [10–12, 16, 17] and the references therein) use the operational matrix of derivatives such as Taylor in a huge number of research works. For approximating an arbitrary unknown function, the advantages of the Bernoulli polynomials over orthogonal polynomials as shifted Legendre polynomials are the following.(i)The operational matrix of differentiation in the Bernoulli polynomials has less nonzero elements than for shifted Legendre polynomials. Because for Bernoulli polynomials, the nonzero elements of this matrix are located in below (or above) its diagonal. However for the shifted Legendre polynomials is a strictly lower (or upper) filled triangular matrix.(ii)The Bernoulli polynomials have less terms than shifted Legendre polynomials. For example, 𝐵₆(𝑥) (the 6th Bernoulli polynomial) has 5 terms while 𝑃₆(𝑥) (the 6th shifted Legendre polynomial) has 7 terms, and this difference will increase by increasing the index. Hence, for approximating an arbitrary function we use less CPU time by applying the Bernoulli polynomials as compared to shifted Legendre polynomials; this issue is claimed in [25] and is proved in its examples for solving nonlinear optimal control problems.(iii)The coefficient of individual terms in the Bernoulli polynomials is smaller than the coefficient of individual terms in the shifted Legendre polynomials. Since the computational errors in the product are related to the coefficients of individual terms, the computational errors are less by using the Bernoulli polynomials.

3. Basic Idea

In this section by applying the Bernoulli operational matrix of differentiation and also the collocation scheme, the basic idea of this paper would be constructed. We again consider (1) and its approximated solution 𝑓_𝑁(𝑧) in the form (3). Trivially, 𝑓_𝑁(𝑧) could be rewritten in the vector form 𝑓_𝑁(𝑧)=𝐵(𝑧)𝐹,𝐵(𝑧)=𝐵₀(𝑧)𝐵₁(𝑧)⋯𝐵_𝑁(𝑧),𝐹=𝑓₀𝑓₁⋯𝑓_𝑁^𝑇. (20) By using (19) in the complex form one can conclude that 𝑓_𝑁^(𝑘)(𝑧)=𝐵(𝑧)𝑀^𝑇𝑘𝐹,𝑘≤𝑁, (21) where 𝑀 is introduced in (19).

For the collocation points 𝑧=𝑧_𝑝𝑝(𝑝=0,1,…,𝑁), the matrix relation (21) becomes 𝑓_𝑁^(𝑘)𝑧_𝑝𝑝=𝐵𝑧_𝑝𝑝𝑀^𝑇𝑘𝐹,𝑝=0,1,…,𝑁, (22) where 𝐵(𝑧_𝑝𝑝)=[𝐵₀(𝑧_𝑝𝑝)𝐵₁(𝑧_𝑝𝑝)⋯𝐵_𝑁(𝑧_𝑝𝑝)]. For more details, one can restate (22) as follows: 𝑓_𝑁^(𝑘)𝑧₀₀=𝐵𝑧₀₀𝑀^𝑇𝑘𝑓𝐹,_𝑁(𝑘)𝑧₁₁=𝐵𝑧₁₁𝑀^𝑇𝑘⋮𝑓𝐹,_𝑁(𝑘)𝑧_𝑁𝑁=𝐵𝑧_𝑁𝑁𝑀^𝑇𝑘𝐹. (23) The matrix vector of the above-mentioned equations is 𝐹^(𝑘)=𝑓_𝑁^(𝑘)𝑧₀₀𝑓_𝑁^(𝑘)𝑧₁₁⋮𝑓_𝑁^(𝑘)𝑧_𝑁𝑁=𝐿𝑀^𝑇𝑘𝐹, (24) where 𝐵𝑧𝐿=⁰⁰𝐵𝑧₁₁⋮𝐵𝑧_𝑁𝑁=𝐵₀𝑧₀₀𝐵₁𝑧₀₀𝐵₂𝑧₀₀⋯𝐵_𝑁𝑧₀₀𝐵₀𝑧₁₁𝐵₁𝑧₁₁𝐵₂𝑧₁₁⋯𝐵_𝑁𝑧₁₁𝐵⋮⋮⋮⋱⋮₀𝑧_𝑁𝑁𝐵₁𝑧_𝑁𝑁𝐵₂𝑧_𝑁𝑁⋯𝐵_𝑁𝑧_𝑁𝑁. (25) On the other hand by substituting the collocation points 𝑧=𝑧_𝑝𝑝 defined by (4) into (1), we have 𝑓_𝑁^(𝑚)𝑧_𝑝𝑝+^𝑚−1_𝑘=0𝑃_𝑘𝑧_𝑝𝑝𝑓_𝑁^(𝑘)𝑧_𝑝𝑝=𝐺𝑧_𝑝𝑝,𝑝=0,1,…,𝑁. (26) The associated matrix-vector form of the above equations has the following form by the aid of (24): 𝐹^(𝑚)+^𝑚−1_𝑘=0𝑃_𝑘𝐹^(𝑘)=𝐿𝑀^𝑇𝑚𝐹+^𝑚−1_𝑘=0𝑃_𝑘𝐿𝑀^𝑇𝑘𝐹=𝐺, (27) where 𝐺=𝐺𝑧₀₀𝐺𝑧₁₁⋯𝐺𝑧_𝑁𝑁^𝑇,𝑃_𝑘=𝑃_𝑘𝑧₀₀0⋯00𝑃_𝑘𝑧₁₁⋯0⋮⋮⋱⋮00⋯𝑃_𝑘𝑧_𝑁𝑁. (28) Since the vector 𝐹 is unknown and should be determined, therefore, the matrix-vector equation (27) could be rewritten in the following form: 𝑊𝐹=𝐺,or[𝑊;𝐺]=𝑤_𝑝𝑞;𝑔_𝑝,𝑝,𝑞=0,1,…,𝑁, (29) where 𝑊=𝐿(𝑀^𝑇)^𝑚+∑^𝑚−1_𝑘=0𝑃_𝑘𝐿(𝑀^𝑇)^𝑘.

We now write the vector form of the initial conditions (2) by the aid of (21) as follows: 𝑓^(𝑟)(0)=𝐵(0)𝑀^𝑇𝑟𝐹=𝛽_𝑟,𝑟=0,1,…,𝑚−1. (30) In other words the vector form of the initial conditions could be rewritten as 𝑈_𝑟𝐹=𝛽_𝑟 where 𝑈_𝑟=𝐵(0)(𝑀^𝑇)^𝑟  (𝑟=0,1,…,𝑚−1). Trivially the augmented form of these equations is 𝑈_𝑟;𝛽_𝑟=𝑢_𝑟0𝑢_𝑟1⋯𝑢_𝑟𝑁;𝛽_𝑟,𝑟=0,1,…,𝑚−1. (31) Consequently, to find the unknown Bernoulli coefficients 𝑓_𝑛,𝑛=0,1,…,𝑁, related with the approximate solution of the problem (1) under the initial conditions (2), we need to replace the 𝑚 rows of (31) by the last 𝑚 rows of the augmented matrix (29) and hence we have new augmented matrix 𝐺=𝑤𝑊;⁰⁰𝑤₀₁⋯𝑤_0𝑁;𝑔₀𝑤₁₀𝑤₁₁⋯𝑤_1𝑁;𝑔₁𝑤⋮⋮⋱⋮⋮⋮_{𝑁−𝑚,0}𝑤_{𝑁−𝑚,1}⋯𝑤_{𝑁−𝑚,𝑁};𝑔_𝑁−𝑚𝑢₀₀𝑢₀₁⋯𝑢_0𝑁;𝛽₀𝑢₁₀𝑢₁₁⋯𝑢_1𝑁;𝛽₁𝑢⋮⋮⋱⋮⋮⋮_𝑚−1,0𝑢_𝑚−1,1⋯𝑢_{𝑚−1,𝑁};𝛽_𝑚−1, (32) or the corresponding matrix-vector equation 𝑊𝐹=𝐺. (33) If det(𝑊)≠0 one can rewrite (32) in the form 𝐹=(𝑊)⁻¹𝐺 and the vector 𝐹 is uniquely determined. Thus the 𝑚th order linear complex differential equation with variable coefficients (1) under the conditions (2) has an approximated solution. This solution is given by the truncated Bernoulli series (3). Also we can easily check the accuracy of the obtained solutions as follows [12, 26]. Since the truncated Bernoulli series (3) is an approximate solution of (1), when the solutions 𝑓(𝑧) and its derivatives are substituted in (1), the resulting equation must be satisfied approximately; that is, for 𝑧=𝑧_𝑖∈𝐷,𝑖=0,1,2,…𝐸𝑧_𝑖=𝑓^(𝑚)𝑧_𝑖+^𝑚−1_𝑘=0𝑃_𝑘𝑧_𝑖𝑓^(𝑘)𝑧_𝑖−𝐺𝑧_𝑖𝑧≅0,or𝐸_𝑖≤10^−𝑘𝑖𝑘_𝑖isanypositiveinteger. (34) If max(10^−𝑘𝑖)=10^−𝑘 (𝑘 is any positive integer) is prescribed, then the truncation limit 𝑁 may be increased until the values 𝐸(𝑧_𝑖) at each of the points 𝑧_𝑖 become smaller than the prescribed 10^−𝑘 (for more details see [10–12, 16, 17]).

4. Error Analysis and Accuracy of the Solution

This section is devoted to provide an error bound for the approximated solution which may be obtained by the Bernoulli polynomials. We emphasized that this section is given for showing the efficiency of the Bernoulli polynomial approximation and is independent of the proposed method which is provided for showing how a complex ordinary differential equation (ODE) is equivalent to a system of partial differential equations (PDEs). After conveying this subject, we transform the obtained system of PDEs (together with the initial conditions (2)) to a system of two dimensional Volterra integral equations in a special case. Before presenting the main Theorem of this section, we need to recall some useful corollaries and lemmas. Therefore, the main theorem could be stated which guarantees the convergence of the truncated Bernoulli series to the exact solution under several mild conditions.

Now suppose that 𝐻=𝐿²[0,1] and {𝐵₀(𝑥),𝐵₁(𝑥),…,𝐵_𝑁(𝑥)}⊂𝐻 be the set of the Bernoulli polynomials and 𝑌=span𝐵₀(𝑥),𝐵₁(𝑥),…,𝐵_𝑁(𝑥) (35) and 𝑔 be an arbitrary element in 𝐻. Since 𝑌 is a finite dimensional vector space, 𝑔 has the unique best approximation 𝑔∈𝑌, such that ∀𝑦∈𝑌,𝑔−𝑔_∞≤𝑔−𝑦_∞. (36) Since 𝑔∈𝑌, there exists the unique coefficients 𝑔₀,𝑔₁,…,𝑔_𝑁 such that 𝑔≈𝑔=^𝑁_𝑛=0𝑔_𝑛𝐵_𝑛(𝑥)=𝐵(𝑥)𝐺^𝑇,𝐵(𝑥)=𝐵₀(𝑥),𝐵₁(𝑥),…,𝐵_𝑁(𝑥,𝑔)𝐺=₀,𝑔₁,…,𝑔_𝑁. (37)

Corollary 6. Assume that 𝑔∈𝐻=𝐿²[0,1] be an enough smooth function and also is approximated by the Bernoulli serie ∑^∞_𝑛=0𝑔_𝑛𝐵_𝑛(𝑥), then the coefficients 𝑔_𝑛 for all 𝑛=0,1,…,∞ can be calculated from the following relation: 𝑔_𝑛=1𝑛!¹₀𝑔^(𝑛)(𝑥)𝑑𝑥. (38)

Proof. See [21].

In practice one can use finite terms of the above series. Under the assumptions of Corollary 6, we will provide the error of the associated approximation.

Lemma 7 (see [20]). Suppose that 𝑔(𝑥) be an enough smooth function in the interval and be approximated by the Bernoulli polynomials as done in Corollary 6. With more details assume that 𝑃_𝑁[𝑔](𝑥) is the approximate polynomial of 𝑔(𝑥) in terms of the Bernoulli polynomials and 𝑅_𝑁[𝑔](𝑥) is the remainder term. Then, the associated formulas are stated as follows: 𝑔(𝑥)=𝑃_𝑁𝑔(𝑥)+𝑅_𝑁𝑔𝑥)[],𝑃(,𝑥∈0,1_𝑁𝑔(𝑥)=¹₀𝑔(𝑥)+𝑑𝑥_𝑁𝑗=1𝐵_𝑗(𝑥)𝑔𝑗!^(𝑗−1)(1)−𝑔^(𝑗−1)(0)𝑅_𝑁𝑔𝑥)1(=−𝑁!¹₀𝐵^∗_𝑁(𝑥−𝑡)𝑔^(𝑁)(𝑡)𝑑𝑡, (39) where 𝐵^∗_𝑁(𝑥)=𝐵_𝑁(𝑥−[𝑥])  and  [𝑥] denotes the largest integer not greater than 𝑥.

Proof. See [20].

Lemma 8. Suppose 𝑔(𝑥)∈C^∞[0,1] (with bounded derivatives) and 𝑔_𝑁(𝑥) is the approximated polynomial using Bernoulli polynomials. Then the error bound would be obtained as follows: 𝐸𝑔_𝑁(𝑥)_∞≤𝐶𝐺(2𝜋)^−𝑁,𝑥∈[0,1], (40) where 𝐺 denotes a bound for all the derivatives of function 𝑔(𝑥) (i.e., ‖𝑔^(𝑖)(𝑥)‖_∞≤𝐺, for 𝑖=0,1,…) and 𝐶 is a positive constant.

Proof. By using Lemma 7, we have 𝐸𝑔_𝑁(𝑥)_∞=1𝑁!¹₀𝐵^∗_𝑁(𝑥−𝑡)𝑔^(𝑁)(𝑡)𝑑𝑡_∞≤1𝐺𝐵𝑁!^𝑁(𝑥)_∞. (41)
According to [20] one can write 𝐵_𝑁(𝑥)=^𝑁_𝑛=0𝑁𝑛𝐵_𝑛(0)𝑥^𝑁−𝑛=^[𝑁/2]_𝑙=0𝑁𝐵2𝑙^2𝑙(0)𝑥^𝑁−2𝑙−12𝑛1𝑥^𝑁−1,𝑥∈[0,1]. (42)
Now we use the formulae (1.1.5) in [27] for the even Bernoulli numbers as follows: 𝐵_2𝑙(0)≤2(2𝑙)!(2𝜋)^−2𝑙. (43)
Therefore, 𝐵_𝑁(𝑥)_∞≤2𝑁!^[𝑁/2]_𝑙=0(2𝜋)^−2𝑙()!+𝑁𝑁−2𝑙2=2𝑁!(2𝜋)^{−𝑁[𝑁/2]}_𝑙=0(2𝜋)^𝑁−2𝑙()!+𝑁𝑁−2𝑙2≤2𝑁!(2𝜋)^−𝑁()+𝑁exp2𝜋2. (44)
In other words ‖𝐵_𝑁(𝑥)‖_∞≤𝐶𝑁!(2𝜋)^−𝑁, where 𝐶 is a positive constant independent of 𝑁. This completes the proof.

Corollary 9. Assume that 𝑢(𝑥,𝑦)∈𝐻×𝐻=𝐿²[0,1]×𝐿²[0,1] be an enough smooth function and also is approximated by the two variable Bernoulli series ∑^∞_𝑚=0∑^∞_𝑛=0𝑢_𝑚,𝑛𝐵_𝑚(𝑥)𝐵_𝑛(𝑦), then the coefficients 𝑢_𝑚,𝑛 for all 𝑚,𝑛=0,1,…,∞ can be calculated from the following relation: 𝑢_𝑚,𝑛=1𝑚!𝑛!¹₀¹₀𝜕^𝑚+𝑛𝑢𝑥,𝑦𝜕𝑥^𝑚𝜕𝑦^𝑛𝑑𝑥𝑑𝑦. (45)

Proof. By applying a similar procedure in two variables (which was provided in Corollary 6) we can conclude the desired result.

In [28], a generalization of Lemma 7 can be found. Therefore, we just recall the error of the associated approximation in two dimensional functions.

Lemma 10. Suppose that 𝑢(𝑥,𝑦) be an enough smooth function and 𝑢_𝑁(𝑥,𝑦) be the approximated polynomial of 𝑢(𝑥,𝑦) in terms of linear combination of Bernoulli polynomials by the aid of Corollary 9. Then the error bound would be obtained as follows: 𝐸𝑢_𝑁𝑥,𝑦_∞∶=𝑢𝑥,𝑦−𝑢_𝑁𝑥,𝑦_∞≤𝐶𝑈𝑁(2𝜋)^−𝑁, (46) where 𝑈 is a positive constant independent of 𝑁 and is a bound for all the partial derivatives of 𝑢(𝑥,𝑦).

We now consider the basic equation (1) together with the initial conditions (2). For clarity of presentation we assume that 𝑚=1,𝑃(𝑧)=𝑃₀(𝑧), and 𝛽=𝛽₀. A similar procedure would be applied in the case of larger values of 𝑚. Equation (1) with the above-mentioned assumptions has the following form: 𝑓^′(𝑧)+𝑃(𝑧)𝑓(𝑧)=𝐺(𝑧), (47) where 𝑓(𝑧)(𝑧),𝐺(𝑧).=𝑢𝑥,𝑦+𝑖𝑣𝑥,𝑦,𝑃=𝑝𝑥,𝑦+𝑖𝑞𝑥,𝑦=𝑔𝑥,𝑦+𝑖ℎ𝑥,𝑦 (48) Also the initial condition 𝑓(0)=𝑢(0,0)+𝑖𝑣(0,0)=𝛽=𝛽₁+𝑖𝛽₂ should be considered. Thus, one can write 𝑢(0,0)=𝛽₁ and 𝑣(0,0)=𝛽₂.

According to Remark 5, since 𝑓(𝑧) is a Holomorphic function, then (47) by using assumptions in (48) could be rewritten as follows: 𝑢_𝑥𝑥,𝑦+𝑖𝑣_𝑥+𝑝𝑢.𝑥,𝑦𝑥,𝑦+𝑖𝑞𝑥,𝑦𝑥,𝑦+𝑖𝑣𝑥,𝑦=𝑔𝑥,𝑦+𝑖ℎ𝑥,𝑦 (49) In other words 𝑢_𝑥+𝑝𝑢𝑣,𝑣𝑥,𝑦𝑥,𝑦𝑥,𝑦−𝑞𝑥,𝑦𝑥,𝑦=𝑔𝑥,𝑦_𝑥𝑥,𝑦+𝑞𝑥,𝑦𝑢𝑥,𝑦+𝑝𝑥,𝑦𝑣𝑥,𝑦=ℎ𝑥,𝑦. (50) For imposing the initial conditions 𝑢(0,0)=𝛽₁ and 𝑣(0,0)=𝛽₂, we need to differentiate the above equations with respect to 𝑦 and then integrate with respect to 𝑥 and 𝑦 in the rectangular [0,𝑥]×[0,𝑦]. Therefore, by differentiating both of the equations of (50) with respect to 𝑦 we have 𝑢_𝑥𝑦+𝜕𝑥,𝑦𝑝𝑢𝑣𝜕𝑦𝑥,𝑦𝑥,𝑦−𝑞𝑥,𝑦𝑥,𝑦=𝑔_𝑦,𝑣𝑥,𝑦_𝑥𝑦+𝜕𝑥𝑦𝑞𝑢𝑣𝜕𝑦𝑥,𝑦𝑥,𝑦+𝑝𝑥,𝑦𝑥,𝑦=ℎ_𝑦𝑥,𝑦. (51) Integrating from the above equations in the rectangular [0,𝑥]×[0,𝑦] yields ^𝑥₀^𝑦₀𝑢_𝑥𝑦+𝑥,𝑦𝑑𝑦𝑑𝑥_𝑥0𝑦0𝜕𝑝𝑢𝑣=𝜕𝑦𝑥,𝑦𝑥,𝑦−𝑞𝑥,𝑦𝑥,𝑦𝑑𝑦𝑑𝑥_𝑥0𝑦0𝑔_𝑦𝑥,𝑦𝑑𝑦𝑑𝑥,_𝑥0𝑦0𝑣_𝑥𝑦+𝑥,𝑦𝑑𝑦𝑑𝑥_𝑥0𝑦0𝜕𝑞𝑢𝑣=𝜕𝑦𝑥,𝑦𝑥,𝑦+𝑝𝑥,𝑦𝑥,𝑦𝑑𝑦𝑑𝑥_𝑥0𝑦0ℎ_𝑦𝑥,𝑦𝑑𝑦𝑑𝑥. (52) By imposing the the initial conditions 𝑢(0,0)=𝛽₁ and 𝑣(0,0)=𝛽₂ to these equations we reach to 𝑢𝑥,𝑦+^𝑥₀𝑝𝑢𝑣−𝑝()𝑢()()𝑣()=𝑔,𝑣+𝑥,𝑦𝑥,𝑦−𝑞𝑥,𝑦𝑥,𝑦𝑥,0𝑥,0−𝑞𝑥,0𝑥,0𝑑𝑥𝑥,𝑦𝑥,𝑦^𝑥₀𝑞𝑢𝑣−𝑞()𝑢()()𝑣(=ℎ,𝑥,𝑦𝑥,𝑦+𝑝𝑥,𝑦𝑥,𝑦𝑥,0𝑥,0+𝑝𝑥,0𝑥,0)𝑑𝑥𝑥,𝑦 (53) where 𝑔𝑥,𝑦=^𝑥₀^𝑦₀𝑔_𝑦𝑥,𝑦𝑑𝑦𝑑𝑥+𝛽₁,ℎ𝑥,𝑦=^𝑥₀^𝑦₀ℎ_𝑦𝑥,𝑦𝑑𝑦𝑑𝑥+𝛽₂. (54) Considering 𝑈(𝑥,𝑦)=^{𝑢(𝑥,𝑦)}_{𝑣(𝑥,𝑦)},𝑔𝐺(𝑥,𝑦)=^(𝑥,𝑦)ℎ_(𝑥,𝑦) and 𝐻(𝑥,𝑦)=^{𝑝(𝑥,𝑦)−𝑞(𝑥,𝑦)}_{𝑞(𝑥,𝑦)𝑝(𝑥,𝑦)}, (53) could be restated in the following matrix vector form: 𝑈𝑥,𝑦+^𝑥₀𝐻𝑈−𝑥,𝑦𝑥,𝑦𝑑𝑥_𝑥0𝐻(𝑥,0)𝑈(𝑥,0)𝑑𝑥=𝐺𝑥,𝑦. (55) Our aim is to show that lim_𝑁→∞𝑈_𝑁(𝑥,𝑦)=𝑈(𝑥,𝑦) or lim_𝑁→∞𝑓_𝑁(𝑧)=𝑓(𝑧), where 𝑈_𝑁(𝑥,𝑦)=^𝑢𝑁_𝑣^(𝑥,𝑦)_{𝑁(𝑥,𝑦)}, 𝑓_𝑁(𝑧)=𝑢_𝑁(𝑥,𝑦)+𝑖𝑣_𝑁(𝑥,𝑦) and 𝑓_𝑁(𝑧) was introduced in (3).

In the following lines, the main theorem of this section would be provided. However, some mild conditions should be assumed. These conditions are as follows:(i)‖𝐻(𝑥,𝑦)‖_∞≤𝐻≪1,where𝐻isarealconstant, (ii)lim_𝑁→∞‖𝑈(𝑥,0)−𝑈_𝑁(𝑥,0)‖_∞≤lim_𝑁→∞(1/𝑁)‖𝑈(𝑥,𝑦)−𝑈_𝑁(𝑥,𝑦)‖_∞.

It should be noted that the second condition is based upon Lemmas 8 and 10.

Theorem 11. Assume that 𝑈 and 𝑈₀∶=𝑈(𝑥,0) be approximated by 𝑈_𝑁 and 𝑈_0,𝑁, respectively, by the aid of the Bernoulli polynomials in (55) and also we use a collocation scheme for providing the numerical solution of (55). In other words 𝑈_𝑁𝑥,𝑦+^𝑥₀𝐻𝑥,𝑦𝑈_𝑁−𝑥,𝑦𝑑𝑥_𝑥0𝐻(𝑥,0)𝑈_𝑁()=𝐺𝑥,0𝑑𝑥𝑥,𝑦+𝑅_𝑁+1𝑥,𝑦, (56) where 𝑅_𝑁+1(𝑥,𝑦) is the residual function that is zero at the (𝑁+1) collocation nodes. Also suppose that [𝑎,𝑏]=[𝑐,𝑑]=[0,1]. Then, under the above-mentioned assumptions 𝑈−𝑈_𝑁∞≤𝑅_𝑁+1𝑥,𝑦_∞𝐻(),1−1+1/𝑁 (57) and lim_𝑁→∞𝑈_𝑁(𝑥,𝑦)=𝑈(𝑥,𝑦).

Proof. By subtracting (56) from (55) we have 𝑈−𝑈_𝑁∞=^𝑥₀−𝐻𝑈+𝐻₀𝑈₀+𝐻𝑈_𝑁−𝐻₀𝑈_0,𝑁𝑑𝑥+𝑅_𝑁+1𝑥,𝑦_∞≤^𝑥₀𝐻𝑈−𝐻𝑈_𝑁𝑑𝑥_∞+^𝑥₀𝐻₀𝑈₀−𝐻₀𝑈_0,𝑁𝑑𝑥_∞+𝑅_𝑁+1𝑥,𝑦_∞≤^𝑥₀𝐻𝑈−𝑈_𝑁∞+𝐻𝑈₀−𝑈_0,𝑁∞+𝑅𝑑𝑥_𝑁+1𝑥,𝑦_∞≤𝐻𝑈−𝑈_𝑁∞11+𝑁+𝑅_𝑁+1𝑥,𝑦_∞. (58) Therefore, 𝑈−𝑈_𝑁∞≤𝑅_𝑁+1𝑥,𝑦_∞𝐻().1−1+1/𝑁 (59)
Since lim_𝑁→∞‖𝑅_𝑁+1(𝑥,𝑦)‖_∞=0and𝐻≪1, we have lim_𝑁→∞‖𝑈−𝑈_𝑁‖_∞=0; in other words, lim_𝑁→∞𝑓_𝑁(𝑧)=𝑓(𝑧) and this completes the proof.

5. Numerical Examples

In this section, several numerical examples are given to illustrate the accuracy and effectiveness of the proposed method and all of them are performed on a computer using programs written in MATLAB 7.12.0 (v2011a) (The Mathworks Inc., Natick, MA, USA). In this regard, we have reported in the tables and figures the values of the exact solution 𝑓(𝑧), the polynomial approximate solution 𝑓_𝑁(𝑧), and the absolute error function 𝑒_𝑁(𝑧)=|𝑓(𝑧)−𝑓_𝑁(𝑧)| at the selected points of the given domains. It should be noted that in the first example we consider a complex differential equation with an exact polynomial solution. Our method obtains such exact polynomial solutions readily by solving the associated linear algebraic system.

Example 12 (see [26]). As the first example, we consider the first-order complex differential equation with variable coefficients 𝑓^′(𝑧)+𝑧𝑓(𝑧)=2𝑧²[][],−𝑧+2,𝑧=𝑥+𝑖𝑦,𝑥∈0,1,𝑦∈0,1 (60) with the initial condition 𝑓(0)=−1 and the exact solution 2𝑧−1. We suppose that 𝑁=3. Therefore, the collocation points are 𝑧₀₀=0,𝑧₁₁=(1+𝑖)/3,𝑧₂₂=(2+2𝑖)/3, and 𝑧₃₃=1+𝑖. According to (3), the approximate solution has the following form: 𝑓₃(𝑧)=³_𝑛=0𝑓_𝑛𝐵_𝑛(𝑧),𝑧∈𝐷, (61) where our aim is to find the unknown Bernoulli coefficients 𝑓₀,𝑓₁,𝑓₂,and𝑓₃. Since 𝑃₀(𝑧)=𝑧, then 𝑃₀=.000000.3333+0.3333𝑖00000.6667+0.6667𝑖00001.0000+1.0000𝑖 (62) Also the matrix 𝐿 (with the assumption 𝑁=3) has the following structure:𝐵𝑧𝐿=⁰⁰𝐵𝑧₁₁𝐵𝑧₂₂𝐵𝑧₃₃=𝐵₀𝑧₀₀𝐵₁𝑧₀₀𝐵₂𝑧₀₀𝐵₃𝑧₀₀𝐵₀𝑧₁₁𝐵₁𝑧₁₁𝐵₂𝑧₁₁𝐵₃𝑧₁₁𝐵₀𝑧₂₂𝐵₁𝑧₂₂𝐵₂𝑧₂₂𝐵₃𝑧₂₂𝐵₀𝑧₃₃𝐵₁𝑧₃₃𝐵₂𝑧₃₃𝐵₃𝑧₃=,1.0000−0.50000.166701.0000−0.1667+0.3333𝑖−0.1667−0.1111𝑖0.0926−0.0926𝑖1.00000.1667+0.6667𝑖−0.5000+0.2222𝑖−0.2593−0.4074𝑖1.00000.5000+1.0000𝑖−0.8333+1.0000𝑖−1.5000−0.5000𝑖 (63)where 𝐵₀(𝑧)=1,𝐵₁(𝑧)=𝑧−1/2,𝐵₂(𝑧)=𝑧²−𝑧+1/6and𝐵₃(𝑧)=𝑧³−(3/2)𝑧²+𝑧/2.
According to (29), the matrix coefficients 𝑊 are as follows:𝑊=𝐿𝑀^𝑇+𝑃₀𝐿=.01.0000−1.00000.50000.3333+0.3333𝑖0.8333+0.0556𝑖−0.3519+0.5741𝑖−0.4383−0.3333𝑖0.6667+0.6667𝑖0.6667+0.5556𝑖−0.1481+1.1481𝑖−1.4012+0.2222𝑖1.0000+1.0000𝑖0.5000+1.5000𝑖−0.8333+2.1667𝑖−3.5000+1.0000𝑖 (64)Since 𝐺(𝑧)=2𝑧²−𝑧+2, the right-hand side vector 𝐺 has the following form: 𝐺𝑧𝐺=⁰⁰𝐺𝑧₁₁𝐺𝑧₂₂𝐺𝑧₃₃=2.00001.6667+0.1111𝑖1.3333+1.1111𝑖1.0000+3.0000𝑖. (65) The associated form of the initial condition 𝑓(0)=−1 is 𝑈₀;𝛽₀=[𝐵(0)(0=.;𝐺)]1.0000−0.50000.16670;−1 (66) Imposing the above initial condition to the matrix 𝑊 and vector 𝐺 yields,.𝑊=01.0000−1.00000.50000.3333+0.3333𝑖0.8333+0.0556𝑖−0.3519+0.5741𝑖−0.4383−0.3333𝑖0.6667+0.6667𝑖0.6667+0.5556𝑖−0.1481+1.1481𝑖−1.4012+0.2222𝑖1.0000−0.50000.16670𝐺=2.00001.6667+0.1111𝑖1.3333+1.1111𝑖−1 (67)Therefore, the solution of the system 𝑊𝐹=𝐺 is as follows: 𝐹=0.0000−0.0000𝑖2.0000−0.0000𝑖0.0000−0.0000𝑖0.0000+0.0000𝑖. (68) Thus, we obtain the approximate solution 𝑓₃(𝑧)=𝐵(𝑧)𝐹=2𝑧−1 which is the exact solution. We recall that 𝐵(𝑧)=11𝑧−2𝑧²1−𝑧+6𝑧³−32𝑧²+𝑧2. (69)

Example 13 (see [26]). As the second example, we consider the following second-order complex differential equation 𝑓^′′(𝑧)+𝑧𝑓^′(𝑧)(𝑧)(𝑧)(𝑧),[][],+𝑧𝑓=exp+2𝑧exp𝑧=𝑥+𝑖𝑦,𝑥∈0,1,𝑦∈0,1 (70) with the initial conditions 𝑓(0)=𝑓^′(0)=1 and the exact solution 𝑓(𝑧)=exp(𝑧). In this equation, we have 𝑃₀(𝑧)=𝑃₁(𝑧)=𝑧 and 𝐺(𝑧)=exp(𝑧)+2𝑧exp(𝑧). Then, for 𝑁=5 the collocation points are 𝑧₀₀=0,𝑧₁₁=(1+𝑖)/5,𝑧₂₂=(2+2𝑖)/5,𝑧₃₃=(3+3𝑖)/5,𝑧₄₄=(4+4𝑖)/5,and𝑧₅₅=1+𝑖.
According to (29), the fundamental matrix equation is 𝑊=𝐿𝑀^𝑇2+𝑃₁𝐿𝑀^𝑇+𝑃₀𝐿, (71) where 𝑀^𝑇 and 𝐿 are introduced in (18) and (25), respectively. Also 𝑃₀=𝑃₁=.000000.3333+0.3333𝑖00000.6667+0.6667𝑖00001.0000+1.0000𝑖 (72) The right-hand side vector 𝐺 is 𝐺=1.00001.5788+0.8185𝑖2.0086+2.1449𝑖2.0739+4.0681𝑖1.4770+6.6318𝑖−0.1686+9.7995𝑖. (73) The augmented matrix forms of the initial conditions for 𝑁=5 are [𝐵(0)]=,𝐵(0)𝑀;11.0000−0.50000.16670−0.03330;1_𝑇=.;101.0000−1.00000.50000−0.1667;1 (74) By replacing the above augmented vectors to the last two rows of [𝑊;𝐺], we reach to [𝑊;𝐺]. The solution of the matrix-vector equation 𝑊𝐹=𝐺 is 𝐹=1.7184−0.0000𝑖1.7189+0.0003𝑖0.8603+0.0021𝑖0.2864+0.0046𝑖0.0693+0.0048𝑖0.0107+0.0034𝑖. (75) We also solve this equation by assumptions 𝑁=7 and 𝑁=9. Since the exact solution of the equation is exp(𝑧)=exp(𝑥+𝑖𝑦)=exp(𝑥)(cos(𝑦)+𝑖sin(𝑦)). Therefore, the real and imaginary parts of the exact solution are Re(exp(𝑧))=exp(𝑥)cos(𝑦) and Im(exp(𝑧))=exp(𝑥)sin(𝑦), respectively. The values of the approximate solution in the case of 𝑁=5,7,and9 for both parts of real and imaginary together with the exact solution are provided in Tables 1 and 2. Also an interesting comparison between the presented method (PM) and the Taylor method [26] (TM) has been provided in Figure 1 for 𝑁=7and9. From this figure, one can see the efficiency of our methods with respect to the method of [26].

Table 1 

Numerical results of Example 13, for Re(exp(𝑧))=exp(𝑥)cos(𝑦).

Table 2 

Numerical results of Example 13, for Im(exp(𝑧)=exp(𝑧)sin(𝑦).

Figure 1 

Comparisons of the presented method (PM) and the Taylor method (TM) of Example 13.

Example 14 (see [26]). As the final example, we consider the following second-order complex differential equation 𝑓^′′(𝑧)+𝑧𝑓^′(𝑧)(𝑧)(𝑧)(𝑧)(𝑧)[][],+2𝑧𝑓=2𝑧sin+𝑧cos−sin,𝑧=𝑥+𝑖𝑦,𝑥∈0,1,𝑦∈0,1 (76) with the initial conditions 𝑓(0)=0 and 𝑓^′(0)=1 and also the exact solution 𝑓(𝑧)=sin(𝑧). In this equation, we have 𝑃₀(𝑧)=2𝑧,𝑃₁(𝑧)=𝑧, and 𝐺(𝑧)=2𝑧sin(𝑧)+𝑧cos(𝑧)−sin(𝑧). Similar to the previous two examples, we solve this equation in the case of 𝑁=9 and 11. In Figure 2, we provide a comparison between our method and the Taylor method (TM) [26]. According to this figure, one can see that not only our method is superior in results but also the behaviour of the error of the presented method has a stable manner with respect to the Taylor method [26] during the computational interval. Moreover, since the Bessel method and the presented method can be considered as an preconditioned solution of any linear algebraic system which is originated from the above differential equation, they have the same accuracy. But the condition number of the matrix coefficients 𝑊 of the Bessel method is very larger than the matrix coefficient 𝑊 of our method. This subject is illustrated in Figure 3. It should be noted that this figure is depicted according to the diagonal collocation methods and the square collocation method was not used.

Figure 2 

Comparisons of the presented method (PM) and the Taylor method (TM) of Example 14.

Figure 3 

Condition number comparisons between our method and the Bessel method of Example 14.

6. Conclusions

High-order linear complex differential equations are usually difficult to solve analytically. Then, it is required to obtain the approximate solutions. For this reason, a new technique using the Bernoulli polynomials to numerically solve such equations is proposed. This method is based on computing the coefficients in the Bernoulli series expansion of the solution of a linear complex differential equation and is valid when the functions 𝑃_𝑘(𝑧) and 𝐺(𝑧) are defined in the rectangular domain. An interesting feature of this method is to find the analytical solutions if the equation has an exact solution that is a polynomial of degree 𝑁 or less than 𝑁. Shorter computation time and lower operation count results in reduction of cumulative truncation errors and improvement of overall accuracy are some of the advantages of our method. In addition, the method can also be extended to the system of linear complex equations with variable coefficients, but some modifications are required.

Conflict of Interests

The authors declare that they do not have any conflict of interests in their submitted paper.

Acknowledgments

The authors thank the Editor and both of the reviewers for their constructive comments and suggestions to improve the quality of the paper.