\(^\ast \)e-mail: scheiber@unitbv.ro.

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

where \(t\in [t_0,t_f]\) with \(t_0{\lt}t_f{\lt}\infty .\)

We shall use the notations

Thus, any equation of (1) may be rewritten as

The following hypothesis are introduced:

• The functions \(f_1,\ldots ,f_m\) are continuous and have first and second order partial derivatives in \(x_1,\ldots ,x_m.\)

• There exists \(L{\gt}0\) such that for any \(i\in \{ 1,\ldots ,m\} \)

  \[ |f_i(t,\mathbf{x})-f_i(t,\mathbf{y})|\le L\sum _{j=1}^m|x_j-y_j|=L\| \mathbf{x}-\mathbf{y}\| _1, \quad \forall \ \mathbf{x},\mathbf{y}\in \mathbb {R}^m. \]

  As a consequence

  \[ \Big\vert \frac{\partial f_i(t,\mathbf{x})}{\partial x_j}\Big\vert =|f_{i_{x_j}}(t,\mathbf{x})|\le L, \quad \forall (t,\mathbf{x})\in [t_0,t_f]\times \mathbb {R}^m,\ \forall i,j\in \{ 1,\ldots ,m\} . \]

According to the VIM the sequences of approximations are

\(i\in \{ 1,\ldots ,m\} \) and where \(\mathbf{u}_n=(u_{n,1},\ldots ,u_{n,m}).\)

It is supposed that \(u_{n,i}(t_0)=x^0_i\) and that \(u_{n,i}\) is a continuous differentiable function for any \(i\in \{ 1,\ldots ,m\} .\)

In this case the VIM is a little trickier: the Lagrange multiplier attached to the \(i\)-th equation will act only on \(x_i\) [ 5 ] .

Denoting \(\mathbf{x}(t)=(x_1(t),\ldots ,x_m(t))\) the solution of the initial value problem (1), if \(u_{n,i}(t)=x_i(t)+\delta u_{n,i}(t)\) and \(u_{n+1,i}(t)=x_i(t)+\delta u_{n+1,i}(t)\) but \(u_{n,j}(t)=x_j(t),\) for \(j\not=i,\) then (2) implies

After the integration by parts the above equality becomes

In order that \(u_{n+1,i}\) be a better approximation than \(u_{n,i},\) it is required that \(\lambda _i\) is the solution of the following initial value problem

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

with the solution denoted \(\lambda _{n,i}(s,t).\) The solution is

and

The recurrence formula (2) becomes

for any \(i\in \{ 1,\ldots ,m\} .\)

The convergence result is:

from (7) leads to

where \(e_{n,i}(t)=u_{n,i}(t)-x_i(t),\ i\in \{ 1,\ldots ,m\} \) and \(\mathbf{e}_n(t)= (e_{n,1}(t),\ldots ,e_{n,m}(t))=\mathbf{u}_n(t)-\mathbf{x}(t),\ n\in \mathbb {N}.\)

Again, an integration by parts gives

The hypothesis on \(f_i\) implies the inequality

and consequently

Summing these inequalities, for \(i=1:m,\) we find

Let \(M=(m+1)Le^{LT}.\) From (8) we obtain successively:

For \(n=0\)

For \(n=1\)

Inductively, it results that

and hence \(\lim _{n\rightarrow \infty }\| \mathbf{e}_n\| _{\infty }=0.\)

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

3 Computational results

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.

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 experiments 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 the Eureqa software [ 10 ] , [ 7 ] . A better symbolic implementation would be useful.

A Mathematica code

The Mathematica procedure used to solve an initial value problem.

In[1]:=
 VIM[f_, U0_, m_] := Module[{V, U = U0, df, Lambda},
  df[t_, x_] := D[f[t, x], x];
  Lambda[U_] := -Exp[
     Integrate[df[w, x] /. {x -> U, t -> w}, {w, s, t}]];
  For[i = 0, i < m, i++,
   V = U +
     Integrate[Lambda[U] ((D[U, t] - f[t, U]) /. t -> s), {s, 0, t}];
   U = V; Clear[V]]; U]

In[2]:= f[t_, x_] := {2 x[[1]] + t}
In[3]:= U0 = {0}
In[4]:= VIM[f,U0,1]
Out[4]={(1/4)*(-1 + E^(2*t) - 2*t)}

B Scilab code

