# A Numerical Method for Partial Differential Algebraic Equations Based on Differential Transform Method

Authors: Murat Osmanoglu; Mustafa Bayram
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We have considered linear partial differential algebraic equations (LPDAEs) of the formAut(t,x)+Buxx(t,x)+Cu(t,x)=f(t,x), which has at least one singular matrix of A,B∈ℝn×n. We have first introduced a uniform differential time index and a differential space index. The initial conditions and boundary conditions of the given system cannot be prescribed for all components of the solution vector u here. To overcome this, we introduced these indexes. Furthermore, differential transform method has been given to solve LPDAEs. We have applied this method to a test problem, and numerical solution of the problem has been compared with analytical solution.


---

## Body

## 1. Introduction

The partial differential algebraic equation was first studied by Marszalek. He also studied the analysis of the partial differential algebraic equations [1]. Lucht et al. [2–4] studied the numerical solution and indexes of the linear partial differential equations with constant coefficients. A study about characteristics analysis and differential index of the partial differential algebraic equations was given by Martinson and Barton [5, 6]. Debrabant and Strehmel investigated the convergence of Runge-Kutta method for linear partial differential algebraic equations [7].There are numerous LPDAEs applications in scientific areas given, for instance,   in the field of Navier-Stokes equations, in chemical engineering, in magnetohydrodynamics, and in the theory of elastic multibody systems [4, 8–12].On the other hand, the differential transform method was used by Zhou [13] to solve linear and nonlinear initial value problems in electric circuit analysis. Analysis of nonlinear circuits by using differential Taylor transform was given by Köksal and Herdem [14]. Using one-dimensional differential transform, Abdel-Halim Hassan [15] proposed a method to solve eigenvalue problems. The two-dimensional differential transform methods have been applied to the partial differential equations [16–19]. The differential transform method extended to solve differential-difference equations by Arikoglu and Ozkol [20]. Jang et al. have used differential transform method to solve initial value problems [21]. The numerical solution of the differential-algebraic equation systems has been studied by using differential transform method [22, 23].In this paper, we have considered linear partial differential equations with constant coefficients of the form(1)Aut(t,x)+Buxx(t,x)+Cu(t,x)=f(t,x),(t,x)∈J×Ω,
where J=[0,∞), Ω=[-l,l],l>0, and A,B,C∈ℝn×n.  In (1) at least one of the matrices            A,B∈ℝn×n should be singular.  If A=0 or B=0, then (1)  becomes ordinary differential equation or differential algebraic equation, so we  assume that none of the matrices A  or  B is the zero matrix.


## 2. Indexes of Partial Differential Algebraic Equation

