The Dynamical Functional Particle Method

# The Dynamical Functional Particle Method

Mårten Gulliksson, Sverker Edvardsson, and Andreas Lind Division of Computational Mathematics and Physics, Mid Sweden University, SE-85170 Sundsvall, Sweden
###### Abstract.

We present a new algorithm which is named the Dynamical Functional Particle Method, DFPM. It is based on the idea of formulating a finite dimensional damped dynamical system whose stationary points are the solution to the original equations. The resulting Hamiltonian dynamical system makes it possible to apply efficient symplectic integrators. Other attractive properties of DFPM are that it has an exponential convergence rate, automatically includes a sparse formulation and in many cases can solve nonlinear problems without any special treatment. We study the convergence and convergence rate of DFPM. It is shown that for the discretized symmetric eigenvalue problems the computational complexity is given by , where d is the dimension of the problem and N is the vector size. An illustrative example of this is made for the 2-dimensional Schrödinger equation. Comparisons are made with the standard numerical libraries ARPACK and LAPACK. The conjugated gradient method and shifted power method are tested as well. It is concluded that DFPM is both versatile and efficient.

###### Key words and phrases:
Dynamical systems, Linear eigenvalue problems, ARPACK, Particle methods, DFPM, Lyapunov function, Hamiltonian dynamics

## 1. Introduction

### 1.1. The Dynamical Functional Particle Method

The goal of this paper is to present an idea for solving equations by formulating a dynamical system whose stationary solution is the solution of the original equations. Examples of equations that can be solved are ordinary and partial differential equations, linear or nonlinear system of equations, and particularly eigenvalue problems. In this section we begin by formulating the equation and the dynamical system in an abstract setting. We then give the corresponding finite dimensional formulation by discretizing the infinite dimensional problem. This discretized problem will then be analyzed and studied throughout the paper.

Let be an operator and , where is a Banach space that will be defined by the actual problem setting. We consider the abstract equation

 (1.1) F(v)=0

that could be, e.g., a differential equation. Further, a parameter is introduced, interpreted as artificial time, which belongs to the interval . A related equation in is formulated as

 (1.2) μutt+ηut=F(u).

The parameters are the mass and damping parameters. The idea in the infinite dimensional setting is to solve (1.1) by solving (1.2) in such a way that when , i.e., In addition, the two initial conditions and are applied.

Both (1.1) and (1.2) need to be discretized to attain a numerical solution. For simplicity, we exemplify by applying finite differences but it is possible to use, e.g., finite elements, basis sets or any other method of discretization. We define a grid and approximate by and assume that the discretized version of (1.1) can be written as

 (1.3) Fi(v1…,vn)=0,i=1,…,n

where .

Turning to the dynamical system (1.2) it is discretized such that approximates and for . Further, is discretized as in (1.3) and we approximate (1.2) with the system of ordinary differential equations

 (1.4) μi¨ui+ηi˙ui=Fi(u1,…,un),i=1,…,n.

with initial conditions . Our idea in the discrete setting is to solve (1.3) by solving (1.4) such that when , i.e., The overall approach for solving (1.1) using (1.4) is named the Dynamical Functional Particle Method, DFPM.

### 1.2. Related work and topics

In a recent mechanics oriented article the connection between classical particle methods and differential equations were studied [19]. This work had a clear focus on the physical understanding and mechanical properties. The present work, however, turns the focus towards the mathematical aspects in the attempt to answer questions related to convergence, rate of convergence and the underlying reasons why DFPM is seen to be efficient for some mathematical problems. The idea of studying dynamical particle systems certainly has its origin in basic physics and astronomy. The assumption there is that all matter consists of particles. Their interactions are known and they follow the equations of motion. The basic idea DFPM, however, is that the “forces” and “particles” instead are viewed as mathematical abstract objects rather than physical. Thus mathematical quasi-particles are formed. Their interactions are determined by the functional equation at hand. From the mechanical point of view the quasi particles in (1.3) have masses and all follow a dissipated motion governed by . Such Hamiltonian systems have many interesting and useful properties, see, e.g., [21]. In Hamiltonian systems it is well known that symplectic integration techniques are especially attractive [30, 25].

The idea of solving a time dependent problem to get the stationary solution has also previously been applied in mathematics. A simple example is the solution of an elliptic PDE such as in the heat equation, see Sincovec and Madsen [32]. Indeed, steady state solutions are often the objective when simulating time-dependent PDE systems. Since the stationary state is seeked, the evolution of the system is considered to take place in artificial time. The concept of artificial time is further discussed and analyzed in [6].

A general approach is that of continuation, see [4] for an introduction, where (1.1) is embedded in a family of problems depending on a parameter , i.e.,

 (1.5) F(u;s)=0

where Thus, solving (1.5) for is equivalent to solving (1.1) and it is assumed that solving (1.5) for some say , is computationally cheap. The solution to (1.5) is found by solving a sequence of problems for values of decreasing from to . A general package for continuation methods with additional references may be found in Watson et al. [34]. Further, see Nocedal and Wright [27] for a discussion in the context of optimization and Ascher, Mattheij and Russell [7] for boundary value ODEs. DFPM can in principle be viewed as a sub-method to the group of continuation methods. However, as far as the authors know, the concrete application of a second order system (Hamiltonian dynamics) to solve equations and the corresponding analysis as presented here is novel.

