View of On the convergence of the variational iteration method

(1)

ON THE CONVERGENCE

OF THE VARIATIONAL ITERATION METHOD

ERNEST SCHEIBER^∗

Abstract. Convergence results are stated for the variational iteration method applied to solve an initial value problem for a system of ordinary differential equations.

MSC 2010. 65L20, 65L05, 97N80.

Keywords. variational iteration method, ordinary differential equations, initial value problem.

1. INTRODUCTION

The Ji-Huan He’s Variational Iteration Method (VIM) was applied to a large range of problems for both ordinary and partial differential equations.

The main ingredient of the VIM is the Lagrange multiplier used to improve an approximation of the solution of the differential problem [2].

The purpose of this paper is to prove a convergence theorem for VIM applied to solve an initial value problem for a system of ordinary differential equations.

The convergence of the VIM for the initial value problem of an ordinary differential equation may be found in D.K. Salkuyeh, A. Tavakoli [6]. For a system of linear differential equations a convergence result is given by D.K.

Salkuyeh [5].

A particularity of the VIM is that it may be implemented both in symbolic (Computer Algebra System) and numerical programming environments. In the last section there are presented some results of our computational expe- riences. To make the results reproducible we provide some code. In [1] there is a pertinent presentation of the issues concerning the publishing of scientific computations.

2. THE CONVERGENCE OF VIM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS

Consider the following system of ordinary differential equations

∗e-mail: [email protected].

(2)

(1)







x⁰₁(t) =f₁(t, x₁(t), . . . , xm(t)) x₁(t₀) =x⁰₁ ...

x⁰_m(t) =fm(t, x₁(t), . . . , xm(t)) xm(t₀) =x⁰_m wheret∈[t₀, t_f] with t₀ < t_f <∞.

We shall use the notations

x= (x₁, . . . , x_m), kxk₁=

m

X

j=1

|x_j|

x=x(t), kxk_∞= max

t kx(t)k₁. Thus, any equation of (1) may be rewritten as

x⁰_i(t) =f_i(t,x(t)), i∈ {1, . . . , m}.

The following hypothesis are introduced:

• The functionsf₁, . . . , fmare continuous and have first and second order partial derivatives inx₁, . . . , xm.

• There existsL >0 such that for anyi∈ {1, . . . , m}

|f_i(t,x)−f_i(t,y)| ≤L

m

X

j=1

|x_j−y_j|=Lkx−yk₁, ∀ x,y∈R^m.

As a consequence

∂fi(t,x)

∂xj

=|f_i_xj(t,x)| ≤L, ∀(t,x)∈[t₀, tf]×R^m, ∀i, j∈ {1, . . . , m}.

According to the VIM the sequences of approximations are (2) un+1,i(t) =un,i(t) +

Z t t0

λi(s)(u⁰_n,i(s)−fi(s,u_n(s)))ds, n∈N, i∈ {1, . . . , m}and where u_n= (u_n,₁, . . . , u_n,m).

It is supposed thatu_n,i(t₀) =x⁰_i and thatu_n,iis a continuous differentiable function for anyi∈ {1, . . . , m}.

In this case the VIM is a little trickier: the Lagrange multiplier attached to thei-th equation will act only onxi [5].

Denoting x(t) = (x₁(t), . . . , xm(t)) the solution of the initial value problem (1), ifu_n,i(t) =x_i(t) +δu_n,i(t) andu_n₊₁_,i(t) =x_i(t) +δu_n₊₁_,i(t) butu_n,j(t) =

(3)

δu_n₊₁_,i(t) =

=δun,i(t) + Z t

t0

λi(s)x⁰_i(s) +δu⁰_n,i(s)−

−f_i(s, x₁(s), . . . , xi−1(s), x_i(s) +δu_n,i(s), x_i₊₁(s), . . . , x_m(s))) ds=

=δu_n,i(t) + Z t

t0

λ_i(s)x⁰_i(s) +δu⁰_n,i(s)−f_i(s,x(s))−f_i_xi(s,x(s))δu_n,i(s)ds+

+O((δu_n,i)²) =

=δu_n,i(t) + Z t

t0

λ_i(s)δu⁰_n,i(s)−f_i_xi(s,x(s))δu_n,i(s)ds+O((δu_n,i)²).

After the integration by parts the above equality becomes δu_n₊₁_,i(t) =

(1 +λ_i(t))δu_n,i(t)− Z t

t0

