Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.mccme.ru/~ansobol/otarie/slides/zheligovsky.pdf Äàòà èçìåíåíèÿ: Fri May 21 21:15:08 2010 Äàòà èíäåêñèðîâàíèÿ: Sat Jun 26 04:49:01 2010 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: massive stars

` THE MONGE-AMPERE EQUATION: NOVEL FORMS AND NUMERICAL SOLUTION
V.A. Zheligovsky (vlad@mitp.ru), O.M. Po dvigina International Institute of Earthquake Prediction Theory and Mathematical Geophysics, 84/32 Profsoyuznaya St., Moscow, Russian Federation; U. Frisch Observatoire de la C^ d'Azur, ote U.M.R. 6529, BP 4229, 06304 Nice Cedex 4, France

Zheligovsky V., Podvigina O., Frisch U. The Monge­Amp`re equation: e various forms and numerical methods. J. Computational Physics, 229, 2010, 5043-5061 [http://arxiv.org/abs/0910.1301]


Monge-Amp`re equation: det u e
· · · · ·

xi ,xj

= f (x)

The "2nd order divergence" form The Fourier integral form The convolution integral form A test problem with a cosmological flavour Numerical solution: ­ Fixed point algorithm for the regular part of the MAE ­ Algorithm with continuation in a parameter and discrepancy minimisation ­ Algorithm with improvement of convexity and discrepancy minimisation · Results of solution of the test problem


Existence and regularity of solutions to the MAE: · A.V. Pogorelov. The Minkowski multidimensional problem. Nauka, Moscow, 1975. (Geometric approach, convexity.) · I.J. Bakelman. Convex analysis and nonlinear geometric elliptic equations. Springer-Verlag, 1994. · L.A. Caffarelli, X. Cabr´. Fully nonlinear elliptic equations. American Mathematical e Society colloquium publications, vol. 43. Amer. Math. Soc., Providence, Rhode Island, 1995. (The PDE approach, viscous solutions.) Methods for discrete optimal transportation problem: ask Andrei Sobolevsky. Numerical algorithms linked to the geometric interpretation of the MAE: 2 2z 2z 2z · V.I. Oliker, L.D. Prussner. On the numerical solution of the equation 2 2 - =f x y x y and its discretizations, I. Numerische Mathematik, 54 (1988) 271­293. (A mesh comprised of 25 points.) · D. Michaelis, S. Kudaev, R. Steinkopf, A. Gebhardt, P. Schreiber, A. Br¨uer. Incoherent a beam shaping with freeform mirror. Nonimaging optics and efficient illumination systems V. Eds. R. Winston, R.J. Koshel. Proc. of SPIE, vol. 7059 (2008) 705905. (A 55 â 55 mesh, 15 min. of a Pentium 4.)


Application of algorithms for sadd le-point optimisation to the two-dimensional MAE: · J.-D. Benamou, Y. Brenier. A computational fluid mechanics solution to the Monge­ Kantorovich mass transfer problem. Numerische Mathematik, 84 (2000) 375­393. · E.J. Dean, R. Glowinski. Numerical solution of the two-dimensional elliptic Monge­Amp`re e equation with Dirichlet boundary conditions: an augmented Lagrangian approach. C. R. Acad. Sci. Paris, Ser. I, 336 (2003) 779­784. Various numerical approaches: · E.J. Dean, R. Glowinski. Numerical solution of the two-dimensional elliptic Monge­Amp`re e equation with Dirichlet boundary conditions: a least-squares approach. C. R. Acad. Sci. Paris, Ser. I, 339 (2004) 887­892. · E.J. Dean, R. Glowinski. Numerical methods for fully nonlinear elliptic equations of the Monge­Amp`re type. Comput. Methods Appl. Mech. Engrg. 195 (2006) 1344­1386. e · X. Feng, M. Neilan. Galerkin methods for the fully nonlinear Monge­Amp`re equation e (arXiv:0712.1240). · X. Feng, M. Neilan. Mixed finite element methods for the fully nonlinear Monge­Amp`re e equation based on the vanishing moment method (arXiv:0712.1241). · G. Loeper, F. Rapetti. Numerical solution of the Monge­Amp`re equation by a Newton's e algorithm. C. R. Acad. Sci. Paris, Ser. I, 340 (2005) 319­324. (A pseudospectral Newton's algorithm.)


· J.-D. Benamou, B.D. Froese, A.M. Oberman. Two numerical methods for the elliptic Monge­Amp`re equation. Preprint, 2009 [www.divbyzero.ca/froese/w/images/4/40/MA.pdf]. e u=
-2

u2 x

1 x1

+ u2 x

2 x2

+2u

2 x1 x

2

+2f

An iterative Newton­Krylov solver with preconditioning (finite differences, modest accuracy required, discrepancies the order of 10-3 - 10-4 acceptable): · G.L. Delzanno, L. Chac´n, J.M. Finn, Y. Chung, G. Lapenta. An optimal robust equidiso tribution method for two-dimensional grid adaptation based on Monge­Kantorovich optimization. J. Comput. Physics, 227 (2008) 9841­9864. (256 â 256 grid, contrast ratio= 8886, discrepancy=7.78 â 10-5 , 70 s. of a 2.4 GHz Intel Xeon processor.) · J.M. Finn, G.L. Delzanno, L. Chacon. Grid generation and adaptation by Monge­Kantorovich optimization in two and three dimensions. Proceedings of the 17th International Meshing Roundtable (2008) 551­568.


The "2nd order divergence" form
A Fourier integral solution in RN : u = N 1 det aij i ...i j ...j ai j N ! i1 ,...,iN , 1 N 1 N n=1 n
j1 ,...,jN

R
n

N

u ( ) e i ·x d .
1

(p

...pN

is the unit antisymmetric tensor of rank N )
N



det ux (-1) = N!
N

i xj

=

(-1) N!

N


i1 ,...,iN , j1 ,...,jN

i1 ...iN

j
N

1

...j

N

n=1 R

N

u( )in jn ei ·x d
N

R

N

...
N n=1 N

R

N



i
i1 ,...,i
N

1

...i

N

n=1

n in



j
j1 ,...,j
N

1

...j

N

n=1

n jn



â



(-1) = N! â ( 1 ,...,
N N -1 n=1

u( n ) exp i
R
N



N n=1

n · x d 1 ... d
N -1

N

...

is the N â N matrix comprised of columnar vectors 1 ,..., N )

u( n ) u -



R

N

det2 1 ,...,
N -1 n=1

, -

N -1 n=1



n



n i e



·x d 1 ... d

N -1

d


(-1) = N!

N R
N

...

R

N

det2 1 .., ,.

N -1

,



N -1 n=1

u( n ) u -





N -1 n=1

nei ·x d 1 ... d



N -1

d .

"Reverse engineering" yields the "second-order divergence" form of the MAE in RN :

1 N!

i1 ,...,iN ,j1 ,...,jN

i

1 ... iN



j1 ... jN

u

x i1 x j 1

... u

xiN - 1 xj N - 1

u

xiN xj N

= f.

If u is space-periodic, and is a smooth function with a finite support, 1 N!
i1 ,...,iN , j1 ,...,jN

i1 ... iN j1 ... jN RN

u

xi1 xj 1

... u

xi N - 1 xj N - 1

u

xi N x j N

dx =

R3

f d x.

2 u WN -1 (T N ) the integrals are ::::::::::::::::: well-defined :::: (by the Sobolev embedding theorem, u L2(

N -1)

(T N ) u L (T N )).

By contrast, integrals in the similar identity obtained by one integration by parts 2 are not well-defined for u WN -1 (T N ). ::::::::::::::::::::::::::


The Fourier integral form of the MAE
For a space-periodic solution, 0 =

To accommodate f > 0 (of interest in cosmology), let c2 u = |x| + u , u = 0. 2
· denotes the average: f = lim
L

T3

f dx.

1 (2L)

3 [-L,L]

3

f (x ) d x .

cN = f ;
2

f/cN =

RN

f ( )ei ·x d .
RN

u=

RN

( )ei ·x d ,

u=

u ( )ei ·x d

u ( ) = - ( )/ | | 2 .


det u â
m

(ia is a unit vector in the direction of a). Consider the term of order m < N in u : (-1) N! (the sum symbol) i
i1 ,...,iN ,j1 ,...,j | | = m
N 1

xi xj N -1 n=1

=

1 N!

( n ) -


R

N

...

R

N

det2 i 1 ,..., i
n i e

N -1

, i

-

N -1 n=1



n

N -1 n=1

·x d 1 ... d

N -1

d

...i

N

j1 ...j

N

R | |=m n: in ,jn

N

u ( )in jn ei ·x d




n: i
n n

in j

n

or jn /

is over all subsets {1,..., N } of cardinality m;
m 1p1 <... N -m

in j

is the Kronecker
2 n jn

(-1) = m!

(
R

N n )

R

N

...


R

N



m

j
j1 ,...,j n
m

1

...jm p1 ...pN

-m


n=1

m

exp i

n=1

· x d 1 ... d
m-1



m

= â
m- 1 n=1

R

...

N

A


m

i 1 ,..., i
m-1 n=1

, i

( n ) -



n i e



-

m- 1 n=1



n

·x d 1 ... d

m- 1

d ,


where Am (i1 ,..., im )

1 m!

M
1p1 <... -m

N

2 p1 ...pN

-m

(i1 ,..., im )

is the sum of squares of all minors of rank m, M
p1 ...pN
-m

(i1 ,..., im )

j
j1 ,...,j
m

1

...jm p1 ...p

N -m

(i1 )j1 ... (im )jm ,

obtained by crossing out rows of numbers p1 < ... < p Mm i1 ,..., i comprised of m columnar vectors i1 ,..., im .
m

N -m

from the N â m matrix

,

The Fourier integral form of the MAE:
N

( )+ â
m- 1 n=1

( n ) -

m=2 R

N

...

R

N

A

m

i 1 ,..., i d 1 ... d

m-1

, i

-

m- 1 n=1



n

m- 1 n=1



n

m- 1

= f ( ),

= 0 .


Sharp bounds for the kernels Am in the Fourier integral form: 0 Am(i1,..., im) 1 m!

By the Gram­Schmidt orthogonalisation process, change is js such that (i) j1 = i1 , (ii) s, js differs from is by a linear combination of is for s < s ( Am is unaltered), (iii) s, js is orthogonal to all js for s < s |js | = sin s , where s is the angle between is and the subspace spanned by {is | s < s}. Denote vs = ijs . Denote Mm v1 ,..., v det aij =
j1 ,...,jm m m

Am (i ,..., i ) = Am ( v ,..., v )
s=2

1

m

1

m

m

sin2

s

and t Mm the transpose of Mm . a
ij
i



j1 ...j

m

i=1



M
1p1 <... -m

N

2 p1 ...p

N -m

(v1 ,..., vm ) = det(t Mm Mm )

(the orthogonality of {v1 ,..., vm } is not required).


Enlarge the set {v1 ,..., vm } by vectors vs for s > m to a complete orthonormal basis in R V = v1 ,..., vN is an orthogonal matrix. Let E be an N â m matrix, such that Eps = 0 except for Ess = 1 s, 1 s m. det(t Mm Mm ) = det(t E t VVE ) = det(t EE ) = det Im = 1 QED. Am (i1 ,..., im ) = 1 m!
m s=2

N



Mm = VE ,

(Im is the identity matrix of size m).

sin2 s .


Solution of the MAE for a weakly fluctuating r.h.s.
If f - f is small (relative the mean f ), then ( ) f ( ).
K +1 N
N m=2 R

( ) = f ( ) - â


...


RN

Am i 1 ,..., i
m- 1 n=1

m- 1

, i

- m- 1 n=1



n

m -1 n=1

Theorem 1. Suppose C0 > 0 , and for K = 0

K ( n) K -



n



d 1... d

m- 1

.

C1 > 0,

N m=2

(C0 + C1)m C1 , m!

RN

|f ( )| d C0, ()

RN

Then () holds true K > 0.

| K ( ) - f ( )| d C 1 .


2. Under the same conditions,
RN

|

K +1

( ) - K ( )| d C

2 RN

| K ( ) -

K -1

( )| d , ( )| ,

max |K +1( ) - K ( )| C2 max |K ( ) - R N R N hold true K > 0, where C2
N -1 m=1

K -1

(C0 + C1)m . m!

If C2 < 1, the iterations converge to a solution, unique in the bal l
RN

For N = 3, C0 =



|( ) - f ( )| d C1.

Note that C0 < 1 and hence RN |f ( )| d C0 implies f > 0 everywhere. Note that the domain of the solution is non-compact (the entire RN ).

3 - 4/3, C1 = 1/3.


The "convolution" form of the MAE
Set arg( ( )) = arg(u ( ))/2, if | arg(u ( ))| < ; (x ) ( ) u( )/(2 )
N

( (0) = 0),

R

N

arg( ( )) = - arg( (- )), if arg(u ( )) = . ( )ei ·x d is a real-valued function (not uniquely defined).
R
N

Consider the term of order m in u (note (-1)
m R
N

ei ·x d = (2 )N (x)): u ( n ) expi
m
1

...
R

R

N

A

m

1 ,..., A
m

= -(2 ) = (2 )

Nm
N

...

R

N

( 1 ) 1 ,..., ( m )
R
N

m

m

n=1





m n=1

n · x d 1 ... d
m



m

-Nm

R

N

...

R

N

A

m

( x ) e - i
m n=1

·x

1

dx1 ,...,

expi
R m
N

n=1

( x ) e - i

n · x d 1 ... d
m

m

·x

m

dx

m

=

1 m!

â expi
R
N



n · x d 1 ... d



...

R

N

det

t

(x1 ),..., (xm )

(x - x1 ),..., (x - xm )

dx1 ... dxm .


The "convolution" form of the MAE: 1+ â
N m=1

1 m!

RN

...

RN

det

t

(x1),..., (xm) dx1... dxm = f/cN .

(x - x1),..., (x - xm)

A test problem with a cosmological flavour
The reconstruction of the dynamical history of the Universe from present observations of the spatial distribution of masses: · The distribution of masses at the epoch of matter-radiation decoupling (380 000 years after the Big Bang) is uniform. · For the dynamics of matter on sufficiently large scales (of the order of a few million light years) the Lagrangian map from initial to current mass locations is approximately the gradient of a convex potential. This holds exactly for the Zel'dovich approximation and for its refinement, the first-order Lagrangian perturbation approximation. · The potential of the map satisfies the MAE, whose r.h.s. is the (known) ratio of densities at the present and initial positions. · The r.h.s. being positive, the solution is a convex function.


We solve the MAE for N = 3 and a mass distribution for G "localised ob jects" ( is small): f = The g -th "ob ject" of mass m f
(g ) -3 G

f
g =1

(g )



r-r
(g )

(g )



.

> 0islocated at r ( r) =

and has a Gaussian density distribution:

(g )

m( g ) exp(-|r/ (g) |2 ). (g ) ) 3 (

c We seek a solution u = 2 |x|2 + u with a space-periodic zero-mean u , assuming that the space-periodic r.h.s. f of the MAE is the sum of "clones" of f (g) over all periodicity cells.

The total mass is normalised:
G

f ( r) d r =
g =1

m

(g )

=1



c = 1.


An exact solution to the MAE for a spherically symmetric mass distribution
In spherical coordinates centred at r(g) , for spherically symmetric u and f , the MAE is where = |r - r(g) |. A solution, u(, ) =
0 -2

2u u 2

2

= -3 f,

3

r / 0

r f (r) dr

2

1/3

dr ,

has uniformly bounded first derivatives: u xi = 3 xi and O( -1 ) second derivatives: 2u xi x =
m - 1 xi xm 2 / 0

r f (r)dr

2

1/3

= O ( 0 )



2

f

3

/ 2 r 0

f (r)dr

-2/3

xi x m 1 +- 2 3
i m

/ 2 r 0

1/3

f (r)dr ,


since · |xi xm /2 | 1; · for < , f (/ ) is uniformly bounded and c(/ )3
/ 0

r2 f (r)dr c(/ )

3

for some positive constant c and c; · for , (/ )2 f (/ ) is uniformly bounded and 0 / 0

r2 f (r)dr c


A one-cell solution to the MAE for G > 1 spherically symmetric localised ob jects
A solution for G > 1 ob jects can be expected to admit a power series expansion u(r) =
n=0

u n ( r, ) n ,

()

where n (by analogy with the one-ob ject solution) un (r, ) = O( 0 ) and | un (r, )| = O( 0 ). We can expect interaction of one-ob ject solutions u(g) (|r - r(g) |, ) to be asymptotically unimportant: u0 (r) = u (g ) ( | r - r(g ) | , ) .
g =0

Consider a neighbourhood of an ob ject . In the leading order, the l.h.s. of the MAE is ) det u0 xi xj . For g = , |u(g) j | = O( 0 ) only triple products of u(xj contribute in the xi x xi leading order, O( -3 ); they do match the l.h.s. (Nonlinear interaction of pairs of one-ob ject solutions is O( -2 ), and that of triplets O( -1 ).) Nevertheless, () is NOT the leading order term: it is O(|r|) for |r| , and not O(|r|2 ). Furthermore, we cannot construct a global periodic solution summing up the periodically distributed "clones" of one-cell solutions () ­ the sum is infinite; there would be infinitely many pairwise interactions between ob jects yielding products of the order O( -2 ), and the pairwise interaction does not become weaker when the distance between the interacting ob jects grows.


Numerical solution of the o dd-dimensional MAE A basic solver, FPAR
An observation: If the kernels Am in the Fourier integral form
N

( )+ â
m- 1 n=1

are frozen: Am = am , the new system of equations is just the Fourier form of the equation
N m=1

( n ) -

m=2 R

N

...

R

N

A

m

i 1 ,..., i d 1 ... d

m-1

, i

-

m- 1 n=1



n

m-1 n=1



n

m- 1

= f ( ),

= 0

am m = f/ f

in the physical space.


The Fixed Point Algorithm for the Regular part of the MAE: For =
2 N

u , transform the MAE:
N m=1 m=1

am m = f/ f + F ( ),
j

where F ( )

am m - det
N m=1

-2

(

xi x

)+

ij

.
K -1

Iterate in the physical space:

m am K = f/ f + F (

).

0 Am 1/m! for m > 0 it is practical to set a1 = 1; am = 1/(2m!) for m > 1. For a1 = 1, the linear in u term is treated exactly,
N

F ( ) =

m=2 R

N

...

R

N

(Am (i 1 ,..., i m ) - am )



m n=1

( n ) expi





m n=1

n · x d 1 ... d m - 1,



and for m 2, for our am , |Am (i 1 ,..., i m ) - am | am . The l.h.s. "captures" the nonlinear behaviour of the l.h.s. of the MAE, and the algorithm has chances to converge. Since |Am - am | = am for the extreme values of Am , convergence is not guaranteed.


Lemma. odd N and the chosen am , PN ( ) has a unique root for any r.h.s. Pro of .
N m=1

am m =

At a minimum of PN ( ), PN ( ) = 0 PN > 0 identically. QED. (Similarly, even N the number of roots of PN is 0 or 2.) In an application to a test problem inspired by cosmology produces a sequence of iterations, initially converging, but subsequently blowing up.

1 N -1 PN ( ) - PN ( ) = 1+ . 2 (N - 1)!






A more advanced solver, ACPDM (Algorithm with Continuation in a Parameter and Discrepancy Minimisation) I. Continuation in the parameter p [0, 1]: A generalisation of the MAE: Q( ; p)
:::::::::::::::::::::::::::::::::::::::::::::::::::::::

N m=1

For p = 0, this is a set of polynomial equations of degree N ; for p = 1 this is the MAE. Solve () for some pj increasing from 0 to 1, with application of the FPAR iterations
N m=1 m am K = f/ f + pF ( K -1

am m - f/ f - pF ( ) = 0,

()

).

Obtain an initial approximation of the solution for p = pj by polynomial extrapolation over all nodes pj < pj (use quadruple precision). We monitor the r.m.s. discrepancy d( , p) (Note Q( ;1) = 0 d( , 1) det
-2

Q( ; p ) - Q( ; p )
xi x
j

2

.

+

ij

- f/ f

2

.)


II. Minimisation of the residual in Krylov spaces:
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

(The idea is borrowed from, e.g., GMRES; actually we employ a pJFNK without N.) If and are approximate solutions to the gMAE, Q( ; p) = 0, then Q( ; p) - Q( ; p) = J ( - )+ O( -
2

),

where J is the operator of linearisation of the gMAE around the solution. Let (·, ·) denote a scalar product, e.g., of the Lebesgue space L2 (T 3 ): (u, v) = and · the induced norm of a scalar field: v (v, v).
T
3

u(x )v (x )d x

ACPDM proceeds by sequences of stabilised iterations. (f > 0 is not required.) To make an iteration, an approximate solution K , and two sets of S scalar fields, vs (x) and ws (x), 0 s S , where 0 S Smax (in our runs, Smax = 5), are required. All vs are mutually (·, ·)-orthogonal, and vs = Jws + O( ws 2 ) for small ws . A sequence of stabilised iterations of ACPDM is initialised using the current approximation 0 for computation of an FPAR iteration 1 = 1 , and setting S = 0.


For K > 1, a ::::::::::::::::::::::::::::::::::::::::: consists of 6 steps: stabilised iteration of ACPDM i. Make an FPAR iteration, i.e., at each grid point in the physical space solve
N m=1

am (K )m = f/ f + pF (K ).

ii. Compute the discrepancy field Q(K ). iii. Orthogonalise v Q(K ) - Q(K -1 ) to all vs , 1 s S , and set v
S +1

=v -

S s=1

(v ,vs ) vs , (vs ,vs )

w

S +1

= K -

K -1

-

S s=1

(v ,vs ) ws . (vs ,vs )

iv . Minimise the discrepancy Q(K -
K +1

S s=1 qs

ws ) , setting (Q(K ),vs ) ws . (vs ,vs )

= K -

S +1 s=1

v . If S < Smax , increase S by 1; otherwise (i.e., if S = Smax ), discard v1 and w1 , and decrease by 1 the indices s of the remaining Smax fields vs and ws . vi. Compute the r.h.s. of the equation defining the iterations for K +1 , Q(K +1 ), and d(K +1 ). If d(K +1 ) < , then K +1 is the approximate solution for the present p. If Q(K +1 ) > Q(K ) , the current sequence is terminated (the effect of the nonlinearity is significant; need small Smax ).


Summary. · ACPDM performs continuation in the parameter p. · For a given p, start by carrying out the FPAR iterations. · When the convergence slows down (in our runs d2(K ) > d2( start a sequence of stabilised iterations.
K -1

)/2),

· If the sequence terminates in step vi, make the FPAR iterations while d2(K ) > d2(K -1); upon termination, start a new sequence.
modes to set in; we employed = 1.05).