Other works where (1.2) appear are for example the damped harmonic oscillator in classical mechanics, the damped wave equation, [28] and the heavy ball with friction [5]. These problem settings are specific examples of physical systems and not developed to solve equations in general.

In [13, 14] iterative processes to solve, e.g., eigenvalue problems are considered as (gradient driven) dynamical systems. So called fictitious time is used in [33] where, e.g., Dirichlet boundary value problem of quasilinear elliptic equation is numerically solved by using the concept of fictitious time integration method. The inverse problem of recovering a distributed parameter model is considered in [6] using the first order ODE attained from the necessary optimality conditions of the inverse problem.

First order systems, mainly in the form , have been used to solve different kinds of equations , both as a general approach and intended for specific mathematical problems. It is of interest to briefly consider the difference between the first order differential equation and the second order approach, DFPM. Suppose for simplicity that a discretization is made by finite differences leading to a system of equations . Consider an example where the functional contains a derivative w.r.t. u (A) and other functions of u and x (B). Then . Dimensional analysis then gives that

Given a certain component we have that is not variable, so . We see that the dimension of time is related to the discretization, i.e., . In a similar way it can be shown that for the second order differential equation (DFPM) one instead have that . In a numerical integration, the dimension of t must still be valid. This means that also the maximum timestep (for which a numerical algorithm is stable) is given by for the first order case, and for the second order equation. Consider for example the case of a central difference formula, i.e., p=2. For the first order differential equation we see that will have to be very small as finer meshes are selected. This will lead to a less efficient computational complexity for the first order system. The complexity of DFPM will be further discussed in section 6.

As we shall see, the DFPM algorithm seems attractive due to several reasons. The most interesting points that will be studied are related to computational complexity, easiness of implementation, Hamiltonian dynamics and its relation to the total evolution time, stability and cheapness of symplectic integration, exponential convergence and the existence of potential energy.

### 1.3. The outline of the paper

The outline of the paper is based on illustrating the versatility of DFPM and to analyze the convergence aspects of the dynamical system (1.4). In order to introduce the reader to DFPM, the damped harmonic oscillator is revisited. DFPM clearly has a close relationship to this type of classical system. It is then important to remind the reader of the close connection to Hamiltonian dynamics in particular the existence of a potential function. In such a case where the functional is conservative any extreme value of its potential function is a solution to the original problem (1.3). Specifically if the potential has a minimum the solution of (1.4) will converge asymptotically. This is dealt with in Section 4. A Lyapunov function is applied to show asymptotic convergence in Section 4.1. In Section 4.2 we analyze the linearization of DFPM and give precise statements for local asymptotic convergence valid close to the solution and for linear problems such as systems of linear equations. The rate of convergence is treated in Section 5 with four subsections treating the general problem (1.4) when there is a Lyapunov function, the linearized problem, the choice of damping, and examples, respectively. In Section 5.1 the Lyapunov function is used to state a general theorem that together with additional assumptions on the Lyapunov function gives an exponential convergence rate. This theorem is then specialized to the case when there exists a potential. The linearized problem is analyzed in Section 5.2 where we first treat the case with one scalar damping and then discuss the possibility to choose an optimal general damping matrix. The conclusions drawn from the choice of damping in the linear case are used in Section 5.3 to formulate a local strategy for the choice of optimal damping. To demonstrate the efficiency of DFPM we report in the end of the article several examples. The most noteworthy is the efficiency for treating symmetric eigenvalue problems. It is shown that DFPM is order of magnitudes faster than the standard software ARPACK [1]. Finally, in Section 7 we make some conclusions, discuss open problems as well as suggestions for future works.

## 2. Two illuminating examples

### 2.1. The damped harmonic oscillator

Despite its triviality an important example related to DFPM is the damped harmonic oscillator where which we include for later reference since it illustrates many properties of (1.4). In this case the equation at hand is given by

 (2.1) μ¨u+η˙u=−ku.

In DFPM, as well as here, we set the inital condition . The initial condition for may be set arbitrary. In mechanics the parameters , and correspond to particle mass, damping constant and spring constant, respectively. The time-dependent solution is given explicitly by

 u=c1e−γ1t+c2e−γ2t

where

 (2.2) γ1,2=−12ημ±√14η2μ2−kμ

and are constants given by the initial conditions. Although in mechanics it is clear that all parameters in (2.1) are physically positive, it is worthwhile to make a comment why this is so. Consider the case . The roots are then real with one positive and one negative root so will diverge. The situation for is similar with one positive real root and no convergence. When the roots may be complex but one root will always have a positive real part and the solution will not converge. Thus, the only possible choice for convergence into a stationary solution is to apply positive parameters.