Let us consider (1), with initial values and boundary conditions given as follows:
(2)uj(t,±l)=0fort∈J,ui(0,x)=g(x)forx∈Ω,
where j∈𝔐BC⊆{1,2,…,n},  𝔐BC is the set of indices of components of  u for which boundary conditions can be prescribed arbitrarily, and  i∈𝔐IC⊆{1,2,…,n},  𝔐IC is the set of indices of components of  u for which initial conditions can be prescribed arbitrarily. The initial boundary value problem (IBVP) (1) has only one solution where a function u is a solution of the problem, if it is sufficiently smooth, uniquely determined by its initial values (IVs) and boundary values (BVs), and if it solves the LPDAE point wise.Definition of the indexes can be given using the following assumptions.(i)
Each component of the vectorsu,ut, and f satisfy  the following condition:
(3)|y(t,x)|≤Meαt,α≥0,t≥0,
where M and α are independent  of t and x.(ii)
(B,ξA+C),Re(ξ)>α, called as the matrix pencil,  is regular.(iii)
(A,μkB+C) is regular for all k, where μk is an eigenvalue of the operator ∂2/∂x2 together with prescribed BCs.(iv)
The vectorf(t,x) and the initial vector g(x) are sufficiently smooth.If we use Laplace transform, from assumption (ii), (1) can be transformed into
(4)Buξ′′(x)+(ξA+C)uξ(x)=fξ(x)+Ag(x),Re(ξ)>α,
if B is a singular matrix, then (4) is a DAE depending on the parameter ξ. To characterize  𝔐BC, we introduce j∈𝔐BC(ξ)⊆{1,2,…,n} as the set of indices of components of  uξ for which boundary conditions can be prescribed arbitrarily.In order to define a spatial index, we need the Kronecker normal form of the DAE (4). Assumption (iii) guarantees that there are nonsingular matrices PL,ξ,QL,ξ∈ℂn×n such that
(5)PL,ξBQL,ξ=(Im100NL,ξ),PL,ξ(ξA+C)QL,ξ=(RL,ξ00Im2),
where  RL,ξ∈ℂm1×m2 and NL,ξ∈ℝm2×m2 is a nilpotent Jordan chain matrix with m1+m2=n.Ik is the unit matrix of order k. The Riesz index (or nilpotency) of NL,ξ is denoted by νL,ξ (i.e., NL,ξνL,ξ=0,NL,ξνL,ξ-1≠0).Here, we will assume that there is a real numberα*≥α such that the index set 𝔐BC(ξ) is independent of the Laplace parameter ξ, provided Re(ξ)≥α*.Definition 1.
Letα*∈ℝ+ be a number with α*≥α, such that for all ξ∈ℂ with Re(ξ)≥α*(1)
the matrix pencil(B,ξA+C) is regular,(2)
𝔐BC(ξ) is independent of ξ, i.e., 𝔐BC(ξ)=𝔐BC,(3)
the nilpotency ofNL,ξ is νL≥1.
Thenνd,x=2νL-1 is called the “differential spatial index” of the LPDAE. If νL=0, then the differential spatial index of LPDAE is defined to be zero.If we use Fourier transform, (1) can be transformed into
(6)Au^k′(t)+(μkB+C)u^k(t)=f^k(t)+Bρk(t)
with ρk(t)=(ρk1(t),…,ρkn(t))T and
(7)ρki(t)=0fori∈𝔐BC,ρkj(t)=1l[ϕk′(x)uj(t,x)-ϕk(x)ux,j(t,x)]x=-lx=l
for j∉𝔐BC, which results from partial integration of  the term ∫-lluxx(t,x)ϕk(x)dx.IfA is a singular matrix, then  (6) is a DAE depending on the parameter μk which can be solved uniquely with suitable ICs under the assumptions (iv) and (v). Analogous to the case of the Laplace transform, the above assumption (iv) implies that there exist regular matrices  PF,k,QF,k such that
(8)PF,kAQF,k=(In100NF,k),PF,k(μkB+C)QF,k=(RF,k00In2).
With RF,k∈ℝn1×n1. NF,k∈ℝn2×n2 is again a nilpotent Jordan chain matrix with Riesz index νF,k, where n1+n2=n.To characterize𝔐IC, we introduce 𝔐IC(k)⊆{1,2,…,n} as the set of indices of components of  u^k for which initial  conditions can be prescribed arbitrarily. Therefore, we always assume in the context of a Fourier analysis of u that 𝔐IC(k) is independent of k∈ℕ+, i.e.,  𝔐IC(k)=𝔐IC.Definition 2.
Assume fork=1,2,… that(1)
the matrix pencil(A,μkB+C) is regular,(2)
𝔐IC(k) is independent of k, i.e., 𝔐IC(k)=𝔐IC,(3)
the nilpotency ofNF,k is νF,k=νF.
Then the PDAE  (1)  is said to have uniform differential time indexνd,t=νF.The differential spatial and time indexes are used to decide which initial and boundary values can be taken to solve the problem.


## 3. Two-Dimensional Differential Transform Method

