# Convergence of One-Leg Hybrid Methods for Implicit Mixed Differential Algebraic Systems

Authors: Iman H. Ibrahim; Fatma M. Yousry
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

This paper focuses on a hybrid multistep and its twin one-leg methods and implementing them on implicit mixed differential algebraic equations. The
orders of convergence for the above methods are discussed and numerical
tests are solved.


---

## Body

## 1. Introduction

Consider the ordinary differential system:(1)y′=ft,y,where the linear multistep method (LM) [1, 2] (2)∑i=0kαiym+i=h∑i=0kβiftm+i,ym+iis given, and the generating polynomials(3)ρζ=∑i=0kαiζi,σζ=∑i=0kβiζihave real coefficients and no common divisor. Also assume throughout the normalization that (4)σ1=1.Then the associated one-leg (OL) method is defined by (5)∑i=0kαiym+i=hf∑i=0kβitm+i,∑i=0kβiym+i.The author presents hybrid multistep methods that take the form(6a)yn+s=hμfn+∑j=0k-2γn-jyn-j,(6b)yn+∑j=1kαn-jyn-j=hβsfn+s-β∗fn-1,and the one-leg twin of (6a) and (6b) takes the form(7)yn+∑j=1kαn-jyn-j=hfβstn+s-βsβ∗tn-1,βsyn+s-βsβ∗yn-1,where fn+s=f(tn+s,yn+s), tn+s=tn+sh, -1<s<1, -1⩽β∗<1. βs, and αn-j, j=1,2,…,k, are parameters to be determined as functions of s and β∗. Methods (6a) and (6b) with step k have order p=k and yn+s has order k-1. To evaluate the value of yn+s at the off-step point, that is, tn+s, consider the nodes tn (double node), tn-1,…, and tn-k (simple nodes) [3, 4].Applying Newton’s interpolation formula for this data gives the following scheme:(8)ytn+sh=yn+shfn+s2hfn-∇yn+s2s+12!hfn-∇yn-12∇2yn+s2s+1s+23!hfn-∇yn-12∇2yn-∇3yn3+⋯,where ∇yn=yn-yn-1 and f (or f(t,y)) is considered as a derivative of the solution y(t).The hybrid multistep method and its twin one-leg depend on two parameters,β∗ and s, which control their convergence and stability; also the position of the stage point affects the stability regions of the methods. For optimal values of β∗ and s, the methods have larger stability region compared to the hybrid backward differentiation formulae [3]; the corresponding one-leg twin is G-stable for k=2 and k=3; see [5, 6].Differential-algebraic equations (DAEs) often take place in highly scientific technology domains, such as automatic control engineering, simulation of electrical networks, and chemical reaction kinetics [7, 8]. Some systems can be reduced to ODE systems and can be solved by numerical ODE methods. Reduction to explicit differential system (1) in some other systems can be impossible or impractical because the problem is more naturally posed in the form F(t,y′,y)=0 and a reduction might reduce the sparseness of Jacobian matrices. These systems are then solved directly [9, 10].Here LM (2) and OL (5) are defined for implicit mixed differential algebraic systems of the form(9a)Ft,y′,y,x=0,(9b)Gt,y,x=0,where F, G, x, and y are vectors of the same dimension. Rewrite (2) and (5) in the form(10)1hρyn-σyn′=0,(11)1hρyn-fσyn,τ=0,respectively, where τ=σtn, substituting for yn′ in (9a) and (9b):(12a)Ftn,1βk1hρyn-σ-yn′,yn,xn=0,(12b)Gtn,xn,yn=0,where σ-yn′ acts only on backward data. Equations (12a) and (12b) can be solved for (yn,xn). In the one-leg form, the arguments are changed to τn,y′(τn),y(τn), and x(τn). The implementation of OL (5) to (9a) and (9b) gives the equations (13a)Fτ,1hρyn,σyn,σxn=0,(13b)Gtn,yn,xn=0.As a modification technique that applies the same arguments ofF on G, this implementation can be written as(14a)Fτ,1hρyn,σyn,σxn=0,(14b)Gτ,σyn,σxn=0.The LMS and MOL formulations in (12a), (12b), (14a), and (14b) are easier to implement than OL method in (13a) and (13b) because both equations are evaluated on the same arguments.In the following section, the hybrid multistep (HMS) method in (6a) and (6b) and its twin, hybrid one-leg (7) (HOL) method are defined for (9a) and (9b) and expressions for the local truncation errors of (HMS) and (HOL) are given.