There are three different regimes of the parameters that will effect the convergence: the under critical damping, which shows an oscillatory convergence, the critical damping, giving exponential convergence, and over critical damping, resulting in a slower exponential convergence. The critical damped system is known to be the fastest way for the system to return to its equilibrium (i.e., the stationary solution) [2]. It will be illustrative to return to this example later when considering various aspects of convergence and convergence rate of DFPM in Sections 4 and 5.

### 2.2. A symmetric eigenvalue problem

Symmetric eigenvalue problems are of great importance in physics and technology. It is also a relatively straight forward example to illustrate how the DFPM algorithm is applied. Consider the eigenvalue problem

 (2.3) Av=λv

where is a symmetric matrix with normalization . The DFPM equation (1.4) is not directly applicable since the eigenvalue is unknown. However, it is well known, see [22], that the stationary solutions of the Rayleigh quotient

 ρ(v)=vTAv,∥v∥=1

are eigenvectors of (2.3). Specifically, we have that

 argminρ(v)=λmin

where is the smallest eigenvalue of . Thus, one way to formulate the functional vector is

 F=−Au+(uTAu)u,∥u∥=1

where . The DFPM equation (1.4) is then given by

 (2.4)

where . This procedure will yield and its corresponding eigenvector . We shall see later that by replacing with we instead get and its corresponding eigenvector. There are various strategies to get the other solutions. One possibility is to apply a Gram-Schmidt process [11]. Often in applications, only a few of the lowest solutions are of interest. The reader should note that in practice the matrix never needs to be formulated explicitly. In DFPM one instead works with the components of the functional vector . The mechanical interpretation is of course that this is the force acting on particle i. The formulation therefore automatically becomes sparse. This is later illustrated in Section 6.

## 3. The Potential and Total Energy

DFPM can be considered as a many-particle system in classical mechanics where is the force acting on particle i, is its point mass, and is the damping force [21]. In this section we revisit the concept of a conservative force field and thus the existence of a many-particle potential. In DFPM the functional is not necessarily conservative, but if it is, the analysis is greatly simplified. By using the results in this section, we shall see in Section 4.1 that for a convex potential, the stationary solution to (1.4) corresponds to a minimum of the many-particle potential.

We start by taking the view that (1.3) is a vector field in :

 (3.1) F=(F1,F2,…,Fn),Fi=Fi(u1,u2,…,un).
###### Definition 1.

The vector field in (3.1) is conservative if there exists a potential function such that .

For any conservative field we have that

 F(v)=0⇔∇V(v)=0.

Thus, any solution to (1.3) is an extreme value of the potential , i.e., a minimum, maximum or saddle point. In other words, solving (1.3) by DFPM is equivalent to finding the extreme points of the potential . We will explore this fact further when analyzing the convergence of DFPM in Section 4.

By differentiating and assuming that is at least twice continuously differentiable, we get

 (3.2) ∂Fi∂uj−∂Fj∂ui=0,1≤i,j≤n.

as a necessary and sufficient condition for to be conservative. Note that one good example fulfilling the condition (3.2) is the force field where is a symmetric matrix. We shall see later that this fact is very useful for symmetric eigenproblems.

It is possible to derive the condition (3.2) that is interesting in its own since it contains the possibility to consider equations on manifolds. Consider the (work) -form , see [24] for a definition of -forms. The -form is said to be closed if and Poincarés Lemma [15] implies that any closed form is exact, i.e., in our context has a potential, say that is . There is no ambiguity to say that is a potential as in Definition 1. We have the following results for the vector field in (3.1).

###### Theorem 2.

The following statements are equivalent:

1. is conservative, that is

2.

3. is independent of the path

4. for all closed paths .

###### Proof.

The equivalence of 1 and 2 is trivial. The equivalence of 1 and 3 can be found in [24]. Now, assume that 3. is true. Let and be two arbitrary points, and let and be two piecewise smooth paths from to . Define as the closed path which first goes from to , via , and then from to via . Then, since , we get that

Since and are arbitrary points, and and are arbitrary, the closed path is arbitrary, and therefore the implication 3 to 4 is proved. The implication 4 to 3 is similar. ∎

If a potential exists it can be derived from the discretized equations by calculating the work, , simply integrating along any path, say, from to . For example it is possible to use coordinate directions as

 W = u1∫0F1(s1,0,…,0)ds1+u2∫0F2(u1,s2,0,…,0)ds2+… +un∫0Fn(u1,u2,…,un−1,sn)dsn=−V.

#### 3.0.1. A revisit to the symmetric eigenvalue problem

Recall that DFPM equation for the symmetric eigenvalue problem (2.4). The corresponding vector field

 F=−Au+(uTAu)u,∥u∥=1

is conservative with the potential

 (3.4) V(u)=12uTAu,∥u∥=1

To prove this it would at a first glance seem natural to find the gradient of . However, the normalization complicates this somewhat. This can be treated by investigate the gradient on the sphere . Denote the tangent space to the sphere at a point as . By using the Euclidean metric (the 2-norm), the gradient of at is the unique vector such that

 ∇V(u)Tt=∇Sn−1V(u)Tt