λ⁰_i(s) +f_i_xi(s,x(s))λ_i(s)δu_n,i(s)ds+O((δu_n)²).

In order thatu_n₊₁_,i be a better approximation thanu_n,i,it is required thatλ_i is the solution of the following initial value problem

λ⁰(s) = −f_i_xi(s,x(s))λ(s), s∈[t₀, t], (3)

λ(t) = −1.

(4)

Because x(s) is an unknown function, the following problem is considered instead of (3)–(4)

λ⁰(s) = −f_x(s,u_n(s))λ(s), s∈[t₀, t], (5)

λ(t) = −1

(6)

with the solution denotedλ_n,i(s, t).The solution is λ_n,i(s, t) =−e^R

t

sf_ixi(τ,un(τ))dτ

and

|λ_n,i(s, t)| ≤e^L⁽^t−s⁾≤e^LT, ∀t₀ ≤s≤t≤t_f and T =t_f −t₀. The recurrence formula (2) becomes

(7) u_n₊₁_,i(t) =u_n,i(t) + Z t

t0

λ_n,i(s, t)(u⁰_n,i(s)−f_i(s,u_n(s)))ds, n∈N, for anyi∈ {1, . . . , m}.

The convergence result is:

Theorem 1. If the hypotheses stated above are valid, then the sequence (u_n)_n∈_N defined by (7) converges uniformly to x(t),the solution of the initial value problem (1).

(4)

Proof. Subtracting the equality xi(t) =xi(t) +

Z t

0 λn,i(s, t) x⁰_i(s)−fi(s,x(s))ds from (7) leads to

e_n₊₁_,i(t) =e_n,i(t) + Z t

t0

λ_n,i(s, t)e⁰_n,i(s)−(f_i(s,u_n(s))−f_i(s,x(s)))ds, whereen,i(t) =un,i(t)−x_i(t), i∈ {1, . . . , m}ande_n(t) = (en,1(t), . . . , en,m(t)) = u_n(t)−x(t), n∈N.

Again, an integration by parts gives e_n₊₁_,i(t) =

Z t t0

λ_n,i(s, t)f_i_xi(s,u_n(s))e_n,i(s)−(f_i(s,u_n(s))−f_i(s,x(s)))ds.

The hypothesis onfi implies the inequality

