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

Authors: Faezeh Toutounian; Emran Tohidi; Stanford Shateyi
Journal: Abstract and Applied Analysis
 (2013)

---

## 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.


---

## Body

## 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)f(m)(z)+∑k=0m-1‍Pk(z)f(k)(z)=G(z),m≥1,z=x+iy,x∈[a,b],y∈[c,d],
with the initial conditions
(2)f(r)(0)=βr,r=0,1,…,m-1.
Evidently, we will let z denote a variable point in the complex plane. Its real and imaginary parts will be denoted by x and y, respectively.Here, the coefficientsPk(z) and the known function G(z) together with the unknown function f(z) are holomorphic (or analytic) functions in the domain D={z∈ℂ:z=x+iy,  a≤x≤b,  c≤y≤d:a,b,c,d∈ℝ} where the coefficients βr are appropriate complex constants.We assume that the solution of (1) under the conditions (2) is approximated in the form
(3)f(z)≈fN(z)=∑n=0N‍fnBn(z),z∈D,
which is the truncated Bernoulli series of the unknown function f(z), where all of fn(n=0,1,…,N) are the Bernoulli coefficients to be determined. We also use the collocation points
(4)zpp=xp+iyp,xp=a+b-aNp,zyp=c+d-cNp,p=0,1,…,N.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 Section2, 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 functionf(z)=z) 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 n, the function f(z)=zn is differentiable everywhere, with derivative nzn-1. This, in conjunction with the sum and product rules, implies that every polynomial is everywhere differentiable: if f(z)=cnzn+cn-1zn-1+⋯+c1z+c0, where c0,…,cn are complex constants, then f′(z)=ncnzn-1+(n-1)cn-1zn-2+⋯+c1.Remark 1.
The functionf=u+iv is differentiable (in the complex sense) at z0 if and only if u and v are differentiable (in the real sense) at z0 and their first partial derivatives satisfy the relations ∂u(z0)/∂x=∂v(z0)/∂y,∂u(z0)/∂y=-∂v(z0)/∂x. In that case
(5)f′(z0)=∂u(z0)∂x+i∂v(z0)∂x=∂v(z0)∂y-i∂u(z0)∂y.Remark 2.
Two partial differential equations(6)∂u∂x=∂v∂y,∂u∂y=-∂v∂x
are called the Cauchy-Riemann equations for the pair of functions u,v. 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 functionf=u+iv be defined in the open subset G of the complex plane ℂ and assume that u and v have first partial derivatives in G. Then f 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 subsetG of the complex plane ℂ and differentiable at every point of G is said to be holomorphic (or analytic) in G. 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∂/∂x and ∂/∂y are applied to a complex-valued function f=u+iv in the natural way:
(7)∂f∂x=∂u∂x+i∂v∂x,∂f∂y=∂u∂y+i∂v∂y.
We define the complex partial differential operators ∂/∂z and ∂/∂z- by
(8)∂∂z=12(∂∂x-i∂∂y),∂∂z-=12(∂∂x+i∂∂y).
Thus, ∂/∂x=∂/∂z+∂/∂z-and∂/∂y=i(∂/∂z-∂/∂z-).
Intuitively one can think of a holomorphic function as a complex-valued function in an open subset ofℂ that depends only on z, that is, independent of z-. We can make this notion precisely as follows. Suppose the function f=u+iv is defined and differentiable in an open set. One then has
(9)∂f∂z=12(∂u∂x+∂v∂y)+i2(∂v∂x-∂u∂y),∂f∂z-=12(∂u∂x-∂v∂y)+i2(∂v∂x+∂u∂y).The Cauchy-Riemann equations thus can be written∂f/∂z-=0. As this is the condition for f to be holomorphic, it provides a precise meaning for the statement: a holomorphic function is one that is independent ofz-. If f is holomorphic, then (not surprisingly) f′=∂f/∂z, as the following calculation shows:
(10)f′=∂f∂x=∂f∂z+∂f∂z-=∂f∂z.