for all tangent vectors where is the usual gradient in . Solving this equation for by realizing that is the projection of on we get

 ∇Sn−1V(u)=(I−n(u)n(u)T)∇V(u)=∇V(u)−(n(u)T∇V(u))n(u)

where is the normal to . Since, for we have and we get

 ∇Sn−1V(u)=Au−(uTAu)u=−F(u)

showing that in fact is a potential to the vector field .

## 4. Asymptotic convergence

In this section we investigate the convergence properties of the solution given by (1.4). Since we are interested in the asymptotic solution we will use stability theory for dynamical systems, see [26], namely the use of a Lyapunov function and local linear analysis. However, it is generally difficult to find a Lyapunov function. Consequently, we start with the case where there exists a potential and where the Lyapunov function can be chosen as the total energy. If no potential exists we are left with the linear stability analysis in Section 4.2.

For simplicity and without loss of generality we may assume that and that the solution of (1.3) is .

### 4.1. Using the Potential

In this section we assume that there exists a potential, , to the given equations in (1.3). Then (1.4) may be written as

 (4.1) ¨u+N˙u=−∇V

where . The energy functional (the Lyapunov function) is given by

 (4.2) E=T+V

where

 T=12∑˙u2i

is the kinetic energy. We then have the following important result to be used in the analysis of the asymptotic convergence analysis.

###### Lemma 3.

Assume that . For the given solution of the dynamical system (1.4) the energy functional defined in (4.2) is a non-increasing function.

###### Proof.

From the definition of in (4.2) and we have

 dEdt=dTdt+dVdt=∑˙ui¨ui+∂V∂ui˙ui=∑˙ui¨ui−Fi˙ui

and since we get

 (4.3) dEdt=−∑iηi˙u2i≤0.

Lemma 3 tells us that the energy is non-increasing which is not surprising from a mechanical point of view since the damping will decrease the total amount of energy and there are no additional sources of energy in the system.

The next theorem is taken from [26]. The proof is omitted.

###### Theorem 4.

Consider the autonomous system

 (4.4) ˙w=G(w)

where is a continuous function defined on a domain in containing the origin and . Assume that the Lyapunov function is non-negative with respect to (4.4) for all and such that for some constant the set is a closed and bounded component of the set . Let be the largest invariant set in the set

 Z={w∈Ω:dLdt(w)=0}.

Then every solution of the autonomous system (4.4) with approaches the set as .

Using Theorem 4 we are now able to show our main result in this section for the asymptotic convergence of (4.1).

###### Theorem 5.

Assume that there exists a solution, , to (1.3) and a potential to . If is globally convex, i.e., convex in the whole of , then for any initial starting point in the solution of (4.1) will be asymptotically convergent to . If is locally convex in a neighbourhood of then for any initial starting point in the solution of (4.1) will be asymptotically convergent to .

###### Proof.

Rewrite DFPM (4.1) as a system by letting as

 (4.5) [˙u˙v]=[v−Nv−∇uV]

We want to use Theorem 4 to prove asymptotic convergence of DFPM. From Lemma 3 we have that and therefore if and only if Define

 w=[uv].

Let be the invariant set in Theorem 4. Then, by Theorem 4 again, the solution to the system (4.5) approaches the set as If then , so from (4.5) we have that if and can not remain in the set if . We need to verify that as but this follows from the fact that is non-increasing. Hence as

### 4.2. Linear Convergence Analysis

Without a Lyapunov function and with no existing potential one is left with a local linear stability analysis. Such an analysis is based on the linearization of the dynamical system (1.4) at a solution , , to (1.3). Define as the Jacobian of . From the Taylor expansion we define the linearized problem to (1.4) as

 (4.6) ^M¨u+^N˙u=F(0)+J(0)u

where

 ^M=diag(^μ1,…,^μm),^N=diag(^η1,…,^ηm)

and are the values of at . Thus, the local convergence can be analyzed by analyzing the linear system (4.6) which for notational convenience is written

 (4.7) M¨u+N˙u+Au=b

where . In [17] sufficient conditions are given for (4.7) to have asymptotically stable solutions. In order to state these conditions we need some additional notation. Consider the homogeneous equation, i.e.,

 (4.8) M¨u+N˙u+Au=0.

By inserting the eigensolution into (4.8) we get the equation

 (4.9) (ξ2iM+ξiN+A)vi=0

for the eigenvalue and eigenvector of . Let , denote the real and imaginary part of , respectively. Introduce the symmetric and anti-symmetric parts of as . and define

 si(M)=Re(vi)TMSRe(vi)+Im(vi)TMSIm(vi),ai(M)=2Re(vi)TMAIm(vi)

with corresponding definitions for . By applying the general Hurwitz criterion, see e.g. [20], to (4.9) we get the following theorem and corollary. For a detailed presentation of the proofs we refer to [17].

###### Theorem 6.

The solution to (4.7) will converge asymptotically if and only if

 si(M)si(N)+ai(M)ai(N)>0