The code used to obtain numerical results

function [t,u,info]=vim(f,df,x0,t0,tf,n,nmi,tol)
// Computes the solution of an initial value problem
// of a system of ordinary differential equations
// using the Variational Iteration Method -
// Dispatcher function.
//
// Calling Sequence
//    [t,u,info]=vim(f,df,x0,t0,tf,n,nmi,tol)
//
// Output arguments
//    t   : a real vector, the times at which the solution is computed.
//    u   : a real vector or matrix, the numerical solution.
//    info: an integer, error code (0 - OK).
// Input arguments
//    f   : a Scilab function, the right hand size of the
//          differential system.
//    df  : a Scilab function, df=(df_1/dx_1,df_2/dx_2,...).
//    x0  : a real vector, the initial conditions.
//    t0  : a real scalar, the initial time.
//    tf  : a real scalar, the final time.
//    n   : a positive integer, the global discretization parameter.
//    nmi : maximum allowed number of local iterations iteration.
//    tol : a positive real number, a tolerance.
//
    t=linspace(t0,tf,n)
    d=length(x0)
    u=zeros(d,n)
    u(:,1)=x0'
    m=10   // the discretization parameter on [t_i,t_{i+1}]
    u0=x0
    errorMarker=%t
    for i=1:n-1 do  // the passage from t_0 to t_f
        // Computes u0:=u_{i+1} from u_i
        [u0,eM,iter]=vimstep(f,df,u0,t(i),t(i+1),m,nmi,tol)
        u(:,i+1)=u0'
        errorMarker=errorMarker & eM
    end
    if errorMarker then
        info=0
    else
        info=1
    end
endfunction

function [u,eM,iter]=vimstep(f,df,x0,t0,tf,n,nmi,tol)
// Variational Iteration Method.
//
// Output arguments
//    u    : a real vector or matrix, the numerical solution.
//    eM   : a boolean value, error marker, (%t - OK).
//    iter : an integer, the number of local performed iterations.
// Input arguments
//    f   : a Scilab function, the right hand size of the
            differential system.
//    df  : a Scilab function, df=(df_1/dx_1,df_2/dx_2,...).
//    x0  : a real vector, the initial conditions.
//    t0  : a real scalar, the initial time.
//    tf  : a real scalar, the final time.
//    n   : a positive integer, the local discretization parameter.
//    nmi : maximum allowed number of local iterations iteration.
//    tol : a positive real number, a tolerance.
//
    t=linspace(t0,tf,n)
    h=(tf-t0)/(n-1)
    d=length(x0)
    u_old=zeros(d,n)
    u_old(:,1)=x0'
    sw=%t
    iter=0
    while sw do
        iter=iter+1
        u_new=zeros(d,n)
        u_new(:,1)=x0'
        f0=zeros(d,n)
        df0=zeros(d,n)
        for j=1:n do
            p=u_old(:,j)
            q=f(t(j),p)
            dq=df(t(j),p)
            f0(:,j)=q
            df0(:,j)=dq
        end
        for i=2:n do
            z=zeros(d,n)
            for j=i-1:-1:1 do
                z(:,j)=0.5*h*(df0(:,j)+df0(:,j+1))+z(:,j+1)
            end
            w=(f0-df0.*u_old).*exp(z)
            s=zeros(d,1)
            s=w(:,1)+w(:,i)
            if i>2 then
                for j=2:i-1 do
                    s=s+2*w(:,j)
                end
            end
            u_new(:,i)=0.5*h*s+exp(z(:,1)).*x0'
        end
        nrm=max(abs(u_new-u_old))
        u_old=u_new
        if nrm<tol | iter>=nmi then
            sw=%f
        end
    end
    u=u_old(:,n)'
    if nrm<tol then
        eM=%t
    else
        eM=%f
    end
endfunction

deff('q=f(t,p)',
    ['x1=p(1)','x2=p(2)','q(1)=x2','q(2)=x1.*x2.^2.0./(x1.^2-1)'])
deff('q=df(t,p)',
    ['x1=p(1)','x2=p(2)','q(1)=0','q(2)=2*x1.*x2./(x1.^2-1)'])
x0=[0.5,0.5*sqrt(3)]
t0=0
tf=%pi/3
n=100
nmi=50
tol=1e-5
[t,u,info]=vim(f,df,x0,t0,tf,n,nmi,tol)