f_i_xi(s,u_n(s))e_n,i(s)−(f_i(s,u_n(s))−f_i(s,x(s))≤L|e_n,i(s)|+Lke_n(s)k₁ and consequently

|e_n₊₁_,i(t)| ≤Le^LT Z t

t0

(|e_n,i(s)|+ke_n(s)k₁)ds.

Summing these inequalities, fori= 1 :m,we find (8) ke_n₊₁(t)k₁ ≤(m+ 1)Le^LT

Z t t0

ke_n(s)k₁ds.

LetM = (m+ 1)Le^LT.From (8) we obtain successively:

Forn= 0 ke₁(t)k₁≤M

Z t t0

ke₀(s)k₁ds≤M(t−t₀)ke₀k_∞ ⇒ ke₁k_∞≤M Tke₀k_∞. Forn= 1

ke₂(t)k₁ ≤M Z t

t0

ke₁(s)k₁ds≤ ^M²⁽^t−t₂ ⁰⁾²ke₀k_∞ ⇒ ke₂k_∞≤ ^M²₂^T²ke₀k_∞. Inductively, it results that

ke_n(t)k₁≤M Z t

t0

ke_n−₁(s)k₁ds≤ ^Mⁿ⁽^t−t_n_! ⁰⁾ⁿke₀k_∞ ⇒ ke_nk_∞≤ ^Mⁿ_n^T_!ⁿke₀k_∞.

and hence limn→∞ke_nk∞= 0.

A numerical implementation requires the usage of a quadrature method to compute the integral in (7) and the Lagrange multipliers.

(5)

The target of the given examples is twofold: to exemplify the convergence of the VIM and to obtain some clues about the usage of numerical vs. symbolical computations of VIM.

Example 1.

x⁰(t) = 2x(t) +t x(0) = 0

The solutionx(t) = ¹₄(e²^t−2t−1) is obtained in an iteration with theMath-

ematica code provided in Appendix A.

Example 2.

x⁰(t) = 1−x²(t) x(0) = 0

The initial value problem has the solution x(t) = ^e_e^2t2t⁻+1¹. The code used previously does not give an acceptable result in a reasonable time.

Moving on numerical computation we obtain practical results. The relation (7) is transformed into

u_n₊₁_,i(t) = Z t

t0

f_i(s,u_n(s))−f_i_xi(s,u_n(s))u_n,i(s)e Rt

sf_ixi(τ,un(τ))dτ

ds+

(9)

+e Rt

t0f_ixi(τ,un(τ))dτ

x⁰.

Let (t_i)₀≤i≤I be an equidistant grid on [t₀, t_f] and denote byu_i an approximation ofx(ti).Furthermore, the recurrence relation (9) is used only to compute u_i₊₁ fromu_i,i.e. on a [t_i, t_i₊₁] interval. The integrals in (9) will be computed with the trapezoidal rule using a local equidistant grid.

The iterations are done until the distance between two consecutive approximation of u_i₊₁ is less then a given tolerance.

The final approximate solution is a first order spline function defined by the points (t_i,u_i)₀≤i≤I.

This procedure requires a single passage fromt₀ tot_f.

All our numerical results were computed using theScilabcode presented in Appendix B [11].

Fort_f = 1, I = 100 we obtained max₀_≤i≤I|u_i−x(t_i)| ≈0.4×10⁻⁶. Example 3. [5]

x˙₁ = 4x₁+ 6x₂+ 6x₃ x₁(0) = 7 x˙₂ =x₁+ 3x₂+ 2x₃ x₂(0) = 2 x˙₃ =−x₁−5x₂−2x₃ x₃(0) =−⁹₂

The solution is x₁(t) = 4e^t+ 3(1 +t)e²^t, x₂(t) = e^t+ (1 +t)e²^t, x₃(t) =

−3e^t−(³₂ + 2t)e²^t.

Fort_f = 1, I = 100 we found max₀_≤i≤Iku_i−x(t_i)k ≈0.8363×10⁻³.

(6)

Example 4.

x⁰₁(t) =x₂(t), x₁(0) = ¹₂ x⁰₂(t) = ^x_x¹⁽2^t⁾^x²²⁽^t⁾

1(t)−1 , x₂(0) =

√3 2

⇔

(x²(t)−1)x⁰⁰(t) =x(t)x⁰²(t) x(0) = ¹₂

x⁰(0) =

√3 2

The solution isx₁(t) = sin (t+^π₆), x₂(t) = cos (t+^π₆).

Fort_f = ^π₃, I = 100 the result is max₀_≤i≤Iku_i−x(t_i)k ≈0.62×10⁻⁵. Example 5. Van der Pol equation

x⁰₁(t) =x₂(t), x₁(0) = 0.5 x⁰₂(t) = (1−x²₁(t))x₂(t)−x₁(t), x₂(0) = 0 ⇔

⇔

x⁰⁰(t)−(1−x²(t))x⁰(t) +x(t) = 0 x(0) = 0.5

x⁰(0) = 0

In this case we do not have a closed form of the solution. We compare the VIM approximation with the solutionvobtained withode, aScilabnumerical integration function.

The obtained results are given in Table 1.

t_f I max₀_≤i≤Iku_i−v_ik

10 100 0.0001988

20 100 0.0033857

30 100 0.0142215

40 100 0.0384137

50 100 0.0777549

100 100 0.6410885 100 1000 0.0069729

Table 1. Numerical results for the Van der Pol equation.

4. CONCLUSIONS

Despite the convergence properties of the method the amount of the computation is greater than of the usual methods (e.g. Runge-Kutta, Adams type methods). Even so the numerical solution can be taken into consideration.

The numerical implementation can be improved by an adaptive approach and using some parallel techniques (e.g. OpenCL / CUDA) in an appropriate environment.

Although the VIM may be implemented for symbolic computation our ex- periments show disappointing results.

The VIM offers a way to obtain a symbolic approximation of the solution of the initial value problem. But such an approximation may also be obtained from a numerical solution with theEureqasoftware [10], [7]. A better symbolic implementation would be useful.

(7)

suggestions for the improvement of this paper.

REFERENCES

[1] T. Daly,Publishing Computational Mathematics, Notices of the AMS,59(2012) no.

2, pp. 320–321.

[2] M. Inokuti, H. Sekine, T. Mura,General use of the Lagrange multiplier in Nonlinear Mathematical Physics, In Variational Methods in Mechanics and Solids, ed. Nemat- Nasser S., Pergamon Press, pp. 156–162, 1980.

[3] J.H. He,Variational iteration method - Some recent results and new interpretations, J. Comput. Appl. Math., 207(2007), pp. 3–17.

[4] Z.M. Odibat,A study on the convergence of variational iteration method, Math. Com- puter Modelling,51(2010), pp. 1181–1192.

[5] D. K. Salkuyeh,Convergence of the variational iteration method for solving linear systems of ODE with constant coefficients, Comp. Math. Appl., 56(2008), pp. 2027–

2033.

[6] D. K. Salkuyeh, A. Tavakoli, Interpolated variational iteration method for initial value problems, arXiv:1507.01306v1, 2015.

[7] E. Scheiber, From the numerical solution to the symbolic form, Bull. Transilvania University of Bra¸sov, Series III, Mathematics, Informatics, Physics,8(57) (2015) no. 1, pp. 129–137.

[8] M. Tatari, M. Dehghan,On the convergence of He’s variational iteration method.J.

Comput. Appl. Math.,207(2007), pp. 121–128.

[9] M. Torvattanabun, S. Koonprasert, Convergence of the variational iteration method for solving a first-order linear system of PDEs with constant coefficients, Thai J. of Mathematics, Special Issue, pp. 1–13, 2009.

[10] * * *,www.nutonian.com [11] * * *,www.scilab.org

Appendix ^A. MATHEMATICACODE

The Mathematicaprocedure used to solve an initial value problem.

1 In[ 1 ] : =

2 VIM[ f , U0 , m ] := Module[{V, U = U0 , df , Lambda}, 3 df [ t , x ] := D[ f [ t , x ] , x ] ;

4 Lambda [ U ] := −Exp[

5 Integrate[ df [w, x ] / . {x −> U, t −> w}, {w, s , t}] ] ; 6 For[ i = 0 , i < m, i ++,

7 V = U +

8 Integrate[ Lambda [U] ( (D[U, t ] − f [ t , U] ) / . t −> s ) , {s , 0 , t}] ; 9 U = V; Clear[V ] ] ; U]

For Example 1 the calling sequence is

1 In[ 2 ] : = f [ t , x ] := {2 x [ [ 1 ] ] + t} 2 In[ 3 ] : = U0 = {0}

3 In[ 4 ] : = VIM[ f , U0 , 1 ]

4 Out[4]={(1/4)∗(−1 + Eˆ(2∗ t ) − 2∗ t )}

(8)

Appendix ^B. SCILABCODE

The code used to obtain numerical results

1 function [ t , u ,i n f o]=vim ( f , df , x0 , t0 , t f , n , nmi , t o l ) 2 // Computes the s o l u t i o n o f an i n i t i a l value problem 3 // o f a system o f or din ary d i f f e r e n t i a l e q u a t i o n s 4 // using the V a r i a t i o n a l I t e r a t i o n Method − 5 // Dispatcher function.

6 //

7 // C a l l i n g Sequence

8 // [ t , u ,i n f o]=vim ( f , df , x0 , t0 , t f , n , nmi , t o l ) 9 //

10 // Output arguments

11 // t : a r e a l vector , the times at which the s o l u t i o n i s computed . 12 // u : a r e a l v e c t o r or matrix , the numerical s o l u t i o n .

13 // i n f o: an i n t e g e r , error code (0 −OK) . 14 // Input arguments

15 // f : a S c i l a b function, the r i g h t hand s i z e o f the 16 // d i f f e r e n t i a l system .

17 // df : a S c i l a b function, df =( d f 1 / dx 1 , d f 2 / dx 2 , . . . ) . 18 // x0 : a r e a l vector , the i n i t i a l c o n d i t i o n s .

19 // t0 : a r e a l s c a l a r , the i n i t i a l time . 20 // t f : a r e a l s c a l a r , the f i n a l time .

21 // n : a p o s i t i v e i n t e g e r , the g l o b a l d i s c r e t i z a t i o n parameter . 22 // nmi : maximum allowed number o f l o c a l i t e r a t i o n s i t e r a t i o n . 23 // t o l : a p o s i t i v e r e a l number , a t o l e r a n c e .

24 //

25 t=l i n s p a c e( t0 , t f , n ) 26 d=length( x0 )

27 u=zeros(d , n ) 28 u ( : , 1 ) = x0 ’

29 m=10 // the d i s c r e t i z a t i o n parameter on [ t i , t {i +1}]

30 u0=x0

31 errorMarker=%t

32 f o r i =1:n−1 do // the passage from t 0 to t f 33 // Computes u0:= u {i +1} from u i

34 [ u0 ,eM, i t e r ]= vimstep ( f , df , u0 , t ( i ) , t ( i +1) ,m, nmi , t o l ) 35 u ( : , i +1)=u0 ’

36 errorMarker=errorMarker & eM

37 end

38 i f errorMarker then 39 i n f o=0

40 e l s e 41 i n f o=1

42 end

43 en dfu nc tio n

(9)

1 function [ u ,eM, i t e r ]= vimstep ( f , df , x0 , t0 , t f , n , nmi , t o l ) 2 // V a r i a t i o n a l I t e r a t i o n Method .

3 //

4 // Output arguments

5 // u : a r e a l v e c t o r or matrix , the numerical s o l u t i o n . 6 // eM : a boolean value , error marker , (%t − OK) .

7 // i t e r : an i n t e g e r , the number o f l o c a l performed i t e r a t i o n s . 8 // Input arguments

9 // f : a S c i l a b function, the r i g h t hand s i z e o f the 10 d i f f e r e n t i a l system .

11 // df : a S c i l a b function, df =( d f 1 / dx 1 , d f 2 / dx 2 , . . . ) . 12 // x0 : a r e a l vector , the i n i t i a l c o n d i t i o n s .

13 // t0 : a r e a l s c a l a r , the i n i t i a l time . 14 // t f : a r e a l s c a l a r , the f i n a l time .

15 // n : a p o s i t i v e i n t e g e r , the l o c a l d i s c r e t i z a t i o n parameter . 16 // nmi : maximum allowed number o f l o c a l i t e r a t i o n s i t e r a t i o n . 17 // t o l : a p o s i t i v e r e a l number , a t o l e r a n c e .

18 //

19 t=l i n s p a c e( t0 , t f , n ) 20 h=( t f−t0 ) / ( n−1) 21 d=length( x0 ) 22 u o l d=zeros(d , n ) 23 u o l d ( : , 1 ) = x0 ’

24 sw=%t

25 i t e r =0 26 while sw do 27 i t e r=i t e r +1

28 u new=zeros(d , n )

29 u new ( : , 1 ) = x0 ’ 30 f 0=zeros(d , n ) 31 df0=zeros(d , n ) 32 f o r j =1:n do

33 p=u o l d ( : , j )

34 q=f ( t ( j ) , p )

35 dq=df ( t ( j ) , p )

36 f 0 ( : , j )=q

37 df0 ( : , j )=dq

38 end

39 f o r i =2:n do

40 z=zeros(d , n )

41 f o r j=i−1:−1:1 do

42 z ( : , j )=0.5∗h ∗( df0 ( : , j )+df0 ( : , j +1))+z ( : , j +1)

43 end

44 w=(f0−df0 . ∗ u o l d ) . ∗exp( z )

45 s=zeros(d , 1 )

46 s=w( : , 1 ) +w( : , i )

47 i f i>2 then

48 f o r j =2: i−1 do

49 s=s+2∗w( : , j )

50 end

51 end

52 u new ( : , i )=0.5∗h∗ s+exp( z ( : , 1 ) ) . ∗ x0 ’

53 end

54 nrm=max(abs( u new−u o l d ) )

55 u o l d=u new

(10)

56 i f nrm<t o l | i t e r>=nmi then

57 sw=%f

58 end

59 end

60 u=u o l d ( : , n ) ’ 61 i f nrm<t o l then

62 eM=%t

63 e l s e

64 eM=%f

65 end

66 en dfu nction

The calling sequence for Example 4 is

1 d e f f ( ’ q=f ( t , p ) ’ ,

2 [ ’ x1=p ( 1 ) ’ , ’ x2=p ( 2 ) ’ , ’ q(1)= x2 ’ , ’ q(2)= x1 . ∗ x2 . ˆ 2 . 0 . / ( x1 .ˆ2−1) ’ ] ) 3 d e f f ( ’ q=df ( t , p ) ’ ,

4 [ ’ x1=p ( 1 ) ’ , ’ x2=p ( 2 ) ’ , ’ q(1)=0 ’ , ’ q(2)=2∗ x1 . ∗ x2 . / ( x1 .ˆ2−1) ’ ] ) 5 x0 = [ 0 . 5 , 0 . 5 ∗sqrt( 3 ) ]

6 t0=0 7 t f=%p i /3 8 n=100 9 nmi=50 10 t o l =1e−5

11 [ t , u ,i n f o]=vim ( f , df , x0 , t0 , t f , n , nmi , t o l )

Received by the editors: October 5, 2015.