and

 (si(N)si(A)+ai(N)ai(A))(si(M)si(N)+ai(M)ai(N)) − −(si(M)ai(A)−ai(M)si(A))2 >0
###### Corollary 7.

If are positive definite and has eigenvalues with positive real parts then the solution to (4.7) will converge asymptotically.

Let us now return to the question of local convergence for (1.4) and state the following theorem.

###### Theorem 8.

Define as the Jacobian of . Assume that there exists a solution, , to (1.3). Further, assume that

 ^M=diag(^μ1,…,^μm),^N=diag(^η1,…,^ηm)

where are the values of at . Then DFPM will converge asymptotically for any initial starting point of (1.4) close enough to if and only if the conditions in Theorem 6 are fulfilled where . Further, if are positive definite and has eigenvalues with negative real parts then DFPM will converge asymptotically for any initial starting point of (1.4) close enough to .

###### Proof.

The first statement in the theorem follows directly from Theorem 6 and the second from Corollary 7. ∎

## 5. Convergence rate

### 5.1. General results on convergence rate

Sharp estimates of the convergence rate for DFPM in a general case is difficult and not realistic. However, we shall give some important special cases that is relevant for solving equations with DFPM, i.e., to achieve fast exponential convergence. We emphasize that this is crucial to attain an efficient method for solving (1.3). We again assume without loss of generality that the solution of (1.3) is .

###### Definition 9.

The solution, , to (1.3) is locally exponentially stable and the solution to (1.4) has a local exponential convergence rate if there exists an and for every and every , there exists a such that for all solutions of (1.4) for all whenever . The solution, , to (1.3) is globally exponentially stable and the solution to (1.4) has a global exponential convergence rate if for there exists such that whenever .

We begin by stating a general theorem from [26] giving one possible formulation of the requirements for exponential convergence based on the existence of a Lyapunov function .

###### Theorem 10.

Assume that there exist a Lyapunov function , and four positive constants , such that

 c1∥∥∥[u˙u]∥∥∥p≤L(u,˙u)≤c2∥∥∥[u˙u]∥∥∥p

and

 dLdt≤−c3∥∥∥[u˙u]∥∥∥p

for all . Then the solution, , to (1.3) is locally exponentially stable and the solution to (1.4) has a local exponential convergence rate. If is the whole the solution to (1.3) is globally exponentially stable and the solution to (1.4) has a global exponential convergence rate.

In the case that there exists a potential we can choose the Lyapunov function as the total energy in (4.2). We will state a theorem that proves exponential convergence in this case that is slightly different from Theorem 10.

###### Theorem 11.

Assume that there exists a potential which is locally convex in a neighbourhood of the solution, , to (1.3) and satisfies the bound

 (5.1) c1∥u∥p≤V(u)≤c2∥u∥p

where and are positive constants. Then is locally exponentially stable and the solution to (1.4) has a local exponential convergence rate. If the potential is convex and satisfies the bound (5.1) in the whole the solution to (1.4) has a global exponential convergence rate.

###### Proof.

From (5.1) and the definition of kinetic energy we have that the total energy is bounded as

 (5.2) c1∥u∥p+12μmax∥˙u∥2≤E≤c2∥u∥p+12μmax∥˙u∥2

Thus, from Lemma 3 we have

 dEdtE≤−ηmax∥˙u∥2c2∥u∥p+12μmax∥˙u∥2

If is convex we see from Theorem 5 that unless at the solution and therefore there exists a constant such that

 dEdt≤−γE

and thus . From (5.2) we see that

 c1∥u∥p+12μmax∥˙u∥2≤E≤E0e−γ(t−t0)

and this implies

 ∥u∥≤ce−1pγ(t−t0)

for some positive constant depending only on the initial conditions. ∎

### 5.2. Convergence rate for linear problems

Consider again the damped harmonic oscillator (2.1). Obviously, the convergence rate is exponential for all different dampings. However, the fastest convergence rate is achieved for the case where the damping is chosen as critical damping which can be seen in (2.2): The negative real part below critical damping is and therefore should be as large as possible. However, as soon as critical damping is exceeded, one of the roots will be real and the negative real part will increase for the larger real root.

This property is inherited for more general linear problems defined by (4.7) which we now shall investigate further. We will assume that and are diagonal matrices with positive elements and then we can, without loss of generality, assume that and consider the system

 (5.3) ¨u+N˙u+Au=0,N=diag(η1,…,ηn)

where . This linear system can be restated as

 [˙u˙v]=[0I−A−N][uv]

and since this is a linear autonomous system of first order differential equations, it is well known that any convergent solution will have exponential convergence rate. However, we are interested in solving the linear equations as fast as possible and then it is reasonable to try to answer the question: How fast is the exponential rate of convergence? This is a difficult question for a general and because of the following result from linear system theory [29].

###### Theorem 12.

Any two square matrices that are diagonalizable have the same similarity transformation if and only if they commute.