## 2. The Hybrid Method

In the case ofk=2, the method in (6a) and (6b) takes the form(15a)αnyn+αn-1yn-1+αn-2yn-2=hβsfn+s-β∗fn-1,(15b)yn+s=yn+shfn,where(16)αn=3+2s-β∗21-β∗,αn-1=-21+s1-β∗,αn-2=--1-2s-β∗21-β∗,βs=11-β∗.Method (15a) has order 2 and its truncation error is 2+3s(2+s)+β∗/6(-1+β∗)h3y′′′(η), and yn+s has order one and its truncation error is T2=s2h2y′′(η) due to (8), where η∈(tn-2,tn+1). Applying the method in (15a) and (15b) on implicit mixed differential algebraic equations (9a) and (9b) obtains the following:(17a)Ftn+s,un+s,yn+s,xn+s=0,(17b)Gtn,xn,yn=0,where(18)un+s=1hβsαnyn+αn-1yn-1+αn-2yn-2+β∗fn-1.Letx(t) and y(t) be the exact solution of (9a) and (9b). The residues of (17a) and (17b) are the values of the left sides evaluated on x(t) and y(t). Using Taytor expansion for the second argument of F evaluated on y(t),(19)utn+s=1hβsαnytn+αn-1ytn-1+αn-2ytn-2+β∗ftn-1,leads to (20)utn+s≃y′tn+s+T1,where (21)T1=16-2-6s-β∗h2y′′′η.Theorem 1.
The order of convergence of the second-order hybrid method in (15a) and (15b) when applied to implicit mixed DAEs (9a) and (9b) is two.Proof.
Let the local truncation errors be defined byϵn=yn-y(tn), ϕn=xn-x(tn), ϵ~n=yn′-y′(tn), and ϕ~n=xn′-x′(tn), where {xn+s,yn+s} satisfies (17a) and (17b) with exact backward data:(22)0=Ftn+s,yn+s′,yn+s,xn+s≃Ftn+s,1hβsαnyn+αn-1yn-1+αn-2yn-2+β∗fn-1,yn+shyn′,xn+shxn′≃Ftn+s,1hβsαnytn+ϵn+αn-1yn-1+αn-2yn-2+β∗fn-1,ytn+ϵn+shy′tn+ϵ~n,xtn+ϕn+shx′tn+ϕ~n≃Ftn+s,T1+ϵnαnhβs+y′tn+s,T2+ϵn+shϵ~n+ytn+s,T3+ϕn+shϕ~n+xtn+s,where T3=s2h2x′′η; expanding around (tn+s,y′(tn+s),y(tn+s),x(tn+s)) implies that(23)Ftn+s,y′tn+s,ytn+s,xtn+s+T1+ϵnαnhβsFy′+T2+ϵn+shϵ~nFy+T3+ϕn+shϕ~nFx≃0,where the arguments of F, Fy′, Fy, and Fx are (tn+s,y′(tn+s),y(tn+s),x(tn+s))(24)ϵnαnhβsFy′+Fy+T1Fy′+T2Fy+shϵ~nFy+T3+ϕn+shϕ~nFx≃0.The errors ϵn and ϕn satisfy the following equations:(25)ϵnαnhβsFy′+Fy+ϕnFx≃-T1Fy′+T2Fy+shϵ~nFy+T3+shϕ~nFx,(26)0=Gtn,yn,xn≃G+Gyϵn+Gxϕn,where G=F=0 and the arguments of G,Gy,Gx are (tn,y(tn),x(tn)). If Gx-1 and Fy′-1 exist, then(27)ϕn≃-Gx-1Gyϵn.Therefore, ϕn has the same order as ϵn. The substitution for ϕn in (25) implies that(28)ϵnI-hβsαnGx-1GyFx+hβsαnFy=γ1,where(29)γ1=-hβsαn·Fy′-1T1Fy′+T2Fy+shϵ~nFy+shϕ~Fx+T3Fx=Oh3,ϵnI+Oh=Oh3,ϵn=Oh3.Then ϵn and ϕn are third order; thus the method is of second order.
Therefore,ϵn is third order small and in accordance with classical theory we conclude that the global error ϵ-n in yn must be second order small. Furthermore, if ϕ-n denotes the global error in xn—note that the global errors (ϵ-n,ϕ-n) satisfy the same algebraic constraints as the local errors (ϵn,ϕn), namely, (26)—consequently, ϕ-n is also second order small and thus method in (15a) and (15b) is second order small and accurate with respect to both y and x.
Modified Technique for Hybrid Method. It is noticed that the arguments of (17a) are (tn+s,yn+s′,yn+s,xn+s) and that of (17b) is (tn,yn,xn). The arguments of (17a) and (17b) can be taken as(30)Ftn+s,yn+s′,yn+s,xn+s=0,Gtn+s,yn+s,xn+s=0,which is called the modified technique for hybrid method.
In this case,(31)0=Gtn+s,yn+s,xn+s≃Gtn+s,yn+shyn′,xn+shxn′≃Gtn+s,T2+ϵn+shϵ~n+ytn+s,T3+ϕn+shϕ~n+xtn+s.Expanding G around (tn+s,y(tn+s),x(tn+s)) gives the following:(32)G+T2+ϵn+shϵ~nGy+T3+ϕn+shϕ~nGx≃0,T3+ϕn+shϕ~n≃-T2+ϵn+shϵ~nGx-1Gy.Substitute for (T3+ϕn+shϕ~n) in (24):(33)ϵnI+hβsαnFy′-1Fy-Gx-1GyFx≃-hβsαnFy′-1T1Fy′+Fy-Gx-1GyFxT2+shϵ~,where(34)γ2=-hβsαnFy′-1T1Fy′+Fy-Gx-1GyFxT2+shϵ~=Oh3,ϵnI+Oh≃Oh3,ϵn=Oh3,and the global error ϵ-n affecting yn is O(h2).
Since(35)0=Gtn+s,yn+s,xn+s≃Gtn+s,ytn+ϵ-n+shy′tn+ϵ-~,xtn+ϕ-n+shx′tn+ϕ-~≃ytn+ϵ-n+shy′tn+ϵ-~-ytn+sGy+xtn+ϕ-n+shx′tn+ϕ-~-xtn+sGx≃ϵ-n+shϵ-~+T2Gy+ϕ-n+shϕ-~+T3Gx,ϕ-n≃-ϵ-n+shϵ-~+T2Gx-1Gy-shϕ-~+T3=Oh2,thus the method is of second order.