The two-dimensional differential transform of functionw(x,y) is defined as
(9)W(k,h)=1k!h![∂k+hw(x,y)∂xk∂yh]x=0y=0,
where it is noted that upper case symbol W(k,h) is used to denote the two-dimensional differential transform of a function represented by a corresponding lower case symbol w(x,y). The differential inverse transform of W(k,h) is defined as
(10)w(x,y)=∑k=0∞∑h=0∞W(k,h)xkyh.
From (9) and  (10),  we obtain
(11)w(x,y)=∑k=0∞∑h=0∞xkyhk!h![∂k+hw(x,y)∂xk∂yh]x=0y=0.
The concept of two-dimensional differential transform is derived from two-dimensional Taylor series expansion, but the method doesn’t evaluate the derivatives symbolically.Theorem 3.
Differential transform of the functionw(x,y)=u(x,y)±v(x,y) is
(12)W(k,h)=U(k,h)±V(k,h),
see [17].Theorem 4.
Differential transform of the functionw(x,y)=λu(x,y) is
(13)W(k,h)=λU(k,h),
see [17].Theorem 5.
Differential transform of the functionw(x,y)=∂u(x,y)/∂x is
(14)W(k,h)=(k+1)U(k+1,h),
see [17].Theorem 6.
Differential transform of the functionw(x,y)=∂u(x,y)/∂y is
(15)W(k,h)=(h+1)U(k,h+1),
see [17].Theorem 7.
Differential transform of the functionw(x,y)=∂r+su(x,y)/∂xr∂ys  is
(16)W(k,h)=(k+1)(k+2)⋯(k+r)(h+1)×(h+2)⋯(h+s)U(k+r,h+s),
see [17].Theorem 8.
Differential transform of the functionw(x,y)=u(x,y)·v(x,y) is
(17)W(k,h)=∑r=0k∑s=0hU(r,h-s)V(k-r,s),
see [17].Theorem 9.
Differential transform of the functionw(x,y)=xmyn is
(18)W(k,h)=δ(k-m,h-n)=δ(k-m)δ(h-n),
see [17], where
(19)δ(k-m)={1,k=m0,k≠m,δ(h-n)={1,h=n0,h≠n.Theorem 10.
Differential transform of the functionw(x,y)=g(x+a,y)  is
(20)W(k,h)=∑p=kN(pk)ap-kG(p,h).Proof.
From Definition1, we can write
(21)w(x,y)=∑h=0∞∑k=0∞G(k,h)(x+a)kyh=G(0,0)+G(0,1)y+aG(1,0)+G(1,0)x+G(0,2)y2+G(1,1)xy+aG(1,1)y+a2G(2,0)+2aG(2,0)x+G(2,0)x2+G(0,3)y3+G(1,2)xy2+aG(1,2)y2+G(2,1)x2y+2aG(2,1)xy+a2G(2,1)y+a3G(3,0)+3a2G(3,0)x+3aG(3,0)x2+G(3,0)x3+⋯,w(x,y)=[a3G(3,0)G(0,0)+aG(1,0)+a2G(2,0)+a3G(3,0)+⋯]+x[G(1,0)+2aG(2,0)+3a2G(3,0)+⋯]+y[G(0,1)+aG(1,1)+a2G(2,1)+⋯]+x2[G(2,0)+3aG(3,0)+⋯]+xy[G(1,1)+2aG(2,1)+⋯]+y2[G(0,2)+aG(1,2)+⋯]+⋯=∑p=0∞apG(p,0)+x∑p=1∞pap-1G(p,0)+y∑p=0∞apG(p,1)+x2∑p=2∞p!(p-2)!2!×ap-2G(p,0)+xy∑p=1∞pap-1G(p,1)+⋯,
where
(22)w(x,y)=∑h=0∞∑k=0∞∑p=k∞(pk)ap-kG(p,h)xkyh
hence,
(23)W(k,h)=∑p=kN(pk)ap-kG(p,h).Theorem 11.
Differential transform of the functionw(x,y)=g(x+a,y+b)  is
(24)W(k,h)=∑p=kN∑q=hN(qh)(pk)ap-kbq-hG(p,q).Proof.
From Definition2, we  can write
(25)w(x,y)=∑h=0∞∑k=0∞G(k,h)(x+a)k(y+b)h=G(0,0)+G(1,0)x+aG(1,0)+G(0,1)y+bG(0,1)+G(2,0)x2+2aG(2,0)x+a2G(2,0)+G(1,1)xy+bG(1,1)x+aG(1,1)y+abG(1,1)+G(0,2)y2+2bG(0,2)y+b2G(0,2)+G(3,0)x3+3aG(3,0)x2+3a2G(3,0)x+a3G(3,0)+G(2,1)x2y+⋯w(x,y)=[a2G(2,0)G(0,0)+aG(1,0)+bG(0,1)+a2G(2,0)+abG(1,1)+⋯]+[3a2G(3,0)G(1,0)+2aG(2,0)+bG(1,1)+3a2G(3,0)+2abG(2,1)+⋯]x+[a2G(2,1)G(0,1)+aG(1,1)+2bG(0,2)+a2G(2,1)+⋯]y+⋯=∑p=0N∑q=0NapbqG(p,q)+x∑p=1N∑q=0Npap-1bqG(p,q)+y∑p=0N∑q=1Nqapbq-1G(p,q)+x2∑p=2N∑q=0Np!(p-2)!2!ap-2bqG(p,q)+xy∑p=1N∑q=1Npqap-1bq-1G(p,q)+⋯.
Hence, we can write
(26)w(x,y)=∑h=0∞∑k=0∞∑p=kN∑q=hN(qh)(pk)×ap-kbq-hG(p,q)xkyh.
Using Definition 2, we obtain
(27)W(k,h)=∑p=kN∑q=hN(qh)(pk)ap-kbq-hG(p,q).Theorem 12.
Differential transform of the functionw(x,y)=∂r+sg(x+a,y+b)/∂xr∂ys  is
(28)W(k,h)=(k+r)!k!(h+s)!h!×∑p=k+rN∑q=h+sN(qh+s)(pk+r)×ap-k-rbq-h-sG(p,q).Proof.
LetC(k,h) be differential transform of the function g(x+a,y+b). From Theorem 7, we can write that differential transform of the function w(x,y) is
(29)W(k,h)=(k+r)!k!(h+s)!h!C(k+r,h+s),
from Theorem 4, we can write
(30)C(k+r,h+s)=∑p=k+rN∑q=h+sN(qh+s)(pk+r)×ap-k-rbq-h-sG(p,q).
If we substitute (30) into (29), we find
(31)W(k,h)=(k+r)!k!(h+s)!h!×∑p=k+rN∑q=h+sN(qh+s)(pk+r)×ap-k-rbq-h-sG(p,q).


## 4. Application

We have considered the following PDAE as a test problem:(32)(1100)ut+(-10-10)uxx+(11-10)u=f,t∈[0,∞),x∈[-1,1],
with initial values and boundary values
(33)u1(0,x)=x3-x,u2(0,x)=x4-1,u1(t,1)=u1(t,-1)=0,u2(t,1)=u2(t,-1)=0.
The right hand side function f is
(34)f=((x4-1)(cost-sint)-6xe-t,-e-t(x3+5x))T,
and the exact solutions are
(35)u1(t,x)=(x3-x)e-t,u2(t,x)=(x4-1)cost.
If nonsingular matrices PF,k, QF,k, PL,ξ, and QL,ξ are chosen such as
(36)PF,k=(1-μkμk-101μk-1),QF,k=(0-11-1),PL,ξ=(01ξ+11ξ+1-1ξ+1),QL,ξ=(-ξ-10ξ1),
matrices PF,kAQF,k and PL,ξBQL,ξ are found as
(37)PF,kAQF,k=(1000),PL,ξBQL,ξ=(1000).
From (38), we have NL,ξ=0 and NF,k=0. Then the PDAE (32) has differential spatial index 1 and differential time index 1. So, it is enough to take 𝔐BC(ξ)={1} and 𝔐IC(k)={2}  to solve the problem.Taking differential transformation of (32), we obtain
(38)(k+1)U1(k+1,h)+(k+1)U2(k+1,h)-(h+1)(h+2)U1(k,h+2)+U1(k,h)+U2(k,h)=F1(k,h),(39)-(h+1)(h+2)U1(k,h+2)-U1(k,h)=F2(k,h).
The Taylor series of functions f1 and f2 about x=0,t=0 are
(40)f1(t,x)=-1+t-6x+12t2+6xt-16t3-3xt2-124t4+x4+xt3-x4t+1120t5-14xt4+1720t6-12x4t2+120xt5+16x4t3-15040t7-1120xt6+⋯,(41)f2(t,x)=-5x+5xt-52xt2-x3+56xt3+x3t-12x3t2-524xt4+16x3t3+124xt5-1144xt6-124x3t4+⋯.
The values F1(k,h)andF2(k,h) in (39) and (40) are coefficients of polynomials (41) and (42).   If we use Theorem 3 for boundary values, we obtain
(42)∑i=07U1(j,i)=0,j=0,1,…,7,(43)∑i=07(-1)iU1(j,i)=0,j=0,1,…,7.
In order to write k=0 and h=0,1,2,3,4,5   in (40), we have
(44)2U1(0,2)+U1(0,0)=0,20U1(0,5)+U1(0,3)=1,6U1(0,3)+U1(0,1)=5,30U1(0,6)+U1(0,4)=0,12U1(0,4)+U1(0,2)=0,42U1(0,7)+U1(0,5)=0.
If we take j=0 in (43) and (44), we obtain
(45)U1(0,0)+U1(0,1)+U1(0,2)+U1(0,3)+U1(0,4)+U1(0,5)+U1(0,6)+U1(0,7)=0,U1(0,0)-U1(0,1)+U1(0,2)-U1(0,3)+U1(0,4)-U1(0,5)+U1(0,6)-U1(0,7)=0.
From (45) and (46), we find
(46)U1(0,0)=0,U1(0,1)=-1,U1(0,2)=0,U1(0,3)=1,U1(0,4)=0,U1(0,5)=0,U1(0,6)=0,U1(0,7)=0.
In this manner, from (40), (44), and (45), the coefficients of the u1 are obtained as follows:
(47)U1(1,0)=0,U1(1,1)=1,U1(1,2)=0,U1(1,3)=-1,U1(1,4)=0,U1(1,5)=0,U1(1,6)=0,U1(2,0)=0,U1(2,1)=-12,U1(2,2)=0,U1(2,3)=12,U1(2,4)=0,U1(2,5)=0,U1(3,0)=0,U1(3,1)=16,U1(3,2)=0,U1(3,3)=-16,U1(3,4)=0,U1(4,0)=0,U1(4,1)=-124,U1(4,2)=0,U1(4,3)=124,U1(5,0)=0,U1(5,1)=1120,U1(5,2)=0,U1(6,0)=0,U1(6,1)=-1720,U1(7,0)=0.
Using the initial values for the second component, we obtain the following coefficients:
(48)U2(0,0)=-1,U2(0,1)=0,U2(0,2)=0,U2(0,3)=0,U2(0,4)=1,U2(0,5)=0,U2(0,6)=0,U2(0,7)=0.
The coefficients of the u2 can be found using (47), (48), (49), and taking  k=0,1,2,… and h=0,1,2,… in (39) as follows:
(49)U2(1,2)=0,U2(1,3)=0,U2(1,4)=0,U2(1,5)=0,U2(1,6)=0,U2(2,1)=0,U2(2,2)=0,U2(2,3)=0,U2(2,4)=-12,U2(2,5)=0,U2(3,0)=0,U2(3,1)=0,U2(3,2)=0,U2(3,3)=0,U2(3,4)=0,U2(4,1)=0,U2(4,2)=0,U2(4,3)=0,U2(4,0)=-124,U2(5,0)=0,U2(5,1)=0,U2(5,2)=0,U2(6,0)=1720,U2(6,1)=0,U2(7,0)=0.
If we write the above values in (39) and (40), then we have
(50)u1*(t,x)=-x+xt-12xt2+x3+16xt3-x3t+12x3t2-124xt4-16x3t3+1120xt5-1720xt6+124x3t4+⋯,(51)u2*(t,x)=-1+12t2-124t4+x4+1720t6-12x4t2+⋯.
Numerical and exact solution of the given problem has been compared in Tables 1 and 2, and simulations of solutions have been depicted in Figures 1, 2, 3, and 4, respectively.Table 1
The numerical and exact solution of the test problem (32), where  u1(t,x) is the exact solution and u1*(t,x) is the numerical solution, for x=0.1.
t
u
1
(
t
,
x
)
u
1
*
(
t
,
x
)
|
u
1
(
t
,
x
)
-
u
1
*
(
t
,
x
)
|
0.1
−0.0895789043
−0.0895789043
0
0.2
−0.0810543445
−0.0810543422
0.0000000023
0.3
−0.0733410038
−0.0733409887
0.0000000151
0.4
−0.0663616845
−0.0663616355
0.0000000490
0.5
−0.0600465353
−0.0600464409
0.0000000944
0.6
−0.0543323519
−0.0543322800
0.0000000719
0.7
−0.0491619450
−0.0491621943
0.0000002493
0.8
−0.0444835674
−0.0444849422
0.0000013748
0.9
−0.0402503963
−0.0402546487
0.0000042524
1.0
−0.0364200646
−0.0364305555
0.0000104909Table 2
The numerical and exact solution of the test problem (32), where  u2(t,x) is exact solution and u2*(t,x) is numerical solution, for x=0.1.
t
u
2
(
t
,
x
)
u
2
*
(
t
,
x
)
|
u
2
(
t
,
x
)
-
u
2
*
(
t
,
x
)
|
0.1
−0.9949046649
−0.9949046653
0.0000000004
0.2
−0.9799685711
−0.9799685778
0.0000000067
0.3
−0.9552409555
−0.9552409875
0.0000000320
0.4
−0.9209688879
−0.9209689778
0.0000000899
0.5
−0.8774948036
−0.8774949653
0.0000001637
0.6
−0.8252530813
−0.8252532000
0.0000001187
0.7
−0.7647657031
−0.7647652653
0.0000004378
0.8
−0.6966370386
−0.6966345778
0.0000024608
0.9
−0.6215478073
−0.6215398875
0.0000079198
1.0
−0.5402482757
−0.5402277778
0.0000204979Figure 1
The graphic of the functionu1(t,x) in the test problem (32).Figure 2
The graphic of the  functionu1*(t,x) in the test problem (32).Figure 3
The graphic of the functionu2(t,x) in the test problem (32).Figure 4
The graphic of the functionu2*(t,x) in the test problem (32).


## 5. Conclusion

The computations associated with the example discussed above were performed by using Computer Algebra Techniques [24]. We show the results in Tables 1 and 2 for the solution of (32) by numerical method. The numerical values on Tables 1 and 2 obtained above are in full agreement with the exact solutions of (32). This study has shown that the differential transform method often shows superior performance over series approximants, providing a promising tool for using in applied fields.



---

*Source: 705313-2013-03-31.xml*