For the problem (5.3), Theorem 12 means that is a necessary and sufficient condition for having a similarity transformation such that where are diagonal matrices. In other words, the two matrices and has to commute in order to decouple the system (5.3) into one dimensional damped harmonic oscillators where the optimal damping is critical damping as shown earlier. Note that the special case where all damping parameters are the same, trivially commutes with any matrix .

### 5.3. A discussion of optimal damping

The problem of choosing the damping in (1.4) such that the asymptotic solution is attained, to some precision, in minimal time is difficult, see, e.g., [12],[31] for some special cases. Moreover, from a practical point of view this is not very interesting since it requires a priori knowledge of the solution. A more interesting approach is to choose the damping according to a local measure of curvature which we will discuss briefly. Assume that the solution to (1.3) is and that there exists a potential that is convex with and a positive definite Hessian . Consider the case with a single damping parameter . Then a Taylor expansion at in (1.4) gives the approximate problem

 (5.4) ¨u+η(t)˙u=−∇2V(0)u

From the linear case treated in Section 4.2 the optimal damping for (5.4) is where denotes the smallest positive eigenvalue. Now, consider any then we conjecture that a good choice of damping is . To illustrate the possibilities of this choice of damping we given an example with a potential

 (5.5) V=eu21+2u22

that is globally convex with a minimum at .

In Figure 5.1 the eigenvalue distribution is shown for the smallest of the eigenvalues of . Indeed, looking at the trajectories in Figure 5.2 for the choice (solid line) and (curve indicated with ’*’) it is clearly seen how the choice of damping affects the convergence. In fact, the effect is rather striking and further analysis and tests of our conjecture is of great interest in order to improve DFPM.

## 6. DFPM example for the Helium atom

A relevant numerical application is to study the s-limit case of the Helium ground state energy. This example is often used in atomic physics literature as a benchmark to study numerical accuracy and efficiency. The equation at hand is the Schrödinger equation. Due to electronic correlation and consequently discontinuities (“Cato cusps”) this is often considered to be a tough problem. Another complication is the many-particle character leading to extremely high dimensionality. The Helium example here only slightly touches these problems because the full correlation term has been neglected, i.e., that term is replaced by resulting in only a 2D problem. This example is nevertherless sufficient to demonstrate many of the properties of DFPM. More complex examples have already been tested and DFPM remains relevant. However these fall outside the scope of the present work.

Accordingly, consider the Schrödinger equation for the s-limit case of Helium [18]:

 (6.1)

The boundary conditions are given by . The discretization can be made by using central finite differences with equidistant mesh sizes (for both and ) where is an integer chosen to get different problem sizes. The discretized version of of a certain particle at the position becomes

 ^Hv(r1i,r2j)≈Huij=−12ui−1,j+ui+1,j+ui,j−1+ui,j+1−4uijh2−−2uijr1i−2uijr2j+uijmax(r1i,r2j)

Note that the matrix is never explicitly needed so the formulation is automatically sparse. The interaction functional component, i.e., the “force” acting on the particle at the position is given by . This can be compared with the equation (2.4) derived earlier. The notation is the trapezoidal approximation to . The required norm is in the present context given by . The DFPM equation 1.4 for particle is thus given by

 (6.2) uij−Huij=μ¨uij+η˙uij,=1

In this case we apply constant mass and damping parameters ( and ). The boundary is set to which is sufficient for accurate ground state results. The integration method used to solve the ODE (6.2) is the symplectic Euler method, see e.g. [23]. The related Störmer-Verlet method was tested as well but gave no performance advantage. The test results are tabulated in Table 6.1.

The first column shows k which determines the discretization h as mentioned above. In the second column, the corresponding total number of particles, , is listed. Only a triangular domain needs to be computed due to even symmetry of the solution (i.e., ). Then the third column contains the maximum timestep used in the Symplectic Euler method (depends on h). In the fourth column the eigenvalues to 13 significant figures are listed. These values can easily be extrapolated to continuum (i.e., ). This extrapolated value is listed in the last line. In the final three columns we list the total CPU times in seconds to complete the computations to the desired accuracy. DFPM and three other methods are compared.

A single C-code was written where the only difference was whether a function call was made to DFPM, DACG, ARPACK or S-Power. The DACG method (Deflation Accelerated Conjugated Gradient) is an advanced method to find some of the lower eigenvalues of a symmetric positive definite matrix. The algorithm is described in refs. [10, 9]. Unfortunately the DACG algorithm cannot be used directly because there are both negative and positive eigenvalues present in the eigenvalue problem (6.1). Consequently, a shift of the potential (diagonal elements) had to be done. The shift for best performance was identified to be 3.9. The shifted power method (i.e., ’S-Power’) is a simple method that converges to the dominant eigenvalue [3]. This shift was also adjusted to get the maximum performance possible (about ). In the case of ARPACK, it is available at [1]. This iterative numerical library is based on the Implicitly Restarted Lanczos Method (IRLM). All tests were performed on a Linux PC with 1 GB primary memory and the CPU was a AMD Sempron 3600+ 2.0GHz. The compilers used were gcc and g77 version 3.4.4. Both the numerical library and C-code were compiled using the optimization: ’-O’. The CPU times were measured using the Linux command ’time’ three times (the best reading is listed in Table 6.1).