( can slightly exceed 1 in order to allow the transients to die off and the dominant instability


Application to a test problem inspired by cosmology
A poor man's Universe: G = 3 galaxies in the periodicity cell T 3 , -1/2 xi 1/2: m
(1)

: m(2) : m r
(1)

(3)

=

111 : :, 18 27 32 1 =, 6
(3)

= 1, r
(1)

=

(2)

1 = (-1, 1, -1), 4

(2)

1 = (1, -1, -1), 4

1 =, 8 1 r(3) = (-1, -1, 1). 4

· The r.h.s. achieves its maximum at r

(3)

· The contrast number is 1.313 â 107 . · A uniform grid is comprised of 643 points in the periodicity cell T 3 . · det uxi xj + ij is computed by the pseudospectral method with dealiasing (for N = 3 this requires to compute uxi xj in the physical space on the grid of 1283 points). · The fall off of the energy spectrum of = 2 u is fast (11 orders of magnitude).

1 and the minimum at (1, 1, 1); 4


When applied for J = 20 nodes pj = j/J , j = 0,..., J - 1, convergence of ACPDM is fast. A polynomial 20-node extrapolation yields an approximate solution with the accuracy d( , 1) = 0.47 â 10-2 and d ( , 1) = 0.02, where d ( , 1) max det u 3
T xi x
j

+

ij

- f/ f .

When it is used as an initial approximation for a run for p = 1, the pattern of convergence is erratic; after 38 685 evaluations of the determinant of the Hessian, d( , 1) = 0.99 â 10-9 and d ( , 1) = 0.39 â 10-7 , and convergence stalls. ACPDM ceases to stall, after J = 13 nodes are added to the mesh: pj = j/J, j = 0,..., J - 1; p
J +j -1

= 1 - (2j J )-1 , j = 1,..., J . solution u h er 1 8 8 5 e accuracies on a 3. 16 for p = 1 valuations 10-10 and GHz Intel

A polynomial extrapolation over the 33 nodes yields an approximate to the accuracy d( ) = 0.78 â 10-7 and d ( ) = 0.20 â 10-5 . After furt of the determinant (9 814 in total) ACPDM yields a solution with the 0.26 â 10-8 , respectively. The duration of the run of a sequential code Core Duo processor is 1 hr. 46 min. 24 sec. (wallclock time).


Number of evaluations of the determinant of the Hessian (vertical axis, logarithmic scale) performed by ACPDM in successive computations of solutions (pj ) to the gMAE to the accuracy d( , pj ) < 10-10 . Horizontal axis: the index j numbering consecutive nodes pj .


(-0.5,-0.5,-0.5) Isosurfaces of the solution to the test MAE at the levels of a half and 1/8 of the maximum. The periodicity cell T 3 of u is shown.


(-0.5,-0.5,-0.5) Isosurfaces of
2

u for the solution to the test MAE at the level of 1/3 of the maximum.


x

3

x

3

x

2

x

1

x

2

x

1

u for the solution to the test MAE on cross sections of the periodicity cell, parallel to coordinate planes and containing pairs of ob jects (left to right): x2 = -1/4, x1 = -1/4, x3 = -1/4. (Due to the hidden symmetry of u about each of the three planes, components of gradients normal to the planes are zero.) The labels xi refer to the Cartesian coordinate axes, parallel to sides of the cross sections. Stars show locations of the three localised ob jects on the cross sections. Gray-scaling codes the masses of the ob jects (black, gray and white stars: the ob jects at r(g) , g = 1, 2, 3, respectively).


x

3

x

3

x

2

x

1

x

2

x

1

Isolines step 0.02 of n lines: negative values, (left to right) x2 = 0, parallel to sides of the

ormal components of u for the solution to the test MAE (dashed solid lines: zero and positive values) on Cartesian coordinate planes x1 = 0, x3 = 0. The labels xi refer to the Cartesian coordinate axes, shown cross sections of T 3 .


A solver for the MAE with an everywhere positive r.h.s., AICDM The maximum discrepancy between the 20-mode solution and the accurate one, is located at the minimum of the r.h.s., (1/4, 1/4, 1/4). det uxi ,xj < 0 at this point. The (Algorithm with Improvement of Convexity and Discrepancy Minimisation): a single node (p = 1) extension of the ACPDM, incorporating improv::me:::::::::::::v::::t:: of e nt of con :exi:y ::::::::: ::: an approximate solution, which is performed whenever d2 (K ) < d2 (K ). K and K are the numbers of the current iteration and of the iteration of the previous improvement of convexity, respectively; 1, we employed = 0.01. At each grid point in the physical space, compute the eigenvalues of (all real, since the Hessian is a symmetric matrix). Suppose -1 - min , At a po (using the
-2

( K )

xi x

j

If all > -1, the approximate solution 1 |x|2 + u is locally convex no action taken. 2

min < -1. The minimum eigenvalue would become -1, if uxi xi is increased i by and hence the Laplacian K = 2 u is increased by 3(-1 - min ). int, where min < -1 increase K by 6(-1 - min ) factor 6 instead of 3, we "overimprove" K ).

Repeat the procedure till all > -1 - d(K )/2. (The discrepancy d(K ) for the "improved" approximation K is allowed to exceed the initial discrepancy d(K ).)


We have inspected convergence of AICDM employing scalar products with weights: (u, v) =
T
3

u(x)v (x)w(x)dx,

where w(x; q ) = max 1, (f/ f )q .

Duration of the shortest run, for q = -1/2, with a sequential code on the 3.16 GHz Intel Core Duo processor is 25 min. 19 sec. (a gain 4 times compared to the ACPDM). Number of evaluations (NE) of the determinant of the Hessian performed by AICDM with the use of various weights. All runs are terminated as soon as d( , 1) < 10-10 . d( , 1) = 0.66 â 10-10 for q = -1, d( , 1) = 0.95 â 10-10 for q = -3/4 and d( , 1) = 1.0 â 10-10 in all remaining cases. q -1 NE 3169 d ( ) 0.09 â 10 q 1 /2 NE 2568 d ( ) 0.37 â 10 - 3/ 4 2619 0.14 â 10 - 1/ 2 2465 0.21 â 10 - 1/ 4 3743 0.62 â 10 0 2653 0.41 â 10 1/4 2571 0.39 â 10

-8

-8

-8

-8

-8

-8

-8

3 /4 2643 0.39 â 10

-8

1 2523 0.36 â 10

-8

5/4 2489 0.38 â 10

-8

3/2 2481 0.39 â 10

-8

7/4 2505 0.41 â 10

-8

2 2526 0.37 â 10

-8


An enhanced version of AICDM, EAICDM
::::::::::::::::::::::::::

Gradual refinement of solutions with the resolution of (16M )3 Fourier harmonics (1 M 4 is integer). When seeking a 643 -harmonics solution of the accuracy d( , 1) < 10- , the (16M )3 computations are terminated when d( , 1) < 10- min(, 2.5M ) . A further acceleration is likely to be achieved by avoiding dealiasing. Number of evaluations (NE) of the determinant of the Hessian performed by AICDM in the course of computation of solutions of varying accuracy to the test MAE, and wallclock duration of runs (DR, seconds) by EAICDM. d ( ) 10-3 d ( ) 0.64 â 10 NE 64 DR 3.38 d ( ) 10-7 d ( ) 0.28 â 10 NE 1018 DR 81
-2

10-4 0.86 â 10 126 4.69 10-8 0.39 â 10 1515 203

-3

10-5 1.02 â 10 251 11.17 10-9 0.36 â 10 2030 435

-3

10-6 0.23 â 10 439 42.7 10-10 0.41 â 10 2653 774

-4

-5

-6

-7

-8


Conclusion
· Three novel forms of the Monge­Amp`re equation in RN : the second-order divergence, e Fourier integral and convolution forms. · Three novel algorithms for numerical solution of the MAE with the positive r.h.s. for an odd N (the positiveness not required for ACPDM). Open questions: ::::::::::::::::::::: Does the gMAE always have a solution, as long as the MAE does? Does the solution of the gMAE depend analytically on the parameter p, and hence polynomial extrapolation for p = 1 is mathematically sensible? Apparently the algorithms are still applicable for other regions (although the inversion of the Laplacian becomes more difficult). Do we still need the same am 's? What scalar product is optimal for acceleration of convergence of ACPDM? What is the asymptotics of solutions in the small parameter determining the width of the localised ob jects? How does the contrast number measure numerical complexity of the MAE and, in particular, the condition number of the linearisation near the solution (for a given spatial discretisation)? Interesting applications for our algorithms.