## 3. The One-Leg Twin

In the case ofk=2, method (7) takes the form(36a)αnyn+αn-1yn-1+αn-2yn-2=hfβstn+s-βsβ∗tn-1,βsyn+s-βsβ∗yn-1,(36b)yn+s=yn+shfn.Method (36a) has order 2 and its truncation error is 1/6-(1+s)2/2(-1+β∗)2h3y′′′(η) and yn+s has order one and its truncation error is T2=s2h2y′′(η) due to (8).Theorem 2.
The order of convergence of the second-order one-leg twin in (36a) and (36b) when applied to implicit mixed DAE is two.Proof.
Applying the method in (36a) and (36b) on the implicit mixed differential algebraic equations (9a) and (9b) obtains the following:(37a)Fτn,1hρyn,σyn,σxn=0,(37b)Gtn,yn,xn=0.
The residues of (37a) and (37b) are the values of the left sides evaluated on x(t) and y(t):(38)0≃Fτn,1hαnyn+αn-1yn-1+αn-2yn-2,βsyn+s-βsβ∗yn-1,βsxn+s-βsβ∗xn-1,0≃Fτn,1hαnytn+ϵn+αn-1yn-1+αn-2yn-2,βsytn+ϵn+shy′tn+ϵ~n-βsβ∗yn-1,βsxtn+ϕn+shx′tn+ϕ~n-βsβ∗xn-1,0≃Fτn,T-1+αnϵnh+y′τn,T4+βsϵn+βsshϵ~n+yτn,T5+βsϕn+βsshϕ~n+xτn,where(39)T-1=16-1+s22-1+β∗2h2y′′′η,T4=σytn-yτn=1+s2β∗2-1+β∗2h2y′′η,T5=σxtn-xτn=1+s2β∗2-1+β∗2h2x′′η,and expanding around τn,y′(τn),y(τn),x(τn) implies that(40)0≃Fτn,y′τn,yτn,xτn+T-1+αnϵnhFy′+T4+βsϵn+βsshϵ~nFy+T5+βsϕn+βsshϕ~nFx,(41)0≃T-1+αnϵnhFy′+T4+βsϵn+βsshϵ~nFy+T5+βsϕn+βsshϕ~nFx,where the arguments of F,Fy′,Fy,Fx are (τn,y′(τn),y(τn),x(τn)).
Expanding (37b) around (tn,y(tn),x(tn)), if Gx-1 exists, gives the following:(42)ϕn≃-Gx-1Gyϵn;substitute ϕn in (41):(43)ϵnI+hαnFy′-1βsFy-βsGx-1GyFx≃-hαn·Fy′-1T4+βsshϵ~nFy+T5+βsshϕ~nFx+T-1Fy′,ϵn1+Oh≃γ3=Oh3,where (44)γ3=-hαnFy′-1T4+βsshϵ~nFy+T5+βsshϕ~nFx+T-1Fy′,ϵn=Oh3;thus the method is of second order.
Modified Technique for One-Leg Twin Method. Here the arguments of F and G are different if the arguments of G are taken as (τn,y(τn),x(τn)), and (37a) and (37b) become the following:(45)Fτn,1hρyn,σyn,σxn=0,Gτn,σyn,σxn=0,Gτn,βsyn+s-βsβ∗yn-1,βsxn+s-βsβ∗xn-1=0,0=Gτn,T4+βsϵn+βsshϵ~n+yτn,T5+βsϕn+βsshϕ~n+xτn,0≃Gτn,yτn,xτn+T4+βsϵn+βsshϵ~nGy+T5+βsϕn+βsshϕ~nGx,T5+βsϕn+βsshϕ~n≃-T4+βsϵn+βsshϵ~n·Gx-1Gy.
Substitute(T5+βsϕn+βsshϕ~n) in (41):(46)0≃T-1+αnϵnhFy′+T4+βsϵn+βsshϵ~nFy-T4+βsϵn+βsshϵ~nGx-1GyFx,ϵnI+hαnFy′-1βsFy-βsGx-1GyFx≃γ4=Oh3,where(47)γ4=-hαnFy′-1T4+βsshϵ~nFy-T4+βsshϵ~nGx-1GyFx+T-1Fy′,ϵnI+Oh≃γ4=Oh3.Consequently, ϵn=O(h3) and the global error ϵ-n affecting yn is O(h2). However, since (48)0=Gτn,σyn,σxn=Gτn,σytn+σϵ-n,σxtn+σϕ-n≃σytn+σϵ-n-yτnGy+σxtn+σϕ-n-xτnGx≃σϵ-n+T4Gy+σϕ-n+T5Gx,the global error ϕ-n is related to ϵ-n by the difference (49)σϕ-n≃σϵ-n+T4Gx-1Gy-T5=Oh2.
The solutionϕ-n of this difference equation is also O(h2) since σ-1 is a bounded operator; thus the method is of second order.