As can be seen in Table 6.1, the performance of DFPM is quite good considering that no real optimization of DFPM has been attempted yet. DACG performs in parity with DFPM, although it is some 40% less efficient for this example. Both ARPACK and S-Power are far behind. However, ARPACK is seen to be more than twice as efficient as the S-Power method. For the smaller problems in Table 6.1, DFPM is one order of magnitude faster than ARPACK. For the larger problems, the advantage to DFPM is seen to approach two orders of magnitudes. Also LAPACK was put to test (routine DSBEVX). This routine computes selected eigenvalues and, optionally, eigenvectors of a real symmetric band matrix. The parameter ’jobz’ was set to compute eigenvalues only and the range was set to compute only the lowest eigenvalue. Unfortunately, due to how DSBEVX handles matrix memory allocation, large sizes, as applied here, is not realistic to compute. Even the smallest size N in Table 6.1 requires to complete.

There are several reasons for good efficiency of DFPM. Firstly, the symplectic integration method for the second order differential equation allows a relatively large timestep. High numerical accuracy is not necessary during the evolution towards the stationary state. Only stability and speed is desired. Secondly, the cost per iteration is small because the symplectic integration method is as fast as Euler’s method. Due to the damped dynamical system, the convergence rate is exponential in time. This is also theoretically consistent with Section 5.

Further, in Fig. 6.1 it is seen that the computational complexity of ARPACK is given by . The DFPM complexity, however, is found to be . The cost of one iteration is the same for all the methods. Because of the sparsity of , it is . The number of iterations is what separates the various algorithms. For ARPACK and S-Power the number of iterations are both . In fact, in the case of S-Power it is straight forward to prove that the number of iterations for a d-dimensional problem is given by , i.e., for . DFPM and DACG, however, both show complexity for the number of iterations. It is interesting to discuss this further for the DFPM algorithm.

Consider equation (6.2) and let us assume that after only a few iterations, . Most of the iterations are then spent solving:

 E0uij−Huij=μ¨uij+η˙uij

for each one of the particles . The symplectic integration method is applied to approximate on its way to the stationary solution. By applying a linear stability analysis of the symplectic method for the problem at hand, one finds that the maximum step size is

 △tmax=2√E1−E0+√EN−1−E0

where is the dominant eigenvalue. For the present problem, using the central difference formula, it is easy to show that this eigenvalue depends on the mesh size h according to, . Since the mesh size only constitutes a minor correction to the lowest eigenvalues and , we have that . The equidistant mesh size h and the total number of particles N in a d-dimensional problem are related, i.e. , thus . The number of iterations is therefore given by , where is the total evolution time until the stationary solution is achieved. The time does not depend on N. That is, it is assumed that the mesh is fine enough and that the symplectic Euler algorithm approximately follows the true solution curve (i.e., the continuum case). The total complexity for DFPM is then given by . In the present test case, d=2, and the behavior in Fig. 6.1 is thus explained.

The presented benchmark indicates that DFPM is some 40% more efficient than DACG. DFPM can be further optimized by 1. allowing a varying timestep during iterations, 2. preconditioning of the functional F and 3. allowing variation in the damping parameter during the iterations. However, also DACG can be optimized by a preconditioning of the matrix H. In the benchmark tests we experience that DFPM is more robust than DACG. DACG is quite sensitive to the selected shift. If the shift is small it has tendencies to diverge (despite that all eigenvalues are positive). If the shift is larger the convergence rate starts to suffer. The selected shift is which gave the best convergence rate. Since DACG requires that the matrix H is positive definite, it would seem that DFPM is a better choice for Schrödinger problems where there often is a mix of positive and negative eigenvalues. Initially, the lowest (negative) eigenvalue is of course unknown which is why it can be difficult to apply DACG since one cannot assume a priori to know an appropriate shift. Another advantage is that the DFPM algorithm (1.4) remains the same also for nonlinear problems. The idea presented here is therefore quite versatile. The basic idea of DFPM as a dynamical system may also be attractive from a user’s perspective since the algorithm is physically intuitive and very pedagogical.

## 7. Conclusions and Future Work

We have presented a new versatile method to solve equations that we call the Dynamical Functional Particle method, DFPM. It is based on the idea of formulating the equations as a finite dimensional damped dynamical system where the stationary points are the solution to the equations. The attractive properties of the method is that it has an exponential convergence rate, includes a sparse formulation for linear systems of equations and linear eigenvalue problems, and does not require an explicit solver for nonlinear problems.

There are still a lot of interesting questions to be addressed. This includes the details for solving the ODE (1.4) in the most stable and efficient way. Motivated by the numerical tests reported in Section 6 we believe that symplectic solvers, see [8], will be of great importance here. However, the stability properties of the ODE solver is linked to the choice of parameters and especially the damping parameter. The key question is how to find the maximum time step retaining stability and getting a fast exponential convergence rate. We are currently working on these issues using the ideas presented here.