### 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 polynomialsBn(x) are usually defined by means of the exponential generating functions (see [21])
(11)wexwew-1=∑k=0∞‍Bk(x)wkk!,(|w|<2π).
The following familiar expansion (see [20]):
(12)∑k=0n‍(n+1k)Bk(x)=(n+1)xn
is the most primary property of the Bernoulli polynomials. The first few Bernoulli polynomials are
(13)B0(x)=1,B1(x)=x-12,B2(x)=x2-x+16,B3(x)=x3-32x2+12x,B4(x)=x4-2x3+x2-130.
The Bernoulli polynomials satisfy the well-known relations (see [21])
(14)dBn(x)dx=nBn-1(x),(n≥1),∫01Bn(x)dx=0,(n≥1).The Bernoulli polynomials have another specific property that is satisfied in the following linear homogeneous recurrence relation [20]:
(15)Bn(x)=(x-12)Bn-1(x)-1n∑k=0n-2‍(nk)Bn-k(0)Bk(x).
Also, the Bernoulli polynomials satisfy the following interesting property [25]:
(16)∫01Bm(x)Bn(x)dx=(-1)m+nm!n!(m+n)!Bm+n(0).
Moreover, Bn(x) satisfy the differential equation [20]
(17)Bn(0)n!y(n)(x)+Bn-1(0)(n-1)!y(n-1)(x)+⋯+B2(0)2!y′′(x)+(12-x)y′(x)+ny(x)=0.
According to the discussions in [20], the Bernoulli polynomials form a complete basis over the interval [0,1].If we introduce the Bernoulli vectorB(x) in the form B(x)=[B0(x),B1(x),…,BN(x)], then the derivative of B(x), with the aid of the first property of (14), can be expressed in the matrix form by
(18)[B0′(x)B1′(x)B2′(x)⋮BN-1′(x)BN′(x)]︸B′(x)T=[000⋯000100⋯000020⋯000⋮⋮⋮⋱⋮⋮⋮000…N-100000…0N0]︸M[B0(x)B1(x)B2(x)⋮BN-1(x)BN(x)]︸BT(x),
where M is the (N+1)×(N+1) operational matrix of differentiation. Note that if we replace the real variable x by the complex variable z in the above relation, again we reach to the same result, since (zn)′=nzn-1.Accordingly, thekth derivative of B(x) can be given by
(19)B′(x)T=MB(x)T⟹B(1)(x)=B(x)MT,B(2)(x)=B(1)(x)MT=B(x)(MT)2,B(3)(x)=B(1)(x)(MT)2=B(x)(MT)3,⋮B(k)(x)=B(x)(MT)k,ghol;iogll;lghuli
where M 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,B6(x) (the 6th Bernoulli polynomial) has 5 terms while P6(x) (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.


## 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 functionf(z)=z) 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 n, the function f(z)=zn is differentiable everywhere, with derivative nzn-1. This, in conjunction with the sum and product rules, implies that every polynomial is everywhere differentiable: if f(z)=cnzn+cn-1zn-1+⋯+c1z+c0, where c0,…,cn are complex constants, then f′(z)=ncnzn-1+(n-1)cn-1zn-2+⋯+c1.Remark 1.
The functionf=u+iv is differentiable (in the complex sense) at z0 if and only if u and v are differentiable (in the real sense) at z0 and their first partial derivatives satisfy the relations ∂u(z0)/∂x=∂v(z0)/∂y,∂u(z0)/∂y=-∂v(z0)/∂x. In that case
(5)f′(z0)=∂u(z0)∂x+i∂v(z0)∂x=∂v(z0)∂y-i∂u(z0)∂y.Remark 2.
Two partial differential equations(6)∂u∂x=∂v∂y,∂u∂y=-∂v∂x
are called the Cauchy-Riemann equations for the pair of functions u,v. 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 functionf=u+iv be defined in the open subset G of the complex plane ℂ and assume that u and v have first partial derivatives in G. Then f 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 subsetG of the complex plane ℂ and differentiable at every point of G is said to be holomorphic (or analytic) in G. 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∂/∂x and ∂/∂y are applied to a complex-valued function f=u+iv in the natural way:
(7)∂f∂x=∂u∂x+i∂v∂x,∂f∂y=∂u∂y+i∂v∂y.
We define the complex partial differential operators ∂/∂z and ∂/∂z- by
(8)∂∂z=12(∂∂x-i∂∂y),∂∂z-=12(∂∂x+i∂∂y).
Thus, ∂/∂x=∂/∂z+∂/∂z-and∂/∂y=i(∂/∂z-∂/∂z-).
Intuitively one can think of a holomorphic function as a complex-valued function in an open subset ofℂ that depends only on z, that is, independent of z-. We can make this notion precisely as follows. Suppose the function f=u+iv is defined and differentiable in an open set. One then has
(9)∂f∂z=12(∂u∂x+∂v∂y)+i2(∂v∂x-∂u∂y),∂f∂z-=12(∂u∂x-∂v∂y)+i2(∂v∂x+∂u∂y).The Cauchy-Riemann equations thus can be written∂f/∂z-=0. As this is the condition for f to be holomorphic, it provides a precise meaning for the statement: a holomorphic function is one that is independent ofz-. If f is holomorphic, then (not surprisingly) f′=∂f/∂z, as the following calculation shows:
(10)f′=∂f∂x=∂f∂z+∂f∂z-=∂f∂z.


## 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 polynomialsBn(x) are usually defined by means of the exponential generating functions (see [21])
(11)wexwew-1=∑k=0∞‍Bk(x)wkk!,(|w|<2π).
The following familiar expansion (see [20]):
(12)∑k=0n‍(n+1k)Bk(x)=(n+1)xn
is the most primary property of the Bernoulli polynomials. The first few Bernoulli polynomials are
(13)B0(x)=1,B1(x)=x-12,B2(x)=x2-x+16,B3(x)=x3-32x2+12x,B4(x)=x4-2x3+x2-130.
The Bernoulli polynomials satisfy the well-known relations (see [21])
(14)dBn(x)dx=nBn-1(x),(n≥1),∫01Bn(x)dx=0,(n≥1).The Bernoulli polynomials have another specific property that is satisfied in the following linear homogeneous recurrence relation [20]:
(15)Bn(x)=(x-12)Bn-1(x)-1n∑k=0n-2‍(nk)Bn-k(0)Bk(x).
Also, the Bernoulli polynomials satisfy the following interesting property [25]:
(16)∫01Bm(x)Bn(x)dx=(-1)m+nm!n!(m+n)!Bm+n(0).
Moreover, Bn(x) satisfy the differential equation [20]
(17)Bn(0)n!y(n)(x)+Bn-1(0)(n-1)!y(n-1)(x)+⋯+B2(0)2!y′′(x)+(12-x)y′(x)+ny(x)=0.
According to the discussions in [20], the Bernoulli polynomials form a complete basis over the interval [0,1].If we introduce the Bernoulli vectorB(x) in the form B(x)=[B0(x),B1(x),…,BN(x)], then the derivative of B(x), with the aid of the first property of (14), can be expressed in the matrix form by
(18)[B0′(x)B1′(x)B2′(x)⋮BN-1′(x)BN′(x)]︸B′(x)T=[000⋯000100⋯000020⋯000⋮⋮⋮⋱⋮⋮⋮000…N-100000…0N0]︸M[B0(x)B1(x)B2(x)⋮BN-1(x)BN(x)]︸BT(x),
where M is the (N+1)×(N+1) operational matrix of differentiation. Note that if we replace the real variable x by the complex variable z in the above relation, again we reach to the same result, since (zn)′=nzn-1.Accordingly, thekth derivative of B(x) can be given by
(19)B′(x)T=MB(x)T⟹B(1)(x)=B(x)MT,B(2)(x)=B(1)(x)MT=B(x)(MT)2,B(3)(x)=B(1)(x)(MT)2=B(x)(MT)3,⋮B(k)(x)=B(x)(MT)k,ghol;iogll;lghuli
where M 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,B6(x) (the 6th Bernoulli polynomial) has 5 terms while P6(x) (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 fN(z) in the form (3). Trivially, fN(z) could be rewritten in the vector form
(20)fN(z)=B(z)F,B(z)=[B0(z)B1(z)⋯BN(z)],F=[f0f1⋯fN]T.
By using (19) in the complex form one can conclude that
(21)fN(k)(z)=B(z)(MT)kF,k≤N,
where M is introduced in (19).For the collocation pointsz=zpp(p=0,1,…,N), the matrix relation (21) becomes
(22)fN(k)(zpp)=B(zpp)(MT)kF,p=0,1,…,N,
where B(zpp)=[B0(zpp)B1(zpp)⋯BN(zpp)]. For more details, one can restate (22) as follows:
(23)fN(k)(z00)=B(z00)(MT)kF,fN(k)(z11)=B(z11)(MT)kF,⋮fN(k)(zNN)=B(zNN)(MT)kF.
The matrix vector of the above-mentioned equations is
(24)F(k)=[fN(k)(z00)fN(k)(z11)⋮fN(k)(zNN)]=L(MT)kF,
where
(25)L=[B(z00)B(z11)⋮B(zNN)]=[B0(z00)B1(z00)B2(z00)⋯BN(z00)B0(z11)B1(z11)B2(z11)⋯BN(z11)⋮⋮⋮⋱⋮B0(zNN)B1(zNN)B2(zNN)⋯BN(zNN)].
On the other hand by substituting the collocation points z=zpp defined by (4) into (1), we have
(26)fN(m)(zpp)+∑k=0m-1‍Pk(zpp)fN(k)(zpp)=G(zpp),fN(m)(zpp)∑k=0m-1‍(zpp)=G(zpp),p=0,1,…,N.
The associated matrix-vector form of the above equations has the following form by the aid of (24):
(27)F(m)+∑k=0m-1‍PkF(k)=L(MT)mF+∑k=0m-1‍PkL(MT)kF=G,
where
(28)G=[G(z00)G(z11)⋯G(zNN)]T,Pk=[Pk(z00)0⋯00Pk(z11)⋯0⋮⋮⋱⋮00⋯Pk(zNN)].
Since the vector F is unknown and should be determined, therefore, the matrix-vector equation (27) could be rewritten in the following form:
(29)WF=G,or[W;G]=[wpq;gp],WF=G,or[],p,q=0,1,…,N,
where W=L(MT)m+∑k=0m-1‍PkL(MT)k.We now write the vector form of the initial conditions (2) by the aid of (21) as follows:
(30)f(r)(0)=B(0)(MT)rF=βr,r=0,1,…,m-1.
In other words the vector form of the initial conditions could be rewritten as UrF=βr where Ur=B(0)(MT)r  (r=0,1,…,m-1). Trivially the augmented form of these equations is
(31)[Ur;βr]=[ur0ur1⋯urN;βr],r=0,1,…,m-1.
Consequently, to find the unknown Bernoulli coefficients fn,n=0,1,…,N, related with the approximate solution of the problem (1) under the initial conditions (2), we need to replace the m rows of (31) by the last m rows of the augmented matrix (29) and hence we have new augmented matrix
(32)[W^;G^]=[w00w01⋯w0N;g0w10w11⋯w1N;g1⋮⋮⋱⋮⋮⋮wN-m,0wN-m,1⋯wN-m,N;gN-mu00u01⋯u0N;β0u10u11⋯u1N;β1⋮⋮⋱⋮⋮⋮um-1,0um-1,1⋯um-1,N;βm-1],
or the corresponding matrix-vector equation
(33)W^F=G^.
If det(W^)≠0 one can rewrite (32) in the form F=(W^)-1G^ and the vector F is uniquely determined. Thus the mth 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 f(z) and its derivatives are substituted in (1), the resulting equation must be satisfied approximately; that is, for z=zi∈D,i=0,1,2,…(34)E(zi)=|f(m)(zi)+∑k=0m-1‍Pk(zi)f(k)(zi)-G(zi)|≅0,orE(zi)≤10-ki(kiis any positive integer).
If max(10-ki)=10-k (k is any positive integer) is prescribed, then the truncation limit N may be increased until the values E(zi) at each of the points zi become smaller than the prescribed 10-k (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 thatH=L2[0,1] and {B0(x),B1(x),…,BN(x)}⊂H be the set of the Bernoulli polynomials and
(35)Y=span{B0(x),B1(x),…,BN(x)}
and g be an arbitrary element in H. Since Y is a finite dimensional vector space, g has the unique best approximation g^∈Y, such that
(36)∀y∈Y,∥g-g^∥∞≤∥g-y∥∞.
Since g^∈Y, there exists the unique coefficients g0,g1,…,gN such that
(37)g≈g^=∑n=0N‍gnBn(x)=B(x)GT,B(x)=[B0(x),B1(x),…,BN(x)],B(x)[B0(x)]G=[g0,g1,…,gN].Corollary 6.
Assume thatg∈H=L2[0,1] be an enough smooth function and also is approximated by the Bernoulli serie ∑n=0∞‍gnBn(x), then the coefficients gn for all n=0,1,…,∞ can be calculated from the following relation:
(38)gn=1n!∫01‍g(n)(x)dx.Proof.
See [21].In practice one can use finite terms of the above series. Under the assumptions of Corollary6, we will provide the error of the associated approximation.Lemma 7 (see [20]).
Suppose thatg(x) 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 PN[g](x) is the approximate polynomial of g(x) in terms of the Bernoulli polynomials and RN[g](x) is the remainder term. Then, the associated formulas are stated as follows:
(39)g(x)=PN[g](x)+RN[g](x),x∈[0,1],PN[g](x)=∫01‍g(x)dx+∑j=1N‍Bj(x)j!(g(j-1)(1)-g(j-1)(0))RN[g](x)=-1N!∫01‍BN*(x-t)g(N)(t)dt,
where BN*(x)=BN(x-[x])  and  [x] denotes the largest integer not greater than x.Proof.
See [20].Lemma 8.
Supposeg(x)∈C∞[0,1] (with bounded derivatives) and gN(x) is the approximated polynomial using Bernoulli polynomials. Then the error bound would be obtained as follows:
(40)∥E(gN(x))∥∞≤CG^(2π)-N,x∈[0,1],
where G^ denotes a bound for all the derivatives of function g(x) (i.e., ∥g(i)(x)∥∞≤G^, for i=0,1,…) and C is a positive constant.Proof.
By using Lemma7, we have
(41)∥E(gN(x))∥∞=∥1N!∫01‍BN*(x-t)g(N)(t)dt∥∞≤1N!G^∥BN(x)∥∞.
According to [20] one can write
(42)BN(x)=∑n=0N‍(Nn)Bn(0)xN-n=∑l=0[N/2]‍(N2l)B2l(0)xN-2l-12(n1)xN-1,x∈[0,1].
Now we use the formulae (1.1.5) in [27] for the even Bernoulli numbers as follows:
(43)|B2l(0)|≤2(2l)!(2π)-2l.
Therefore,(44)∥BN(x)∥∞≤2N!∑l=0[N/2]‍(2π)-2l(N-2l)!+N2=2N!(2π)-N∑l=0[N/2]‍(2π)N-2l(N-2l)!+N2≤2N!(2π)-Nexp(2π)+N2.
In other words∥BN(x)∥∞≤CN!(2π)-N, where C is a positive constant independent of N. This completes the proof.Corollary 9.
Assume thatu(x,y)∈H×H=L2[0,1]×L2[0,1] be an enough smooth function and also is approximated by the two variable Bernoulli series ∑m=0∞‍∑n=0∞‍um,nBm(x)Bn(y), then the coefficients um,n for all m,n=0,1,…,∞ can be calculated from the following relation:
(45)um,n=1m!n!∫01‍∫01‍∂m+nu(x,y)∂xm∂yndxdy.Proof.
By applying a similar procedure in two variables (which was provided in Corollary6) 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 thatu(x,y) be an enough smooth function and uN(x,y) be the approximated polynomial of u(x,y) in terms of linear combination of Bernoulli polynomials by the aid of Corollary 9. Then the error bound would be obtained as follows:
(46)∥E(uN(x,y))∥∞:=∥u(x,y)-uN(x,y)∥∞≤CU^N(2π)-N,
where U^ is a positive constant independent of N and is a bound for all the partial derivatives of u(x,y).We now consider the basic equation (1) together with the initial conditions (2). For clarity of presentation we assume that m=1,P(z)=P0(z), and β=β0. A similar procedure would be applied in the case of larger values of m. Equation (1) with the above-mentioned assumptions has the following form:
(47)f′(z)+P(z)f(z)=G(z),
where
(48)f(z)=u(x,y)+iv(x,y),P(z)=p(x,y)+iq(x,y),G(z)=g(x,y)+ih(x,y).
Also the initial condition f(0)=u(0,0)+iv(0,0)=β=β1+iβ2 should be considered. Thus, one can write u(0,0)=β1 and v(0,0)=β2.According to Remark5, since f(z) is a Holomorphic function, then (47) by using assumptions in (48) could be rewritten as follows:
(49)ux(x,y)+ivx(x,y)+(p(x,y)+iq(x,y))(u(x,y)+iv(x,y))=g(x,y)+ih(x,y).
In other words
(50)ux(x,y)+(p(x,y)u(x,y)-q(x,y)v(x,y))=g(x,y),vx(x,y)+(q(x,y)u(x,y)+p(x,y)v(x,y))=h(x,y).
For imposing the initial conditions u(0,0)=β1 and v(0,0)=β2, we need to differentiate the above equations with respect to y and then integrate with respect to x and y in the rectangular [0,x]×[0,y]. Therefore, by differentiating both of the equations of (50) with respect to y we have
(51)uxy(x,y)+∂∂y(p(x,y)u(x,y)-q(x,y)v(x,y))=gy(x,y),vxy(xy)+∂∂y(q(x,y)u(x,y)+p(x,y)v(x,y))=hy(x,y).
Integrating from the above equations in the rectangular [0,x]×[0,y] yields
(52)∫0x‍∫0y‍uxy(x,y)dydx+∫0x‍∫0y‍∂∂y(p(x,y)u(x,y)-q(x,y)v(x,y))dydx=∫0x‍∫0y‍gy(x,y)dydx,∫0x‍∫0y‍vxy(x,y)dydx+∫0x‍∫0y‍∂∂y(q(x,y)u(x,y)+p(x,y)v(x,y))dydx=∫0x‍∫0y‍hy(x,y)dydx.
By imposing the the initial conditions u(0,0)=β1 and v(0,0)=β2 to these equations we reach to
(53)u(x,y)+∫0x[(p(x,y)u(x,y)-q(x,y)v(x,y))u(x,y)+∫0x-(p(x,0)u(x,0)-q(x,0)v(x,0))]dx=g^(x,y),v(x,y)+∫0x[(q(x,y)u(x,y)+p(x,y)v(x,y))u(x,y)+∫0x-(q(x,0)u(x,0)+p(x,0)v(x,0))]dx=h^(x,y),
where
(54)g^(x,y)=∫0x‍∫0y‍gy(x,y)dydx+β1,h^(x,y)=∫0x‍∫0y‍hy(x,y)dydx+β2.
Considering U(x,y)=[u(x,y)v(x,y)],G^(x,y)=[g^(x,y)h^(x,y)] and H(x,y)=[p(x,y)-q(x,y)q(x,y)p(x,y)], (53) could be restated in the following matrix vector form:
(55)U(x,y)+∫0x‍H(x,y)U(x,y)dx-∫0x‍H(x,0)U(x,0)dx=G^(x,y).
Our aim is to show that limN→∞UN(x,y)=U(x,y) or limN→∞fN(z)=f(z), where UN(x,y)=[uN(x,y)vN(x,y)],fN(z)=uN(x,y)+ivN(x,y) and fN(z) 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)
∥H(x,y)∥∞≤H^≪1,whereH^is a real constant,(ii)
limN→∞∥U(x,0)-UN(x,0)∥∞≤limN→∞(1/N)∥U(x,y)-UN(x,y)∥∞.It should be noted that the second condition is based upon Lemmas8 and 10.Theorem 11.
Assume thatU and U0:=U(x,0) be approximated by UN and U0,N, 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
(56)UN(x,y)+∫0x‍H(x,y)UN(x,y)dx-∫0xH(x,0)UN(x,0)dx=G^(x,y)+RN+1(x,y),
where RN+1(x,y) is the residual function that is zero at the (N+1) collocation nodes. Also suppose that [a,b]=[c,d]=[0,1]. Then, under the above-mentioned assumptions
(57)∥U-UN∥∞≤∥RN+1(x,y)∥∞1-H^(1+1/N),
and limN→∞UN(x,y)=U(x,y).Proof.
By subtracting (56) from (55) we have
(58)∥U-UN∥∞=∥∫0x(-HU+H0U0+HUN-H0U0,N)dx+RN+1(x,y)∫0x(-HU+H0U0+HUN-H0U0,N)dx∥∞≤∥∫0x(HU-HUN)dx∥∞+∥∫0x(H0U0-H0U0,N)dx∥∞+∥∫0x(H0U0-H0U0,N)dxRN+1(x,y)∥∞≤∫0x‍(H^∥U-UN∥∞+H^∥U0-U0,N∥∞)dx+∥RN+1(x,y)∥∞≤H^(∥U-UN∥∞(1+1N))+∥RN+1(x,y)∥∞.
Therefore,
(59)∥U-UN∥∞≤∥RN+1(x,y)∥∞1-H^(1+1/N).
SincelimN→∞∥RN+1(x,y)∥∞=0andH≪1, we have limN→∞∥U-UN∥∞=0; in other words, limN→∞fN(z)=f(z) 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 solutionf(z), the polynomial approximate solution fN(z), and the absolute error function eN(z)=|f(z)-fN(z)| 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(60)f′(z)+zf(z)=2z2-z+2,z=x+iy,f′(z)+zf(z)2z2-zx∈[0,1],y∈[0,1],
with the initial condition f(0)=-1 and the exact solution 2z-1. We suppose that N=3. Therefore, the collocation points are z00=0,z11=(1+i)/3,z22=(2+2i)/3, and z33=1+i. According to (3), the approximate solution has the following form:
(61)f3(z)=∑n=03‍fnBn(z),z∈D,
where our aim is to find the unknown Bernoulli coefficients f0,f1,f2,andf3. Since P0(z)=z, then
(62)P0=[000000.3333+0.3333i00000.6667+0.6667i00001.0000+1.0000i].
Also the matrix L (with the assumption N=3) has the following structure:(63)L=[B(z00)B(z11)B(z22)B(z33)]=[B0(z00)B1(z00)B2(z00)B3(z00)B0(z11)B1(z11)B2(z11)B3(z11)B0(z22)B1(z22)B2(z22)B3(z22)B0(z33)B1(z33)B2(z33)B3(z3)]=[1.0000-0.50000.166701.0000-0.1667+0.3333i-0.1667-0.1111i0.0926-0.0926i1.00000.1667+0.6667i-0.5000+0.2222i-0.2593-0.4074i1.00000.5000+1.0000i-0.8333+1.0000i-1.5000-0.5000i],where B0(z)=1,B1(z)=z-1/2,B2(z)=z2-z+1/6andB3(z)=z3-(3/2)z2+z/2.
According to (29), the matrix coefficients W are as follows:(64)W=L(MT)+P0L=[01.0000-1.00000.50000.3333+0.3333i0.8333+0.0556i-0.3519+0.5741i-0.4383-0.3333i0.6667+0.6667i0.6667+0.5556i-0.1481+1.1481i-1.4012+0.2222i1.0000+1.0000i0.5000+1.5000i-0.8333+2.1667i-3.5000+1.0000i].Since G(z)=2z2-z+2, the right-hand side vector G has the following form:
(65)G=[G(z00)G(z11)G(z22)G(z33)]=[2.00001.6667+0.1111i1.3333+1.1111i1.0000+3.0000i].
The associated form of the initial condition f(0)=-1 is
(66)[U0;β0]=[B(0);G(0)]=[1.0000-0.50000.16670;-1].
Imposing the above initial condition to the matrix W and vector G yields(67)W^=[01.0000-1.00000.50000.3333+0.3333i0.8333+0.0556i-0.3519+0.5741i-0.4383-0.3333i0.6667+0.6667i0.6667+0.5556i-0.1481+1.1481i-1.4012+0.2222i1.0000-0.50000.16670],G^=[2.00001.6667+0.1111i1.3333+1.1111i-1].Therefore, the solution of the system W^F=G^ is as follows:
(68)F=[0.0000-0.0000i2.0000-0.0000i0.0000-0.0000i0.0000+0.0000i].
Thus, we obtain the approximate solution f3(z)=B(z)F=2z-1 which is the exact solution. We recall that
(69)B(z)=⌊1z-12z2-z+16z3-32z2+z2⌋.Example 13 (see [26]).
As the second example, we consider the following second-order complex differential equation(70)f′′(z)+zf′(z)+zf(z)=exp(z)+2zexp(z),f′′(z)2exp(z),z=x+iy,x∈[0,1],y∈[0,1],
with the initial conditions f(0)=f′(0)=1 and the exact solution f(z)=exp(z). In this equation, we have P0(z)=P1(z)=z and G(z)=exp(z)+2zexp(z). Then, for N=5 the collocation points are z00=0,z11=(1+i)/5,z22=(2+2i)/5,z33=(3+3i)/5,z44=(4+4i)/5,andz55=1+i.
According to (29), the fundamental matrix equation is
(71)W=L(MT)2+P1LMT+P0L,
where MT and L are introduced in (18) and (25), respectively. Also
(72)P0=P1=[000000.3333+0.3333i00000.6667+0.6667i00001.0000+1.0000i].
The right-hand side vector G is
(73)G=[1.00001.5788+0.8185i2.0086+2.1449i2.0739+4.0681i1.4770+6.6318i-0.1686+9.7995i].
The augmented matrix forms of the initial conditions for N=5 are
(74)[B(0);1]=[1.0000-0.50000.16670-0.03330;1],[B(0)MT;1]=[01.0000-1.00000.50000-0.1667;1].
By replacing the above augmented vectors to the last two rows of [W;G], we reach to [W^;G^]. The solution of the matrix-vector equation W^F=G^ is
(75)F=[1.7184-0.0000i1.7189+0.0003i0.8603+0.0021i0.2864+0.0046i0.0693+0.0048i0.0107+0.0034i].
We also solve this equation by assumptions N=7 and N=9. Since the exact solution of the equation is exp(z)=exp(x+iy)=exp(x)(cos(y)+isin(y)). Therefore, the real and imaginary parts of the exact solution are Re(exp(z))=exp(x)cos(y) and Im(exp(z))=exp(x)sin(y), respectively. The values of the approximate solution in the case of N=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 N=7and9. From this figure, one can see the efficiency of our methods with respect to the method of [26].Table 1
Numerical results of Example13, for Re(exp(z))=exp(x)cos(y).
z
=
x
+
i
y
Exact solution (real)
PM forN=5
PM forN=7
PM forN=9
0.0
+
0.0
i
1.000000000000000
1.000000000000002
0.999999999999999
1.000000000000002
0.1
+
0.1
i
1.099649666829409
1.099650127059032
1.099649681790901
1.099649666809790
0.2
+
0.2
i
1.197056021355891
1.197058065414137
1.197056071303251
1.197056021305326
0.3
+
0.3
i
1.289569374044936
1.289572683703872
1.289569451705194
1.289569373971185
0.4
+
0.4
i
1.374061538887522
1.374064757066184
1.374061646628469
1.374061538795035
0.5
+
0.5
i
1.446889036584169
1.446891795476775
1.446889179088898
1.446889036484456
0.6
+
0.6
i
1.503859540558786
1.503862872087464
1.503859710928599
1.503859540462517
0.7
+
0.7
i
1.540203025431780
1.540203451564549
1.540203240474244
1.540203025363573
0.8
+
0.8
i
1.550549296807422
1.550520218427163
1.550549536105736
1.550549296761452
0.9
+
0.9
i
1.528913811884699
1.528765905385641
1.528913283429636
1.528913812578978
1.0
+
1.0
i
1.468693939915885
1.468204121679876
1.468688423751699
1.468693951050675Table 2
Numerical results of Example13, for Im(exp(z)=exp(z)sin(y).
z
=
x
+
i
y
Exact solution (Im)
PM forN=5
PM forN=7
PM forN=9
0.0
+
i
0.0
0
0.000000000000000
−0.000000000000001
0.000000000000000
0.1
+
i
0.1
0.110332988730204
0.110330942667805
0.110332993217700
0.110332988786893
0.2
+
i
0.2
0.242655268594923
0.242645839814175
0.242655282953655
0.242655268747208
0.3
+
i
0.3
0.398910553778490
0.398893389645474
0.398910574081883
0.398910554017485
0.4
+
i
0.4
0.580943900770567
0.580922057122650
0.580943925558894
0.580943901105221
0.5
+
i
0.5
0.790439083213615
0.790412124917929
0.790439110267023
0.790439083643388
0.6
+
i
0.6
1.028845666272092
1.028807744371502
1.028845687924972
1.028845666809275
0.7
+
i
0.7
1.297295111875269
1.297248986448220
1.297295129424920
1.297295112510736
0.8
+
i
0.8
1.596505340600251
1.596503892694281
1.596505330955563
1.596505341401562
0.9
+
i
0.9
1.926673303972717
1.926900526193919
1.926672856270450
1.926673303646987
1.0
+
i
1.0
2.287355287178843
2.288259022526098
2.287352245460983
2.287355267120311Figure 1
Comparisons of the presented method (PM) and the Taylor method (TM) of Example13.Example 14 (see [26]).
As the final example, we consider the following second-order complex differential equation(76)f′′(z)+zf′(z)+2zf(z)=2zsin(z)+zcos(z)-sin(z),z=x+iy,=2zsin(z)+zcos(z)c-,x∈[0,1],y∈[0,1],
with the initial conditions f(0)=0 and f′(0)=1 and also the exact solution f(z)=sin(z). In this equation, we have P0(z)=2z,P1(z)=z, and G(z)=2zsin(z)+zcos(z)-sin(z). Similar to the previous two examples, we solve this equation in the case of N=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 W^ of the Bessel method is very larger than the matrix coefficient W^ 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 Example14.Figure 3
Condition number comparisons between our method and the Bessel method of Example14.


## 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 functionsPk(z) and G(z) 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 N or less than N. 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.



---

*Source: 823098-2013-04-08.xml*