## 4. Numerical Tests

Here, some numerical results are presented to evaluate the performance of the proposed technique [11, 12].Test 1. Consider the differential algebraic equations:(50)x′t=21-ysin⁡y+x1-y,0=x2+y-1cos2⁡y,with the initial conditions x(1)=1, y1=0, and the exact solution is xt=tcos⁡(1-t2), y(t)=1-t2.Test 2. Consider the nonlinear DAEs:(51)y1′-y2-y1y3=0,y2′+π32y1-y2y3=0,-π32y12-y22+1=0,t∈0,1,with the initial conditions y1(0)=0, y2(0)=1, and y3(0)=0.The exact solutions arey1(t)=3/πsin⁡((π/3)t), y2(t)=cos⁡((π/3)t), and y3(t)=0.Test 3. Consider the nonlinear DAEs:(52)x′t=fx,t-Bx,ty,0=gx,y,with(53)f=-x1+x2-sin⁡t-1+2t0,B=0x1,x=x1x2,g=x12+x1x2-sin⁡t-1+2t,subject to the initial conditions x1(0)=1, x2(0)=0, and y(t)=-1.The exact solution isxe=(1-2t,sin(t)),  ye=-cos⁡(t)/(1-2t).Test 4 (practical test). Consider rectifier diode circuit [13] in Figure 1 for transforming an AC voltage source into a DC voltage. It is designed in such a way that it damps the incoming sine-wave.Figure 1
Rectifier circuit.In Figure1, the circuit is a half wave rectifier circuit. The diode permits the flow of the current during the positive half cycle and stops the current flow during the negative half cycle. The capacitor is used to smooth the output voltage such that the output voltage at the load resistor R is close to a DC voltage:(54)vn1t=10sin⁡377t.Kirchhoff’s Current Law (KCL) at n1(55)ist=iDt,iDt=Isexp⁡38.685vn1t-vn2t-1.Assuming Is=10-13A,(56)ist=Isexp⁡38.685vn1t-vn2t-1.KCL at n2:(57)iDt=iRt+ict,Isexp⁡38.685vn1t-vn2t-1=vn2tR+cddtvn2t,where is is the supply current, iD is the diode current, iR is the current in the resistor R, iC is the current in the capacitor C, Is is the reverse current of the diode, n1 and n2 are node 1 and node 2, n0 is the reference node, and v is the supply voltage. We solve this circuit with hybrid formula (6a) and (6b) method and draw the exact and the numerical solutions in Figures 2–4.Figure 2
The exact and approximate solutions for the current of the diode circuit.Figure 3
The exact and numerical solutions for the voltage(vn2(t)) of the diode circuit.Figure 4
The exact and approximate solutions of the voltage(vn1(t)) of the diode circuit.The above tests are solved by the hybrid multistep method in (6a) and (6b) and its hybrid twin one-leg method (7) (with k=2,  β∗=-0.4, and s=-0.4) and compared with HBDF method [3] (with k=2 and s=0.3) at different values of t. In solving the tests, Newton-Raphson method is used and the initial guesses are obtained by an interpolating polynomial. The errors of solutions of tests 1, 2, and 3 are tabulated in Tables 1, 2, and 3.Table 1
The errors of solutions for the first test.
t
h
Er(x(t))
Er(y(t))
HMS
2
0.01
4.83216
E
-
4
9.09651
E
-
4
4
0.01
1.19431
E
-
2
4.42402
E
-
3
2
0.001
1.96883
E
-
6
3.71621
E
-
6
4
0.001
5.36086
E
-
4
1.98856
E
-
4
2
0.000l
1.68493
E
-
7
3.18069
E
-
7
4
0.0001
3.80065
E
-
7
1.40968
E
-
7
HOL
2
0.01
2.69539
E
-
4
5.05238
E
-
4
4
0.01
1.3771
E
-
2
5.10544
E
-
3
2
0.001
4.43816
E
-
7
8.37372
E
-
7
4
0.001
1.40624
E
-
4
5.21576
E
-
5
2
0.0001
3.55151
E
-
8
6.70428
E
-
8
4
0.0001
1.39345
E
-
6
5.16839
E
-
7
HBDF
2
0.01
1.99174
E
-
2
3.53307
E
-
2
4
0.01
2.83319
E
-
1
9.99465
E
-
2
2
0.001
6.02345
E
-
5
1.13683
E
-
4
4
0.001
1.37187
E
-
3
5.08695
E
-
4
2
0.0001
6.21028
E
-
8
1.17233
E
-
7
4
0.0001
4.78321
E
-
6
1.77411
E
-
6Table 2
The errors of solutions for the second test.
t
h
Er(y1(t))
Er(y2(t))
Er(y3(t))
HMS
0.5
0.01
4.44120
E
-
6
2.67822
E
-
6
5.5219
E
-
8
1
0.01
5.17837
E
-
6
9.38061
E
-
6
5.5219
E
-
8
0.5
0.001
4.51402
E
-
8
2.72911
E
-
8
5.54467
E
-
11
1
0.001
5.21721
E
-
8
9.46285
E
-
8
5.52825
E
-
11
0.5
0.0001
4.52156
E
-
10
2.73373
E
-
10
2.00697
E
-
12
1
0.0001
5.22150
E
-
10
9.47077
E
-
10
1.29584
E
-
12
HOL
0.5
0.01
1.10985
E
-
5
6.7033
E
-
6
2.27002
E
-
7
1
0.01
1.29380
E
-
5
2.34557
E
-
5
2.27002
E
-
7
0.5
0.001
1.1285
E
-
7
6.82285
E
-
8
2.27173
E
-
10
1
0.001
1.3043
E
-
7
2.36571
E
-
7
2.27181
E
-
10
0.5
0.0001
1.13038
E
-
9
6.8343
E
-
10
2.18404
E
-
12
1
0.0001
1.30538
E
-
9
2.36771
E
-
9
1.15477
E
-
12
HBDF
0.5
0.01
2.89796
E
-
5
1.75273
E
-
5
9.67763
E
-
7
1
0.01
3.40875
E
-
5
6.18349
E
-
5
9.67763
E
-
7
0.5
0.001
2.99655
E
-
7
1.81172
E
-
7
9.67982
E
-
10
1
0.001
3.46637
E
-
7
6.28731
E
-
7
9.67637
E
-
10
0.5
0.0001
3.00632
E
-
9
1.81762
E
-
9
8.01653
E
-
13
1
0.0001
3.47202
E
-
9
6.29754
E
-
9
5.0231
E
-
13Table 3
The errors of solutions for the third test.
t
h
Er(x1(t))
Er(x2(t))
Er(y(t))
HMS
1
0.01
1.0116
E
-
5
1.011656
E
-
5
5.88302
E
-
5
3
0.01
1.58952
E
-
6
1.58952
E
-
6
2.0603
E
-
6
1
0.001
1.03830
E
-
7
1.03831
E
-
7
6.02688
E
-
7
3
0.001
1.61765
E
-
8
1.61766
E
-
8
2.07402
E
-
8
1
0.0001
1.04072
E
-
9
1.04072
E
-
9
6.0429
E
-
9
3
0.0001
1.64812
E
-
10
1.64812
E
-
10
2.12782
E
-
10
HOL
1
0.01
4.47718
E
-
6
4.47718
E
-
6
1.0808
E
-
4
3
0.01
4.05837
E
-
6
4.05837
E
-
6
5.02559
E
-
6
1
0.001
5.35857
E
-
8
5.35858
E
-
8
8.50777
E
-
8
3
0.001
3.97665
E
-
8
3.97665
E
-
8
5.01833
E
-
8
1
0.0001
9.21768
E
-
10
9.21768
E
-
10
9.21768
E
-
10
3
0.0001
3.94626
E
-
10
3.94624
E
-
10
5.00513
E
-
10
HBDF
1
0.01
5.77901
E
-
7
5.77901
E
-
7
676543
E
-
5
3
0.01
1.07426
E
-
6
1.07427
E
-
6
1.2481
E
-
5
1
0.001
2.5353
E
-
8
2.5353
E
-
8
7.15841
E
-
7
3
0.001
1.1663
E
-
8
1.1663
E
-
8
1.24319
E
-
7
1
0.0001
2.71828
E
-
10
2.71828
E
-
10
7.19898
E
-
9
3
0.0001
1.11241
E
-
10
1.11242
E
-
10
1.23669
E
-
9


## 5. Conclusion

This paper focuses on the implementation of hybrid multistep classes and their twin one-leg classes on implicit mixed differential algebraic equations. The orders of convergence for these classes are discussed. Numerical tests are introduced, which show that the introduced methods give better results than HBDF.



---

*Source: 609015-2015-06-11.xml*