DFPM has proved useful for very large sparse linear eigenvalue problems as indicated in Section 6. It is plausible that DFPM will also be efficient for very large and sparse linear system of equations.

Since DFPM has a local exponential convergence rate it may be that it can be an alternative to the standard methods for nonlinear system of equations such as quasi-Newton or trust-region methods, see [16]. Moreover, if there exists a potential (or a Lyapunov function) DFPM has global convergence properties that can be useful for developing a solver for nonlinear problems with multiple solutions.

## References

• [1]
• [2]
• [3]
• [4] E. L. Allgower and K. Georg. Numerical Continuation Methods: An Introduction. Springer, New York, 1990.
• [5] Felipe Alvarez. On the minimization property of a second order dissipative system in Hilbert spaces. SIAM J. Control Optim., 38:1102–1119, 2000.
• [6] U. Ascher, H. Huang, and K. van den Doel. Artificial time integration. BIT, 47:3–25, 2007.
• [7] U. Ascher, R. Mattheij, and R. Russell. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. SIAM, 1995.
• [8] U. M. Ascher. Numerical Methods for Evolutionary Differential Equations. SIAM, 2008.
• [9] L. Bergamaschi, G. Pini, and F. Sartoretto. Approximate inverse preconditioning in the paralell solution of sparse eigenproblems. Numer. Linear Algebra Appl., 7:99, 2000.
• [10] L. Bergamaschi and M. Putti. Numerical comparison of iterative eigensolvers for large sparse symmetric positive definite matrices. Comput. Methods Appl. Mech. Engrg., 191:5233, 2002.
• [11] Å. Björck. Numerics of Gram-Schmidt orthogonalization. Linear Algebra Appl., 197–198:297, 1994.
• [12] C. Castro and S. J. Cox. Achieving arbitrarily large decay in the damped wave equation. SIAM J. Control Optim, 39(3):1748–1755, 2001.
• [13] Moody T. Chu. On the continuous realization of iterative processes. SIAM Review, 30(3):375–387, 1988.
• [14] Moody T. Chu. Numerical linear algebra algorithms as dynamical systems. Acta Numerica, pages 1–86, 2008.
• [15] L. Conlon. Differentiable Manifolds: A first course. Birkhauser, 1993.
• [16] P. Deuflhard. Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms. Springer, New York, 2004.
• [17] A. M. Diwekar and R. K. Yedavalli. Stability of matrix second-order systems: New conditions and perspectives. IEEE Transactions on Automatic Control, 44(9), 1999.
• [18] S. Edvardsson, D. Aberg, and P. Uddholm. Comp. Phys. Commun., 165:260, 2005.
• [19] S. Edvardsson, M. Gulliksson, and J. Persson. The dynamical functional particle method: An approach for boundary value problems. In press. Journal of Applied Mechanics, 2011.
• [20] F. R. Gantmacher. Applications of Theory of Matrices. New York, Wiley, 1959.
• [21] H. Goldstein. Classical Mechanics, 2nd ed. Addison-Wesley Publishing Company, 1980.
• [22] G. Golub and van Loan. Matrix Computations. John Hopkins Univesity Press, 1996.
• [23] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration, 2nd ed. Springer, 2006.
• [24] John Hubbard and Barbara Burke Hubbard. Vector Calculus, Linear Algebra, and Differential Forms: A Unified Approach. Matrix Editions, 2009.
• [25] B. Leimkuhler and S. Reich. Simulating hamiltonian dynamics. Cambridge university press, 2004.
• [26] A. N. Michel, L. Hou, and D. Liu. Stability of Dynamical Systems. Birkhauser, 2008.
• [27] J. Nocedal and S. Wright. Numerical Optimization. Springer, New York, 1999.
• [28] Vittorino Pata and Marco Squassina. On the strongly damped wave equation. Commun. Math. Phys., 253:511–533, 2005.
• [29] A. Srikantha Phani. On the necessary and sufficient conditions for the existence of classical normal modes in damped lineary dynamical systems. Journal of Sound and Vibration, 264:741–743, 2003.
• [30] J. M. Sanz-Serna and M. P. Calvo. Numerical Hamiltonian Problems. Chapman & Hall, 1994.
• [31] S. M. Shahruz, G. Langari, and M. Tomizuka. Optimal damping ratio for linear second-order systems. JOTA, 73(3):564–576, 1992.
• [32] R. Sincovec and N. Madsen. Software for nonlinear partial differential equations. ACM Trans. Math. Softw., 1:232–260, 1975.
• [33] Chia-Cheng Tsai, Chein-Shan Liu, and Wei-Chung Yeih. Fictious time integration mthod of fundamental solutions with chebyshev polynomials for solving Poisson-type nonlinear pdes. CMES, 56(2):131–151, 2010.
• [34] L. T. Watson, M. Sosonkina, R. C. Melville, A. P. Morgan, and H. F. Walker. Alg 777: Hompack90: A suite of fortran 90 codes for globally convergent homotopy algorithms. ACM Trans. Math. Softw., 23:514–549, 1997.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters