Arbitrary-order functionally fitted energy-diminishing methods for gradient systems
Bin Wang Ting Li Yajun Wu
math.NA 65L05








([ )
⋰⋱  


𝔹
ℂ
𝔻
𝔼
𝔽
𝔾
ℍ
𝕀
𝕁
𝕂
𝕃
𝕄
ℕ
𝕆
ℙ
ℚ
ℝ
𝕊
𝕋
𝕌
𝕍
𝕎
𝕏
𝕐
ℤ

A
B
C
D
E
F
G
H̋
I
J
K
ŁL
M
N
ØO
P
Q
R
S
T
U
V
W
X
Y
Z

d
e
diag
null
rank
rank
span
i
Re
Im

theoTheorem
lemLemma
corCorollary
remRemark
defiDefinition
criteCriterion
Proof  


a
b
c
d
e
f
g
h
i
j
k
l
m
n
o
p
q
r
s
t
u
v
w
x
y
z

0
1


α
β
γ
δ
η
φ
ϕ
λ
ψ
τ
ξ
ζ
ω



D̃
sinc







Arbitrary-order  functionally fitted energy-diminishing methods for gradient systems
    Bin Wang 
[School of Mathematical Sciences, Qufu Normal University,
Qufu 273165, P.R. China; Mathematisches Institut, University of
Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany.
The research is supported in part by the Alexander von Humboldt
Foundation and by the Natural Science Foundation of Shandong
Province (Outstanding Youth Foundation) under Grant ZR2017JL003.
E-mail: wang@na.uni-tuebingen.de ] Ting LiSchool
of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R.
China. E-mail: 1009587520@qq.com Yajun WuSchool
of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R.
China. E-mail: 1921170786@qq.com 
    December 30, 2023
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



 It is well known that for gradient systems in Euclidean space or on
a Riemannian manifold, the energy decreases monotonically along
solutions. In this letter we derive and analyse functionally fitted
energy-diminishing methods to preserve this key property of gradient
systems. It is proved that the novel methods are unconditionally
energy-diminishing and can achieve damping for very stiff gradient
systems. We also show that the methods can be of arbitrarily high
order and discuss their implementations. A numerical test is
reported to illustrate the efficiency of the  new methods in
comparison with three existing numerical methods in the literature.

Keywords:  gradient systems, energy-diminishing methods,
functionally fitted methods, arbitrary-order methods

MSC:65L05










§ INTRODUCTION

In this letter, we investigate   the  following  gradient systems in
coordinates:

    G(y(t))ẏ(t)=- ∇ U(y(t)),  
    y(0)=y_0∈ℝ^d,   t∈[0,T],

where U(y): ℝ^d→ℝ is a potential
function  and the symmetric matrix G(y) is assumed to satisfy

    v^⊺ G(y) v ≥  v^⊺Γ v>0

for all vectors v≠ 0. Here Γ is a fixed positive
definite matrix.

Gradient systems frequently  arise in a wide variety of applications
both  in a finite-dimensional and infinite-dimensional setting.
There are many examples of this system (see, e.g.
<cit.>)
such as models in quantum systems,  in differential geometry,  in
image processing, and in material science. A fundamental and key
property of gradient systems is that along every exact solution of
(<ref>),  one has

    d/dtU(y(t))= ∇
     U(y(t))^⊺ẏ(t)=-ẏ(t)^⊺G(y(t))ẏ(t)≤0,

which shows that U(y(t)) is monotonically decreasing  with t.
The monotonicity is true with strict inequality except at stationary
points of U. The aim of this letter is to formulate and analyse a
novel kind of methods preserving this monotonicity in the numerical
treatment, i.e., after one step of the method starting from y_0
with a time step h one would have 
    U(y_1)≤ U(y_0).



In order to get methods with this property,  Hairer and Lubich
 analysed various energy-diminishing  methods in <cit.>.
 They showed that implicit Euler method has this property but it  is only of order
 one.  Algebraically stable Runge-Kutta methods were proved to reduce the energy in each
step under a mild step-size restriction, which means that
Runge-Kutta methods are not unconditionally energy-diminishing.
 They also showed that discrete-gradient methods,  averaged vector field (AVF) methods and AVF collocation
methods are unconditionally energy-diminishing, but cannot achieve
damping for very stiff gradient systems. In this letter, we will
derive a novel kind of methods which  can  be
 of arbitrarily high order.
 Moreover, the methods will be shown that they
 are unconditionally energy-diminishing  and are strongly damped even
for very stiff gradient systems.

The rest of this letter is organised as follows. In Section
<ref>, we derive the novel methods and prove that they are
unconditionally
 energy-diminishing for gradient systems. The unconditionally damping
property is analysed in  Section <ref>.    We study the
order of the methods in Section <ref>. Section
<ref> is devoted to the implementation issue. In
Section <ref>, a numerical test is carried out to
 demonstrate the   excellent qualitative behavior.
Section <ref> focuses on the concluding remarks.







§ FUNCTIONALLY FITTED ENERGY-DIMINISHING METHODS

In order to formulate the novel methods, we will use the
functionally fitted technology, which  is a popular approach to
constructing  efficient and effective methods in scientific
computing (see, e.g. <cit.>). To this end,
  define a function space
Y=span{φ_0(t),…,φ_r-1(t)} on
[0,T]
  by (see <cit.>)

    Y={w : w(t)=∑_i=0^r-1φ_i(t)W_i, t∈ I,W_i∈ℝ^d},

where {φ_i(t)}_i=0^r-1 are sufficiently smooth and
linearly independent on [0,T].  In this letter, we consider   two
finite-dimensional function spaces  Y and X as follows

    Y=span{φ_0(t),…,φ_r-1(t)},  
    X=span{1,∫_0^tφ_0(s)ds,…,∫_0^tφ_r-1(s)ds}.

Choose a stepsize h and define the function spaces Y_h and
X_h on [0,1] by

    Y_h=span{φ̃_0(τ),…,φ̃_r-1(τ)},  
    X_h=span{1,∫_0^τφ̃_0(s)ds,…,∫_0^τφ̃_r-1(s)ds},

 where φ̃_i(τ)=φ_i(τ h), τ∈[0,1]
for i=0,1,…,r-1. We remark that  for all the functions
throughout this letter,  the notation f̃ (τ) is referred
to f(τ h).

We will use a projection 𝒫_h  in the formulation of
the new methods.  It was defined in <cit.> and we
summarise it here.
 (See <cit.>)
Let 𝒫_hw̃ be a projection of w̃ onto
Y_h, where w̃(τ) is a continuous
ℝ^d-valued function on [0,1].
  The definition of
𝒫_hw̃ is given by

    ⟨ṽ(τ),𝒫_hw̃(τ)⟩=⟨ṽ(τ),w̃(τ)⟩  for any ṽ(τ)∈ Y_h,

where   the inner product ⟨·,·⟩ is defined by

    ⟨w̃_1,w̃_2⟩=⟨w̃_1(τ),w̃_2(τ)⟩_τ=∫_0^1w̃_1(τ)·w̃_2(τ)dτ.

Here w̃_1 and w̃_2 are two integrable
functions (scalar-valued or vector-valued)
 on [0,1], and  `·' denotes the entrywise multiplication
operation if they are  both vector-valued functions.


The  following property of  𝒫_h will also be useful in
this letter, which has been proved in <cit.>.

(See <cit.>) The projection 𝒫_hw̃ can be explicitly expressed as

    𝒫_hw̃(τ)=⟨
    P_τ,σ,w̃(σ)⟩_σ,

where

    P_τ,σ=∑_i=0^r-1ψ̃_i(τ)ψ̃_i(σ),

and {ψ̃_0,…,ψ̃_r-1} is a
standard orthonormal basis of Y_h under the inner product
⟨·,·⟩.


Based on these preliminaries, we are in a position to present the
scheme of functionally fitted energy-diminishing methods.




Choose a stepsize h  and consider a function ũ∈ X_h
with ũ(0)=y_0, satisfying

    G(ũ(τ))ũ^'(τ) = -𝒫_h(
    ∇ U(ũ(τ))).

The numerical approximation  after one step is  defined by y_1=ũ(1). We call this method     as functionally-fitted
energy-diminishing  method and denote it by FFED.



The  FFED method (<ref>) reduces the energy of gradient systems
(<ref>), i.e., U(y_1)≤ U(y_0). This means that our
FFED method (<ref>) is unconditionally energy-diminishing.

According to the definitions of X_h and Y_h, we know that if
ũ(τ)∈ X_h, then one has ũ'(τ)∈ Y_h.
From the definition of 𝒫_h, it follows that

    ∫_0^1ũ'(τ)_i𝒫_h( ∇
    U(ũ(τ)))_idτ=∫_0^1ũ'(τ)_i( ∇ U(ũ(τ)))_idτ,   i=1,2,…,d,

where   (·)_i denotes the ith entry of a vector. Then, we
arrive at

    ∫_0^1ũ'(τ)^⊺𝒫_h( ∇
    U(ũ(τ)))
    dτ=∫_0^1ũ'(τ)^⊺∇
    U(ũ(τ)) dτ.

Therefore, it is obtained that

    U(y_1)-U(y_0)=∫_0^1d/dτU(ũ(τ))dτ
     =h∫_0^1ũ'(τ)^⊺∇
    U(ũ(τ))dτ=h∫_0^1ũ'(τ)^⊺𝒫_h(
    ∇ U(ũ(τ)))dτ.

Inserting the scheme (<ref>) into this formula yields

    U(y_1)-U(y_0)   =-h∫_0^1ũ'(τ)^⊺G(ũ(τ))ũ^'(τ)dτ≤0.

 



§ UNCONDITIONALLY   DAMPING
PROPERTY
 In this section, we consider the
 potential U of  the form

    U(y)=12y^⊺Ay+V(y),

  where the function V  is  twice continuously differentiable  and the matrix A  is   symmetric positive semi-definite
and of arbitrarily large norm. In this case, (<ref>) is a
stiff gradient system.   This kind of stiff systems arise from the
spatial discretization of Cahn–Hilliard and Allen–Cahn partial
differential equations (see, e.g. <cit.>). Many
effective methods have been derived for this stiff gradient system
with a constant matrix G and we refer to
<cit.>
for example.

 The FFED method (<ref>) for solving this stiff gradient system is defined  as
 follows.

We consider a function ũ∈ X_h with
ũ(0)=y_0, satisfying

    G(ũ(τ))
    ũ^'(τ)+Aũ(τ)=-𝒫_h(
    ∇ V(ũ(τ))),     y_1=ũ(1).

This method is denoted  by EFFED.




The EFFED  method (<ref>) reduces the energy of the stiff
gradient system (<ref>) with (<ref>), i.e.,
U(y_1)≤ U(y_0).


This proof is similar to that of Theorem <ref>. For any
w̃ (τ)∈ Y_h, it is easy to prove that

    ∫_0^1w̃(τ)^⊺G(ũ(τ))ũ^'(τ)dτ=-∫_0^1w̃(τ)^⊺Aũ(τ)dτ
    -∫_0^1w̃(τ)^⊺𝒫_h( ∇
    V(ũ(τ))) dτ
       =-∫_0^1w̃(τ)^⊺Aũ(τ)dτ
    -∫_0^1w̃(τ)^⊺∇
    V(ũ(τ)) dτ= -∫_0^1w̃(τ)^⊺∇ U(ũ(τ)) dτ.

Considering ũ^'(τ)∈ Y_h and letting
w̃(τ)= ũ^'(τ) in (<ref>), we
obtain

    ∫_0^1ũ^'(τ)^⊺G(ũ(τ))ũ^'(τ)dτ=- ∫_0^1ũ^'(τ)^⊺∇ U(ũ(τ)) dτ.

Thus, one arrives

    U(y_1)-U(y_0)   =∫_0^1d/dτU(ũ(τ))dτ
    =h∫_0^1ũ'(τ)^⊺∇
    U(ũ(τ))dτ
    =-h∫_0^1ũ^'(τ)^⊺G(ũ(τ))ũ^'(τ)dτ≤0.

 


For G(y)=I and a quadratic potential (i.e., V(y)=0 in
(<ref>)),  our   EFFED method (<ref>) becomes 
    ũ^'(τ)+Aũ(τ)=0,
 which leads to 
    y_1=ũ(1)=e^ -h  Ay_0.
 This scheme has been presented
and researched in <cit.>. Rewrite this scheme as y_1=R(-hA)y_0 with  the stability function R(-hA)=e^ -h  A.
It is noted that the damping property
 |R(∞)| <1 plays an important role in   the  properties of
Runge-Kutta methods when solving  semilinear parabolic equations,
which has been researched in <cit.>. The role of
the condition  |R(∞)| =0 has been well understood in Chapter
VI of <cit.>. It has been shown in <cit.>
that all the unconditionally energy-diminishing methods including
discrete-gradient methods, AVF methods and AVF collocation methods
show no damping for very stiff gradient systems since they do not
satisfy damping property unconditionally.
  However,  it is noted  that  for our EFFED method (<ref>), one  has 
    |R(∞)| =|e^
    -∞|=0.
 This means that our  methods have unconditionally
damping property, which  is a significant feature especially  for
very stiff gradient systems.





§ ALGEBRAIC ORDER

In this section, we analyse the algebraic order of the new methods.


The FFED method (<ref>) and EFFED method (<ref>) are
both of order 2r, which means that

    ũ(1)-y(t_0+h)=𝒪(h^2r+1).


 Denote by y(·,t̃,
ỹ) the solution of y'=-G^-1(y) ∇ U(y) satisfying
the initial condition y(t̃,t̃, ỹ)=ỹ
for any given t̃∈[0,h]. Let Φ(s,t̃,
ỹ)=∂ y(s,t̃, ỹ)/∂ỹ. Recalling the elementary theory of ordinary
differential equations, the following standard result is obtained

    ∂ y(s,t̃, ỹ)/∂t̃=-Φ(s,t̃, ỹ)(-G^-1(ỹ)
    ∇ U(ỹ)).

On the base of this result, we get

    ũ(1)-y(t_0+  h)=y(t_0+  h,t_0+  h,ũ(1))-y(t_0+  h,t_0,y_0)
       =∫_0^1d/dαy(t_0+  h,t_0+α h,ũ(α))dα
       =∫_0^1(h∂ y/∂t̃(t_0+
    h,t_0+α h,ũ(α))+
    ∂ y/∂ỹ(t_0+  h,t_0+α h,ũ(α))hũ^'(α))dα
       =∫_0^1(-h∂ y/∂ỹ(t_0+
    h,t_0+α
    h,ũ(α))G^-1(ũ(α))(Aũ(α)- ∇ U(ũ(α)))
         + ∂ y/∂ỹ(t_0+ h,t_0+α
    h,ũ(α))G^-1(ũ(α))(hA
    ũ(α)-h⟨
    P_1,σ, ∇ U(ũ(α))⟩_α))dα
       =-h∫_0^1Φ^1(α)G^-1(ũ(α))(g(ũ(α))-𝒫_h(g∘ũ)(α))dα,

where Φ^1(α)=∂ y/∂ỹ(t_0+
h,t_0+α h,ũ(α)).
 We denote the matrix Φ^1(α)G^-1(ũ(α))  by K(α) and
 partition it as
K(α)=(K_1(α),…,K_d(α))^⊺. It
follows from Lemma 3.4 of  <cit.> that

    K_i(α)=𝒫_hK_i(α)+𝒪(h^r),  
    i=1,2,…,d,

and

    g(ũ(α))-𝒫_h(g∘ũ)(α)=𝒪(h^r).

 On the other hand, in light of the definition of 𝒫_h, we
 obtain

    ∫_0^1(𝒫_hK_i(α))^⊺g(ũ(α))dα
    =∫_0^1(𝒫_hK_i(α))^⊺𝒫_h(g∘ũ)(α)dα,   i=1,2,…,d.

 Therefore, one has

    ũ(1)-y(t_0+h)
    = -h∫_0^1(([ (𝒫_hK_1(α))^⊺;             ⋮; (𝒫_hK_d(α))^⊺ ])+𝒪(h^r))(g(ũ(α))-𝒫_h(g∘ũ)(α))dα
    
    =   -h∫_0^1([ (𝒫_hK_1(α))^⊺(g(ũ(α))-𝒫_h(g∘ũ)(α));                                  ⋮; (𝒫_hK_d(α))^⊺(g(ũ(α))-𝒫_h(g∘ũ)(α)) ])dα
    -h∫_0^1𝒪(h^r)×𝒪(h^r)dα
    
    =   0+𝒪(h^2r+1)=𝒪(h^2r+1).

 

It follows from this theorem that our new methods  can be of
arbitrarily high order   only by   choosing a large integer r,
which is  very convenient  and simple.





§ IMPLEMENTATIONS

This section considers the implementation issue  of the new methods.


 §.§ For the case that G(y) is a constant matrix

We first  consider the case that G(y) is a constant matrix M.
Under this condition, the FFED method (<ref>) becomes

    ũ^'(τ)=-M^-1∫_0^1 P_τ,σ∇
    U(ũ(σ))dσ,

which can be solved by the variation-of-constants formula  as
follows

    ũ(τ)
    = y_0-τ h  M^-1∫_0^1  ∫_0^1 P_ξτ,σ∇
    U(ũ(σ))dσ dξ
       = y_0-τ h M^-1∑_i=0^r-1(
     ∫_0^1  ψ̃_i(ξτ)dξ) ∫_0^1 ψ̃_i(σ)∇
    U(ũ(σ)) dσ.


 Following
<cit.>, we introduce the generalized Lagrange
interpolation functions l_i(τ)∈ X_h with respect to
(r+1) distinct points {c_i}_i=1^r+1⊆[0,1]:

    (l_1(τ),…,l_r+1(τ))=(Φ_1(τ),Φ_2(τ),…,Φ_r+1(τ))Λ^-1,

where {Φ_i(t)}_i=1^r+1 is a basis of X,
Φ_i(τ)=Φ_i(τ h) and

    Λ=([     Φ_1(c_1)     Φ_2(c_1)            …   Φ_r+1(c_1);     Φ_1(c_2)     Φ_2(c_2)            …   Φ_r+1(c_2);            ⋮            ⋮                         ⋮;   Φ_1(c_r+1)   Φ_2(c_r+1)            … Φ_r+1(c_r+1);  ]).

   Then it follows from
<cit.> that {l_i(τ)}_i=1^r+1 is a basis
of X_h
 satisfying l_i(c_j)=δ_ij. Since ũ(τ)∈
X_h, ũ(τ) can be expressed by the basis of X_h
as

    ũ(τ)=∑_i=1^r+1ũ(c_i)l_i(τ).

Choosing 0= c_1<c_2<⋯<c_r+1= 1 and denoting
y_σ=ũ(σ), we present the following practical
FFED  methods.


The practical  FFED  method  is defined  as follows:

    {   y_c_j= y_0-c_jhM^-1∑_i=0^r-1(
     ∫_0^1  ψ̃_i(ξ c_j)dξ) ∫_0^1ψ̃_i(σ)∇
    U(∑_m=1^r+1y_c_ml_m(σ)) dσ,   j=2,…,r,
       y_1= y_0- h M^-1∑_i=0^r-1(
     ∫_0^1  ψ̃_i(ξ  )dξ) ∫_0^1ψ̃_i(σ)∇
    U(∑_m=1^r+1y_c_ml_m(σ)) dσ.
    .




In a similar way, for the stiff system with the special potential
U of the form (<ref>), we have the following practical EFFED
method for (<ref>).

The practical  EFFED  method  is given  by

    {   y_c_j=e^-c_j h M^-1Ay_0-c_jhM^-1∑_i=0^r-1∫_0^1(
     ∫_0^1 e^-(1-ξ)c_j h M^-1 Aψ̃_i(ξ c_j)dξ) ψ̃_i(σ)∇
    V(∑_m=1^r+1y_c_ml_m(σ)) dσ,
                                                                      j=2,…,r,
       y_1=e^ -h M^-1Ay_0- hM^-1∑_i=0^r-1∫_0^1(
     ∫_0^1 e^-(1-ξ)  h M^-1 Aψ̃_i(ξ  )dξ) ψ̃_i(σ)∇
    V(∑_m=1^r+1y_c_ml_m(σ)) dσ.
    .



It is noted that when M=I, this method has been proposed and
analysed  in <cit.>.




 §.§ For the general case that  G(y) depends on y

Now we pay attention to the general case that G(y) is a  matrix
depending on y. The FFED method (<ref>) becomes

    ũ^'(τ)=-G^-1(ũ(τ))∫_0^1
    P_τ,σ∇ U(ũ(σ))dσ,

which implies

    ũ(x)=y_0 -h∫_0^xG^-1(ũ(τ))∫_0^1
    P_τ,σ∇
    V(ũ(σ))dσ dτ
    
    =   y_0 -h∫_0^xG^-1(ũ(τ))∫_0^1
    ∑_i=0^r-1ψ̃_i(τ)ψ̃_i(σ)∇
    V(ũ(σ))dσ dτ
    
    =   y_0 -h∑_i=0^r-1∫_0^xG^-1(ũ(τ))
    ψ̃_i(τ)dτ∫_0^1ψ̃_i(σ)∇
    V(ũ(σ))dσ.

Choosing x=c_j for j=2,…,r+1   and denoting
y_σ=ũ(σ), we obtain the following practical
FFED method

    {   ũ(σ)=∑_i=1^r+1y_c_il_i(σ),
       y_c_j
    = y_0 -h∑_i=0^r-1∫_0^c_jG^-1(ũ(τ))
    ψ̃_i(τ)dτ∫_0^1ψ̃_i(σ)∇
    V(ũ(σ))dσ,     j=2,…,r,
       y_1
    = y_0 -h∑_i=0^r-1∫_0^1G^-1(ũ(τ))
    ψ̃_i(τ)dτ∫_0^1ψ̃_i(σ)∇
    V(ũ(σ))dσ,
    .

which represents a nonlinear system of equations for the unknowns y_c_j for j=2,…,r+1 and it can be solved by iteration.


We note that the integrals appearing in the methods can be
calculated exactly for many cases.  If they cannot be directly
calculated, it is nature  to consider approximating them by a
quadrature rule.




§ NUMERICAL TEST

As an example, we choose 
    φ_k(t)=t^k,   k=0,1,⋯,r-1

for the function spaces X and Y,  and then take r=2 and
c_1=3-√(3)/6, c_2=3+√(3)/6 for our new
methods. We denote this method as FFED. In order to  show its
efficiency and robustness, we choose the following three methods in
the literature:


  * AVF:  the averaged vector field method studied
in <cit.>;



  * AVFC: the averaged vector field collocation method with s=2 given in
<cit.>;


  * EEI:  the explicit exponential integrator of order four derived in
<cit.>.


 It is noted  that  we approximate
the integrals appearing in the methods by the four-point
Gauss-Legendre's quadrature. For implicit methods, we
  set 10^-16 as the error tolerance and 10 as the maximum number of each fixed-point iteration.

Consider the   gradient system  (<ref>) with 
    G(y)=(
        [  cos(θ)  sin(θ); -sin(θ)  cos(θ);         ])^-1
 and

    U(y)=1/2r(y_1^2+y_2^2)-1/2sin(θ)(y_1y_2^2-1/3y_1^3)+1/2cos(θ)(-y_1^2y_2+1/3y_2^3),

where  θ=π/2-10^-4 and r=20. The  initial value is
chosen as y_1(0)=0, y_2(0)=1. This system has been researched in
<cit.>.   We first solve it in [0,100] with
h=1/10,1/20,1/100 and see Figure
<ref>   for the results of the potential function U(y).  Then
the system is solved   with h=0.1/2^i for  i=1,…,4 and
T=100,500,1000. The global errors  are presented in Figure
<ref>.


From these numerical results, it can be observed  that our FFED
method has a higher accuracy, a better energy-diminishing property,
and a more prominent damping behavior in comparison with the other
three methods.









§ CONCLUSIONS

In this letter we derived a novel kind of functionally fitted
energy-diminishing methods for solving gradient systems. The
properties of the methods have been analysed.  It was shown that the
 arbitrary-order  methods are unconditionally energy-diminishing and  achieve
damping for stiff gradient systems. We also discussed the
implementations of the methods. The remarkable efficiency of the
 methods was demonstrated  by
a numerical test in comparison with three existing numerical methods
in the literature.



§ ACKNOWLEDGEMENT

The research of the first author is supported in part by the
Alexander von Humboldt Foundation and by the  Natural Science
Foundation of Shandong Province (Outstanding Younth Foundation).

9
Bao04W. Bao, Q. Du,   Computing the ground state solution of
Bose-Einstein condensates by a normalized gradient flow, SIAM J.
Sci. Comput. 25 (2004) 1674-1697.



 Barrett02 J. Barrett, J. Blowey,  Finite
element approximation of an Allen-Cahn/Cahn-Hilliard system, IMA J.
Numer. Anal.  22 (2002) 11-71.


Chen02  L. Chen,  Phase-field
models for microstructure evolution, Annual review of materials
research, 32 (2002) 113-140.

Droske04 M.  Droske,    M. Rumpf,    A level set formulation for Willmore
flow, Interfaces Free Bound, 6 (2004) 361-378.

 Feng04 X. Feng,  A. Prohl,   Error analysis
of a mixed finite element method for the Cahn-Hilliard equation,
Numer. Math. 99 (2004) 47-84.


hairer2013 E. Hairer, C. Lubich,  Energy-diminishing integration of gradient
systems, IMA J. Numer. Anal.  34 (2014) 452-461.

hairer2006 E. Hairer, C. Lubich, G. Wanner,   Geometric Numerical
Integration: Structure-Preserving Algorithms for Ordinary
Differential Equations, 2nd edn.  Springer-Verlag, Berlin,
Heidelberg, 2006.

 Hairer93-stiff E. Hairer,  G. Wanner, Solving
Ordinary Differential Equations II. Stiff and Differential-Algebraic
Problems, Springer Series in Computational Mathematics 14, 2nd edn.
 Springer-Verlag, Berlin,
Heidelberg, 1996.

 Hochbruck2009 M. Hochbruck, A. Ostermann,     J.
Schweitzer,  Exponential rosenbrock-type methods,  SIAM J. Numer.
Anal.  47 (2009)  786-803.

Li_Wu(sci2016)  Y.W.   Li, X. Wu,   Exponential integrators preserving
first integrals or Lyapunov functions for conservative or
dissipative systems,  SIAM J. Sci. Comput. 38 (2016) 1876-1895.



 Li_Wu(na2016)   Y.W.   Li, X. Wu, Functionally fitted energy-preserving methods for
solving oscillatory nonlinear Hamiltonian systems,   SIAM J. Numer.
Anal.  54 (2016) 2036-2059.

 Liu_Wu(AML)   C. Li, X. Wu, The boundness of the operator-valued functions for multidimensional
nonlinear wave equations with applications, Appl. Math. Lett. 74
(2017) 60-67.

Lubich93 C. Lubich, A. Ostermann,
Runge-Kutta methods for parabolic equations and convolution
quadrature, Math. Comp.  60 (1993) 105-131.

Lubich96 C. Lubich, A. Ostermann,  Runge-Kutta time
discretization of reaction-diffusion and Navier-Stokes equations:
nonsmooth-data error estimates and applications to long-time
behaviour,  Appl. Numer. Math.  22 (1996) 279-292.

 Mclachlan-98
R. I. Mclachlan, G. R. W. Quispel, N. Robidoux,  A unified approach
to Hamiltonian systems, Poisson systems, gradient systems, and
systems with Lyapunov functions or first integrals, Phys. Rev. Lett.
81 (1998)  2399-2411

McLachlan99  R. I. McLachlan, G. R. W. Quispel,  N. Robidoux,  Geometric
integration using discrete gradient, Philos. Trans. R. Soc. Lond. A
357 (1999)   1021-1045.



Michailovich07O. Michailovich, Y.  Rathi, A. Tannenbaum,   Image
segmentation using active contours driven by the Bhattacharyya
gradient flow, IEEE Trans. Image Process. 16 (2007)  2787-2801.

Otto01 F. Otto,  The geometry of dissipative evolution equations: the
porous medium equation, Comm. Part. Diff. Equa. 26  (2001) 101-174.


Penzler12P. Penzler,   M.  Rumpf,  B. Wirth,    A phase-field model for
compliance shape optimization in nonlinear elasticity,  ESAIM
Control Optim. Calc. Var. 18 (2012) 229-258.

Strzodka04R. Strzodka, M.  Droske,    M. Rumpf,   Image registration by a
regularized gradient flow. A streaming implementation in DX9
graphics hardware,  Computing, 73   (2004) 373-389.



  wang-2016  B. Wang, A. Iserles, X. Wu, Arbitrary–order trigonometric  Fourier collocation methods for
multi-frequency oscillatory  systems,  Found. Comput. Math.  16
(2016) 151-181.


  17-newB. Wang, X. Wu,
Arbitrary-order exponential energy-preserving collocation methods
for solving conservative or dissipative systems, Preprint,
  (2017) arXiv:1712.07830

wu2017-JCAM  B. Wang,   X. Wu, F. Meng, Trigonometric collocation methods based
on Lagrange basis polynomials for multi-frequency oscillatory
second-order differential equations,
 J. Comput. Appl. Math. 313  (2017)   185-201.

 wang2017-JCM  B. Wang, X. Wu, F. Meng,  Y. Fang,  Exponential
Fourier collocation methods
 for solving first-order differential equations, J. Comput.
Math. 35 (2017)   711-736.


wu2013-book  X. Wu, X. You, B. Wang,  Structure-Preserving Algorithms for Oscillatory
Differential Equations, Springer-Verlag, Berlin, Heidelberg, 2013. ]