Vitali Vougalter
email: vitali.vougalter@uct.ac.za

Vitaly Volpert
email: volpert@math.univ-lyon1.fr

On the existence of stationary solutions for some non-Fredholm integro-differential equations

Keywords: solvability conditions, non Fredholm operators, integro-differential equations, stationary solutions AMS subject classification: 35J10, 35P30, 35K57

We show the existence of stationary solutions for some reaction-diffusion type equations in the appropriate H 2 spaces using the fixed point technique when the elliptic problem contains second order differential operators with and without Fredholm property.

Introduction

Let us recall that a linear operator L acting from a Banach space E into another Banach space F satisfies the Fredholm property if its image is closed, the dimension of its kernel and the codimension of its image are finite. As a consequence, the equation Lu = f is solvable if and only if φ i (f ) = 0 for a finite number of functionals φ i from the dual space F * . These properties of Fredholm operators are widely used in many methods of linear and nonlinear analysis.

Elliptic problems in bounded domains with a sufficiently smooth boundary satisfy the Fredholm property if the ellipticity condition, proper ellipticity and Lopatinskii conditions are satisfied (see e.g. [START_REF] Agranovich | Elliptic boundary problems[END_REF], [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF], [START_REF] Volevich | Solubility of boundary problems for general elliptic systems[END_REF]). This is the main result of the theory of linear elliptic problems. In the case of unbounded domains, these conditions may not be sufficient and the Fredholm property may not be satisfied. For example, Laplace operator, Lu = ∆u, in R d does not satisfy the Fredholm property when considered in Hölder spaces, L :

C 2+α (R d ) → C α (R d ), or in Sobolev spaces, L : H 2 (R d ) → L 2 (R d ).
Linear elliptic problems in unbounded domains satisfy the Fredholm property if and only if, in addition to the conditions cited above, limiting operators are invertible (see [START_REF] Volpert | Elliptic partial differential equations[END_REF]). In some simple cases, limiting operators can be explicitly constructed. For example, if

1 Lu = a(x)u ′′ + b(x)u ′ + c(x)u, x ∈ R,
where the coefficients of the operator have limits at infinity,

a ± = lim
L ± u = a ± u ′′ + b ± u ′ + c ± u.
Since the coefficients are constant, the essential spectrum of the operator, that is the set of complex numbers λ for which the operator Lλ does not satisfy the Fredholm property, can be explicitly found by means of the Fourier transform:

λ ± (ξ) = -a ± ξ 2 + b ± iξ + c ± , ξ ∈ R.
Invertibility of limiting operators is equivalent to the condition that the essential spectrum does not contain the origin.

In the case of general elliptic problems, the same assertions hold true. The Fredholm property is satisfied if the essential spectrum does not contain the origin or if the limiting operators are invertible. However, these conditions may not be explicitly written.

In the case of non-Fredholm operators the usual solvability conditions may not be applicable and solvability conditions are, in general, not known. There are some classes of operators for which solvability conditions are obtained. Let us illustrate them with the following example. Consider the equation

Lu ≡ ∆u + au = f (1.1)
in R d , where a is a positive constant. The operator L coincides with its limiting operators.

The homogeneous equation has a nonzero bounded solution. Hence the Fredholm property is not satisfied. However, since the operator has constant coefficients, we can apply the Fourier transform and find the solution explicitly. Solvability conditions can be formulated as follows. If f ∈ L 2 (R d ) and xf ∈ L 1 (R d ), then there exist a solution of this equation in

H 2 (R d ) if and only if f (x), e ipx (2π) d 2 L 2 (R d ) = 0, p ∈ S d √ a
a.e.

(see [START_REF] Volpert | Solvability conditions for some linear and nonlinear non-Fredholm elliptic problems[END_REF]). Here and further down S d r denotes the sphere in R d of radius r centered at the origin. Thus, though the operator does not satisfy the Fredholm property, solvability conditions are formulated in a similar way. However, this similarity is only formal since the range of the operator is not closed.

In the case of the operator with a potential,

Lu ≡ ∆u + a(x)u = f,
Fourier transform is not directly applicable. Nevertheless, solvability conditions in R 3 can be obtained by a rather sophisticated application of the theory of self-adjoint operators (see [START_REF] Vougalter | Solvability conditions for some non Fredholm operators[END_REF]). As before, solvability conditions are formulated in terms of orthogonality to solutions of the homogeneous adjoint equation. There are several other examples of linear elliptic operators without Fredholm property for which solvability conditions can be obtained (see [START_REF] Volpert | Elliptic partial differential equations[END_REF]- [START_REF] Volpert | Solvability conditions for some linear and nonlinear non-Fredholm elliptic problems[END_REF]). Solvability conditions play an important role in the analysis of nonlinear elliptic problems. In the case of non-Fredholm operators, in spite of some progress in understanding of linear problems, there exist only few examples where nonlinear non-Fredholm operators are analyzed (see [START_REF] Ducrot | Systemes de réaction-diffusion sans propriété de Fredholm[END_REF]- [START_REF] Ducrot | Reaction-diffusion waves (with the Lewis number different from 1)[END_REF]). In the present article we consider another class of nonlinear equations, for which the Fredholm property may not be satisfied:

∂u ∂t = ∆u + au + Ω G(x -y)F (u(y), y)dy = 0, a ≥ 0. (1.2)
Here Ω is a domain in R d , d = 1, 2, 3, the more physically interesting dimensions. In population dynamics the integro-differential equations describe models with intra-specific competition and nonlocal consumption of resources (see e.g. [START_REF] Apreutesei | Vougalter: Spatial Structures and Generalized Travelling Waves for an Integro-Differential Equation[END_REF], [START_REF] Berestycki | The non-local Fisher-KPP equation: traveling waves and steady states[END_REF], [START_REF] Genieys | Pattern and waves for a model in population dynamics with nonlocal consumption of resources[END_REF]). The linear part of the corresponding operator is the same as in equation (1.1) above. We will use the explicit form of solvability conditions and will study the existence of stationary solutions of the nonlinear equation.

Formulation of the results

The nonlinear part of equation (1.2) will satisfy the following regularity conditions.

Assumption 1. Function F (u, x) : R × Ω → R is such that |F (u, x)| ≤ k|u| + h(x) f or u ∈ R, x ∈ Ω (2.1)
with a constant k > 0 and h(x

) : Ω → R + , h(x) ∈ L 2 (Ω). Moreover, it is a Lipschitz continuous function, such that |F (u 1 , x) -F (u 2 , x)| ≤ l|u 1 -u 2 | f or any u 1,2 ∈ R, x ∈ Ω (2.2)
with a constant l > 0.

Clearly, the stationary solutions of (1.2), if they exist, will satisfy the nonlocal elliptic equation ∆u + Ω G(xy)F (u(y), y)dy + au = 0, a ≥ 0.

Let us introduce the auxiliary problem

-∆u -au = Ω G(x -y)F (v(y), y)dy.
(2.3)

We denote (f 1 (x), f 2 (x)) L 2 (Ω) := Ω f 1 (x) f2 (x)dx, with a slight abuse of notations when these functions are not square integrable, like for instance those used in the one dimensional Lemma A1 of the Appendix. In the first part of the article we study the case of Ω = R d , such that the appropriate Sobolev space is equipped with the norm

u 2 H 2 (R d ) := u 2 L 2 (R d ) + ∆u 2 L 2 (R d ) .
The main issue for the problem above is that the operator -∆-a :

H 2 (R d ) → L 2 (R d
), a ≥ 0 does not satisfy the Fredholm property, which is the obstacle to solve equation (2.3). The similar situations but in linear problems, both self-adjoint and non self-adjoint involving non Fredholm second or fourth order differential operators or even systems of equations with non Fredholm operators have been studied extensively in recent years (see [START_REF] Vougalter | Solvability conditions for some non Fredholm operators[END_REF]- [START_REF] Vougalter | Solvability conditions for some systems with non Fredholm operators[END_REF]). However, we manage to show that equation (2.3) in this case defines a map T a :

H 2 (R d ) → H 2 (R d ), a ≥ 0,
which is a strict contraction under certain technical conditions.

Theorem 1. Let Ω = R d , G(x) : R d → R, G(x) ∈ L 1 (R d )
and Assumption 1 holds. I) When a > 0 we assume that xG(x) ∈ L 1 (R d ), orthogonality relations (6.4) hold if d = 1 and (6.9) when d = 2, 3 and √ 2(2π)

d 2 N a, d l < 1. Then the map T a v = u on H 2 (R d ) defined by equation (2.
3) has a unique fixed point v a , which is the only stationary solution of problem (1.2) in H 2 (R d ).

II) When a = 0 we assume that x 2 G(x) ∈ L 1 (R d ), orthogonality relations (6.10) hold, d = 1, 2, 3 and √ 2(2π)

d 2 N 0, d l < 1.
Then the map T 0 v = u on H 2 (R d ) defined by equation (2.3) admits a unique fixed point v 0 , which is the only stationary solution of problem (1.2) with a = 0 in H 2 (R d ).

In both cases I) and II) the fixed point v a , a ≥ 0 is nontrivial provided the intersection of supports of the Fourier transforms of functions supp F (0, x) ∩ supp G is a set of nonzero Lebesgue measure in R d .

In the second part of the work we study the analogous problem on the finite interval with periodic boundary conditions, i.e. Ω = I := [0, 2π] and the appropriate functional space is

H 2 (I) = {u(x) : I → R | u(x), u ′′ (x) ∈ L 2 (I), u(0) = u(2π), u ′ (0) = u ′ (2π)}.
Let us introduce the following auxiliary constrained subspaces

H 2 0 (I) := {u ∈ H 2 (I) | u(x), e ±in 0 x √ 2π L 2 (I) = 0}, n 0 ∈ N (2.4) and H 2 0, 0 (I) = {u ∈ H 2 (I) | (u(x), 1) L 2 (I) = 0}, (2.5)
which are Hilbert spaces as well (see e.g. Chapter 2.1 of [START_REF] Hislop | Introduction to spectral theory with applications to Schrödinger operators[END_REF]). We prove that equation (2.3) in this situation defines a map τ a , a ≥ 0 on the above mentioned spaces which will be a strict contraction under our assumptions.

Theorem 2. Let Ω = I, G(x) : I → R, G(x) ∈ L 1 (I), G(0) = G(2π), F (u, 0) = F (u, 2π
) for u ∈ R and Assumption 1 holds. I) When a > 0 and a = n 2 , n ∈ Z we assume that 2 √ πN a l < 1. Then the map τ a v = u on H 2 (I) defined by equation ( 2.3) has a unique fixed point v a , the only stationary solution of problem (1.2) in H 2 (I).

II) When a = n 2 0 , n 0 ∈ N assume that orthogonality relations (6.17) hold and 2

√ πN n 2 0 l < 1. Then the map τ n 2 0 v = u on H 2 0 (I) defined by equation (2.
3) has a unique fixed point v n 2 0 , the only stationary solution of problem (1.2) in H 2 0 (I). III) When a = 0 assume that orthogonality relation (6.18) holds and 2 √ πN 0 l < 1. Then the map τ 0 v = u on H 2 0, 0 (I) defined by equation ( 2.3) has a unique fixed point v 0 , the only stationary solution of problem (1.2) in H 2 0, 0 (I). In all cases I), II) and III) the fixed point v a , a ≥ 0 is nontrivial provided the Fourier coefficients G n F (0, x) n = 0 for some n ∈ Z.

Remark. We use the constrained subspaces H 2 0 (I) and H 2 0, 0 (I) in cases II) and III) respectively, such that the operators -

d 2 dx 2 -n 2 0 : H 2 0 (I) → L 2 (I) and - d 2 dx 2 : H 2 0, 0 (I) → L 2 (I)
, which possess the Fredholm property, have empty kernels.

We conclude the article with the studies of our problem on the product of spaces, where one is the finite interval with periodic boundary conditions as before and another is the whole space of dimension not exceeding two, such that in our notations

Ω = I × R d = [0, 2π] × R d , d = 1, 2 and x = (x 1 , x ⊥ ) with x 1 ∈ I and x ⊥ ∈ R d . The appropriate Sobolev space for the problem is H 2 (Ω) defined as {u(x) : Ω → R | u(x), ∆u(x) ∈ L 2 (Ω), u(0, x ⊥ ) = u(2π, x ⊥ ), u x 1 (0, x ⊥ ) = u x 1 (2π, x ⊥ )},
where x ⊥ ∈ R d a.e. and u x 1 stands for the derivative of u(x) with respect to the first variable x 1 . As in the whole space case covered in Theorem 1, the operator -∆a : H 2 (Ω) → L 2 (Ω), a ≥ 0 does not possess the Fredholm property. Let us show that problem (2.3) in this context defines a map t a : H 2 (Ω) → H 2 (Ω), a ≥ 0, a strict contraction under appropriate technical conditions.

Theorem 3. Let Ω = I × R d , d = 1, 2, G(x) : Ω → R, G(x) ∈ L 1 (Ω), G(0, x ⊥ ) = G(2π, x ⊥ ), F (u, 0, x ⊥ ) = F (u, 2π, x ⊥ ) for x ⊥ ∈ R d a.e. and u ∈ R and Assumption 1 holds. I) When n 2 0 < a < (n 0 + 1) 2 , n 0 ∈ Z + = N ∪ {0} let x ⊥ G(x) ∈ L 1 (Ω), condition (6.29) holds if dimension d = 1 and (6.30) if d = 2 and √ 2(2π) d+1 2 M a l < 1. Then the map t a v = u on H 2 (Ω) defined by equation (2.3) has a unique fixed point v a , the only stationary solution of problem (1.2) in H 2 (Ω). II) When a = n 2 0 , n 0 ∈ N let x 2 ⊥ G(x) ∈ L 1 (Ω)
, conditions (6.25), (6.27) hold if dimension d = 1 and conditions (6.26), (6.27)

hold if d = 2 and √ 2(2π) d+1 2 M n 2 0 l < 1. Then the map t n 2 0 v = u on H 2 (Ω) defined by equation (2.
3) has a unique fixed point v n 2 0 , the only stationary solution of problem (1.2) in H 2 (Ω).

III) When a = 0 let x 2 ⊥ G(x) ∈ L 1 (Ω), conditions (6.23) hold and √ 2(2π)

d+1 2 M 0 l < 1.
Then the map t 0 v = u on H 2 (Ω) defined by equation ( 2.3) has a unique fixed point v 0 , the only stationary solution of problem (1.2) in H 2 (Ω).

In all cases I), II) and III) the fixed point v a , a ≥ 0 is nontrivial provided that for some n ∈ Z the intersection of supports of the Fourier images of functions supp F (0,

x) n ∩ supp G n is a set of nonzero Lebesgue measure in R d .
Remark. Note that the maps discussed above act on real valued functions due to the assumptions on F (u, x) and G(x) involved in the nonlocal term of (2.3).

The Whole Space Case

Proof of Theorem 1. We present the proof of the theorem in case I) and when a = 0 the argument will be similar. Let us first suppose that in the case of

Ω = R d for some v ∈ H 2 (R d ) there exist two solutions u 1,2 ∈ H 2 (R d ) of problem (2.3). Then their difference w := u 1 -u 2 ∈ H 2 (R d
) will satisfy the homogeneous problem -∆w = aw. Since the Laplacian operator acting in the whole space does not have any nontrivial square integrable eigenfunctions, w(x) vanishes a.e. in R d . Let v(x) ∈ H 2 (R d ) be arbitrary. We apply the standard Fourier transform to both sides of (2.3) and arrive at

u(p) = (2π) d 2 G(p) f (p) p 2 -a (3.1)
with f (p) denoting the Fourier image of F (v(x), x). Clearly, we have the upper bounds

| u(p)| ≤ (2π) d 2 N a, d | f (p)| and |p 2 u(p)| ≤ (2π) d 2 N a, d | f (p)|
with N a, d < ∞ by means of Lemma A1 of the Appendix in one dimension and via Lemma A2 for d = 2, 3 under orthogonality relations (6.4) and (6.9) respectively. This enables us to estimate the norm

u 2 H 2 (R d ) = u(p) 2 L 2 (R d ) + p 2 u(p) 2 L 2 (R d ) ≤ 2(2π) d N 2 a, d F (v(x), x) 2 L 2 (R d ) ,
which is finite by means of (2.1) of Assumption 1. Therefore, for any

v(x) ∈ H 2 (R d ) there is a unique solution u(x) ∈ H 2 (R d ) of problem (2.
3) with its Fourier image given by (3.1) and the map T a :

H 2 (R d ) → H 2 (R d ) is well defined. This enables us to choose arbitrarily v 1,2 (x) ∈ H 2 (R d ) such that their images u 1,2 = T a v 1,2 ∈ H 2 (R d ) and estimate | u 1 (p) -u 2 (p)| ≤ (2π) d 2 N a, d | f 1 (p) -f 2 (p)|, |p 2 u 1 (p) -p 2 u 2 (p)| ≤ (2π) d 2 N a, d | f 1 (p) -f 2 (p)|,
where f 1,2 (p) stand for the Fourier images of F (v 1,2 (x), x). For the appropriate norms of functions this yields

u 1 -u 2 2 H 2 (R d ) ≤ 2(2π) d N 2 a, d F (v 1 (x), x) -F (v 2 (x), x) 2 L 2 (R d ) . Note that v 1,2 (x) ∈ H 2 (R d ) ⊂ L ∞ (R d ), d ≤ 3
T a v 1 -T a v 2 H 2 (R d ) ≤ √ 2(2π) d 2 N a, d l v 1 -v 2 H 2 (R d )
with the constant in the right side of this estimate less than one by the assumption of the theorem. Therefore, by means of the Fixed Point Theorem, there exists a unique function

v a ∈ H 2 (R d ) with the property T a v a = v a ,

The Problem on the Finite Interval

Proof of Theorem 2. Let us demonstrate the proof of the theorem in case I) and when a = n 2 0 , n 0 ∈ N or a = 0 the ideas will be similar, using the constrained subspaces (2.4) and (2.5) respectively instead of H 2 (I). First we suppose that for v ∈ H 2 (I) there are two solutions u 1,2 ∈ H 2 (I) of problem (2.3) with Ω = I. Then function w := u 1u 2 ∈ H 2 (I) will be a solution to the problem -w ′′ = aw. But a = n 2 , n ∈ Z and therefore, it is not an eigenvalue of the operator -d 2 dx 2 on L 2 (I) with periodic boundary conditions. Therefore, w(x) vanishes a.e. in I. Suppose v(x) ∈ H 2 (I) is arbitrary. Let us apply the Fourier transform to problem (2.3) considered on the interval I which yields

u n = √ 2π G n f n n 2 -a , n ∈ Z (4.1)
with

f n := F (v(x), x) n .
Clearly for the transform of the second derivative we have

(-u ′′ ) n = √ 2π n 2 G n f n n 2 -a , n ∈ Z,
which enables us to estimate

u 2 H 2 (I) = ∞ n=-∞ |u n | 2 + ∞ n=-∞ |n 2 u n | 2 ≤ 4πN 2 a F (v(x), x) 2 L 2 (I) < ∞ due to (2.
u 1 -u 2 2 H 2 (I) = ∞ n=-∞ |u 1n -u 2n | 2 + ∞ n=-∞ |n 2 (u 1n -u 2n )| 2 ≤ ≤ 4πN 2 a F (v 1 (x), x) -F (v 2 (x), x) 2 L 2 (I) . Obviously v 1,2 (x) ∈ H 2 (I) ⊂ L ∞ (I)
due to the Sobolev embedding. By means of (2.2) we easily obtain

τ a v 1 -τ a v 2 H 2 (I) ≤ 2 √ πN a l v 1 -v 2 H 2 (I) ,
such that the constant in the right side of this upper bound is less than one as assumed. Thus, the Fixed Point Theorem implies the existence and uniqueness of a function v a ∈ H 2 (I) satisfying τ a v a = v a , which is the only stationary solution of problem (1.2) in H 2 (I). Suppose v a (x) = 0 a.e. in I. Then we obtain the contradiction to the assumption that G n F (0, x) n = 0 for some n ∈ Z. Note that in the case of a = n 2 , n ∈ Z the argument does not require any orthogonality conditions.

The Problem on the Product of Spaces

Proof of Theorem 3. We present the proof of the theorem for case II) since when the parameter a vanishes or is located on the open interval between squares of two nonnegative integers the ideas are similar. Suppose there exists v(x) ∈ H 2 (Ω) which generates u 1,2 (x) ∈ H 2 (Ω) solving equation (2.3). Then the difference w := u 1u 2 ∈ H 2 (Ω) will satisfy -∆w = n 2 0 w in our domain Ω. By applying the partial Fourier transform to this equation we easily arrive at

-∆ ⊥ w n (x ⊥ ) = (n 2 0 -n 2 )w n (x ⊥ ). Clearly w 2 L 2 (Ω) = ∞ n=-∞ w n 2 L 2 (R d ) such that w n (x ⊥ ) ∈ L 2 (R d ), n ∈ Z. Since the transversal Laplacian operator -∆ ⊥ on L 2 (R d )
does not have any nontrivial square integrable eigenfunctions, w(x) is vanishing a.e. in Ω. Let v(x) ∈ H 2 (Ω) be arbitrary. We apply the Fourier transform to both sides of problem (2.3) and obtain

u n (p) = (2π) d+1 2 G n (p) f n (p) p 2 + n 2 -n 2 0 , n ∈ Z, p ∈ R d , d = 1, 2, (5.1)
where f n (p) stands for the Fourier image of F (v(x), x). Obviously,

| u n (p)| ≤ (2π) d+1 2 M n 2 0 | f n (p)| and |(p 2 + n 2 ) u n (p)| ≤ (2π) d+1 2 M n 2 0 | f n (p)|, where M n 2
0 < ∞ by means of Lemma A5 of the Appendix under the appropriate orthogonality conditions stated in it. Thus

u 2 H 2 (Ω) = ∞ n=-∞ R d | u n (p)| 2 dp + ∞ n=-∞ R d |(p 2 + n 2 ) u n (p)| 2 dp ≤ ≤ 2(2π) d+1 M n 2 0 2 F (v(x), x) 2 L 2 (Ω)
< ∞ by means of (2.1) of Assumption 1, such that for any v(x) ∈ H 2 (Ω) there exists a unique u(x) ∈ H 2 (Ω) solving equation (2.3) with its Fourier image given by (5.1) and the map t a : H 2 (Ω) → H 2 (Ω) in case II) of the Theorem is well defined. Then we consider arbitrary v 1,2 ∈ H 2 (Ω) such that their images under the map are u 1,2 = t n 2 0 v 1,2 ∈ H 2 (Ω) and obtain

u 1 -u 2 2 H 2 (Ω) = ∞ n=-∞ R d | u 1n (p) -u 2n (p)| 2 dp + ∞ n=-∞ R d |(p 2 + n 2 )( u 1n (p) -u 2n (p))| 2 dp ≤ ≤ 2(2π) d+1 M n 2 0 2 F (v 1 (x), x) -F (v 2 (x), x) 2 L 2 (Ω) . Clearly v 1,2 ∈ H 2 (Ω) ⊂ L ∞ (Ω)
via the Sobolev embedding theorem. Using (2.2) we easily arrive at the estimate

t n 2 0 v 1 -t n 2 0 v 2 H 2 (Ω) ≤ √ 2(2π) d+1 2 M n 2 0 l v 1 -v 2 H 2 (Ω)
with the constant in the right side of it less than one by assumption. Therefore, the Fixed Point Theorem yields the existence and uniqueness of a function v n 2 0 ∈ H 2 (Ω) which satisfies t n 2 0 v n 2 0 = v n 2 0 and is the only stationary solution of problem (1.2) in H 2 (Ω) in case II) of the theorem. Suppose v n 2 0 (x) = 0 a.e. in Ω. This yields the contradiction to the assumption that there exists n ∈ Z for which supp G n ∩ supp F (0, x) n is a set of nonzero Lebesgue measure in R d .

Appendix

Let G(x) be a function, G(x) : R d → R, d ≤ 3 for which we denote its standard Fourier transform using the hat symbol as

G(p) := 1 (2π) d 2 R d G(x)e -ipx dx, p ∈ R d such that G(p) L ∞ (R d ) ≤ 1 (2π) d 2 G L 1 (R d ) (6.1)
and

G(x) = 1 (2π) d 2 R d G(q)e iqx dq, x ∈ R d .
Let us define the auxiliary quantities

N a, d := max G(p) p 2 -a L ∞ (R d ) , p 2 G(p) p 2 -a L ∞ (R d ) (6.2)
for a > 0 and

N 0, d := max G(p) p 2 L ∞ (R d ) , G(p) L ∞ (R d ) (6.3) when a = 0. Lemma A1. Let G(x) ∈ L 1 (R). a) If a > 0 and xG(x) ∈ L 1 (R) then N a, 1 < ∞ if and only if G(x), e ±i √ ax √ 2π L 2 (R) = 0. (6.4) b) If a = 0 and x 2 G(x) ∈ L 1 (R) then N 0, 1 < ∞ if and only if (G(x), 1) L 2 (R) = 0 and (G(x), x) L 2 (R) = 0. (6.5)
Proof. In order to prove part a) of the lemma we write the function

G(p) p 2 -a = G(p) p 2 -a χ I δ + G(p) p 2 -a χ I δ c , (6.6)
where χ A here and further down stands for the characteristic function of a set A, A c for its complement, the set

I δ = I + δ ∪ I - δ with I + δ = {p ∈ R | √ a -δ < p < √ a + δ}, I - δ = {p ∈ R | - √ a -δ < p < - √ a + δ} and 0 < δ < √ a.
The second term in the right side of (6.6) can be easily estimated in absolute value from above using (6.1) as

1 √ 2πδ 2 G L 1 (R) < ∞
and the remaining term in the right side of (6.6) can be written as

G(p) p 2 -a χ I + δ + G(p) p 2 -a χ I - δ .
We will use the expansions near the points of singularity given by

G(p) = G( √ a) + p √ a d G(s) ds ds, G(p) = G(- √ a) + p - √ a d G(s) ds ds with d G(p) dp L ∞ (R) ≤ 1 √ 2π xG L 1 (R) <
∞ by the assumption of the lemma. This enables us to obtain the bound

p √ a d G(s) ds ds p 2 -a χ I + δ ≤ C 2 √ a -δ < ∞, p - √ a d G(s) ds ds p 2 -a χ I - δ ≤ C 2 √ a -δ < ∞.
Therefore it remains to estimate

G( √ a) p 2 -a χ I + δ + G(- √ a) p 2 -a χ I - δ ,
which belongs to L ∞ (R) if and only if G(± √ a) = 0, which is equivalent to the orthogonality relations (6.4). To estimate the second term in the right side of (6.2) under these orthogonality relations we consider the two situations. The first one is when |p| ≤ √ a + δ and we have the bound

p 2 G(p) p 2 -a ≤ ( √ a + δ) 2 G(p) p 2 -a L ∞ (R) < ∞.
In the second one |p| > √ a + δ which yields

p 2 p 2 -a ∈ L ∞ (R) and G(p) is bounded via (6.1),
which completes the proof of part a) of the lemma. Then we turn our attention to the situation of a = 0, such that

G(p) p 2 = G(p) p 2 χ {|p|≤1} + G(p) p 2 χ {|p|>1} . (6.7)
The second term in the right side of the identity above can be easily estimated as

G(p) p 2 χ {|p|>1} ≤ G(p) L ∞ (R) < ∞ (6.8)
due to (6.1). We will make use of the representation

G(p) = G(0) + p d G dp (0) + p 0 s 0 d 2 G(q) dq 2 dq ds. Obviously d 2 G(p) dp 2 ≤ 1 √ 2π x 2 G(x) L 1 (R) <
∞ by the assumption of the lemma. Hence

p 0 s 0 d 2 G(q) dq 2 dq ds p 2 χ {|p|≤1} ≤ C 2 < ∞
and the only expression which remains to estimate is given by

G(0) p 2 + d G dp (0) p χ {|p|≤1} , which is contained in L ∞ (R) if and only if G(0) and d G dp (0)
vanish. This is equivalent to the orthogonality relations (6.5). Note that G(p) L ∞ (R) < ∞ by means of (6.1).

The proposition above can be generalized to higher dimensions in the following statement.

Lemma A2. Let G(x) ∈ L 1 (R d ), d = 2, 3. a) If a > 0 and xG(x) ∈ L 1 (R d ) then N a, d < ∞ if and only if G(x), e ipx (2π)
d 2 L 2 (R d ) = 0 f or p ∈ S d √ a
a.e. (6.9)

b) If a = 0 and x 2 G(x) ∈ L 1 (R d ) then N 0, d < ∞ if and only if (G(x), 1) L 2 (R d ) = 0 and (G(x), x k ) L 2 (R d ) = 0, 1 ≤ k ≤ d. (6.10)
Proof. To prove part a) of the lemma we introduce the auxiliary spherical layer in the space of d = 2, 3 dimensions

A δ := {p ∈ R d | √ a -δ < |p| < √ a + δ}, 0 < δ < √ a and write G(p) p 2 -a = G(p) p 2 -a χ A δ + G(p) p 2 -a χ A c δ . (6.11)
For the second term in the right side of (6.11) we have the upper bound in the absolute value as

G(p) L ∞ (R d ) √ aδ < ∞ due to (6.1). Let us expand G(p) = |p| √ a ∂ G(|s|, σ) ∂|s| d|s| + G( √ a, σ),
where σ stands for the angle variables on the sphere. Using the elementary inequality ∂ G(

d 2 xG(x) L 1 (R d ) p) ∂|p| ≤ 1 (2π)
with its right side finite by the assumption of the lemma we

estimate |p| √ a ∂ G(|s|,σ) ∂|s| d|s| p 2 -a χ Aσ ≤ C √ a < ∞.
The only remaining term G( √ a, σ)

p 2 -a χ A δ ∈ L ∞ (R d ), d = 2, 3 if and only if G( √ a, σ) vanishes
a.e. on the sphere S d √ a , which is equivalent to orthogonality relations (6.9). The proof of the fact that the second norm in the right side of (6.2) under conditions (6.9) is finite is analogous to the one presented in Lemma A1 in one dimension. For the proof of part b) of the lemma we apply the two and three dimensional analog of formula (6.7), such that for the second term in its right side there is a bound analogous to (6.8). Let us use the representation formula

G(p) = G(0) + |p| ∂ G ∂|p| (0, σ) + |p| 0 s 0 ∂ 2 G(|q|, σ) ∂|q| 2 d|q| ds. Apparently ∂ G ∂|p| (0, σ) = - i (2π) d 2 R d G(x)|x|cosθdx, (6.12)
where θ is the angle between vectors p and x in R d and for the second derivative

∂ 2 G(p) ∂|p| 2 ≤ 1 (2π) d 2 x 2 G(x) L 1 (R d ) < ∞
by the assumption of the lemma. This yields

|p| 0 s 0 ∂ 2 G(|q|,σ) ∂|q| 2 d|q| ds p 2 χ {|p|≤1} ≤ C 2 < ∞,
such that the only expression remaining to estimate is given by

G(0) p 2 + ∂ G ∂|p| (0, σ) |p| χ {|p|≤1} (6.13)
with the first derivative (6.12) containing the angular dependence. We consider first the case of d = 2 such that p = (|p|, θ p ), x = (|x|, θ x ) ∈ R 2 and the angle between them θ = θ pθ x . A straightforward computation yields that the right side of (6.12) is given by

- i 2π Q 2 1 + Q 2 2 cos(θ p -α) with Q 1 := R 2 G(x)x 1 dx, Q 2 := R 2 G(x)x 2 dx, tanα := Q 2 Q 1 (6.14)
and x = (x 1 , x 2 ) ∈ R 2 such that (6.13) is equal to

G(0) p 2 - i 2π Q 2 1 + Q 2 2 cos(θ p -α) |p| χ {|p|≤1} .
Note that the situation of Q 1 = 0 and Q 2 = 0 corresponds to the cases of α equal to π 2 or -π 2 .

Obviously, the expression above is contained in L ∞ (R 2 ) if and only if the quantities G(0), Q 1 and Q 2 vanish, which is equivalent to orthogonality relations (6.10) in two dimensions. In the case of d = 3 the argument is quite similar. The coordinates of vectors

x = (x 1 , x 2 , x 3 ) = (|x|sinθ x cosϕ x , |x|sinθ x sinϕ x , |x|cosθ x ) ∈ R 3 and p = (|p|sinθ p cosϕ p , |p|sinθ p sinϕ p , |p|cosθ p ) ∈ R 3
are being used to compute cosθ = (p, x) R 3 |p||x| involved in the right side of (6.12). Here (p, x) R 3 stands for the scalar product of the vectors in three dimensions. An easy calculation shows that when d = 3 the right side of (6.12) can be written as

- i (2π) 3 2 { Q 2 1 + Q 2 2 sinθ p cos(ϕ p -α) + Q 3 cosθ p }
with α given by (6.14) and here Q k = R 3 G(x)x k dx, k = 1, 2, 3, which are the three dimensional generalizations of the correspondent expressions given by (6.14) and term (6.13) will be equal to

G(0) p 2 - i (2π) 3 2 |p| { Q 2 1 + Q 2 2 sinθ p cos(ϕ p -α) + Q 3 cosθ p } χ {|p|≤1}
and will belong to L ∞ (R 3 ) if and only if G(0) along with Q k , k = 1, 2, 3 vanish, which is equivalent to orthogonality conditions (6.10) in three dimensions. The second norm in the right side of (6.3) is finite under relations (6.1).

Let the function G(x) : I → R, G(0) = G(2π) and its Fourier transform on the finite interval is given by

G n := 2π 0 G(x) e -inx √ 2π dx, n ∈ Z and G(x) = ∞ n=-∞ G n e inx √ 2π
. Similarly to the whole space case we define

N a := max G n n 2 -a l ∞ , n 2 G n n 2 -a l ∞ (6.15)
for a > 0. In the situation of a = 0

N 0 := max G n n 2 l ∞ , G n l ∞ . ( 6.16)
We have the following elementary statement. Proof. Clearly we have the bound

G n l ∞ ≤ 1 √ 2π G L 1 (I) < ∞. (6.19)
Thus in case a) when a = n 2 , n ∈ Z the expressions under the norms in the right side of (6.15) do not contain any singularities and the result of the lemma is obvious. When a = n 2 0 for some n 0 ∈ N or a = 0 conditions (6.17) and (6.18) respectively are necessary and sufficient for eliminating the existing singularities by making the corresponding Fourier coefficients equal to zero: G ±n 0 in case b) and G 0 in case c).

Let G(x) be a function on the product of spaces studied in Theorem 3, G(x) : Ω

= I × R d → R, d = 1, 2, G(0, x ⊥ ) = G(2π, x ⊥ ) for x ⊥ ∈ R d a.
e. and its Fourier transform on the product of spaces equals to

G n (p) := 1 (2π) d+1 2 R d dx ⊥ e -ipx ⊥ 2π 0 G(x 1 , x ⊥ )e -inx 1 dx 1 , p ∈ R d , n ∈ Z such that G n (p) L ∞ n,p := sup {p∈R d , n∈Z} | G n (p)| ≤ 1 (2π) d+1 2 G L 1 (Ω) (6.20) and G(x) = 1 (2π) d+1 2 ∞ n=-∞ R d G n (p)e ipx ⊥ e inx 1 dp.
It is also useful to consider the Fourier transform only in the first variable, such that

G n (x ⊥ ) := 2π 0 G(x 1 , x ⊥ ) e -inx 1 √ 2π dx 1 , n ∈ Z. Let us introduce ξ a n (p) := G n (p) p 2 + n 2 -a
and define

M a := max{ ξ a n (p) L ∞ n,p , (p 2 + n 2 )ξ a n (p) L ∞ n,p } (6.21)
for a > 0 and

M 0 := max G n (p) p 2 + n 2 L ∞ n,p , G n (p) L ∞ n,p (6.22)
when a = 0. Here the momentum vector p ∈ R d .

Lemma A4. Let G(x) ∈ L 1 (Ω), x 2 ⊥ G(x) ∈ L 1 (Ω) and G(0, x ⊥ ) = G(2π, x ⊥ ) for x ⊥ ∈ R d a.e., d = 1, 2. Then M 0 < ∞ if and only if (G(x), 1) L 2 (Ω) = 0, (G(x), x ⊥, k ) L 2 (Ω) = 0, 1 ≤ k ≤ d, d = 1, 2. (6.23) Proof. Let us expand G n (p) p 2 + n 2 = G 0 (p) p 2 χ {p∈R d , n=0} + G n (p) p 2 + n 2 χ {p∈R d , n∈Z, n =0} .
The second term in the right side of this identity can be estimated above in the absolute value by means of (6.20) by 1

(2π) d+1 2 G L 1 (Ω) < ∞.
Clearly we have the bound on the norm

x 2 ⊥ G n (x ⊥ ) L 1 (R d ) ≤ 1 √ 2π 2π 0 dx 1 R d dx ⊥ x 2 ⊥ |G(x)| < ∞, n ∈ Z (6.24)
by the assumption of the lemma. Thus the term G 0 (p) p 2 ∈ L ∞ (R d ) if and only if the orthogonality conditions (6.23) hold, which is guaranteed for d = 1 by Lemma A1 and when dimension d = 2 by Lemma A2. Note that the last term in the right side of (6.22) is bounded via (6.20).

Next we turn our attention to the situation when the parameter a is nontrivial.

Lemma A5. Let G(x) ∈ L 1 (Ω), x 2 ⊥ G(x) ∈ L 1 (Ω) and G(0, x ⊥ ) = G(2π, x ⊥ ) for x ⊥ ∈ R d a.e., d = 1, 2 and a = n 2 0 , n 0 ∈ N. Then M a < ∞ if and only if G(x 1 , x ⊥ ), e inx 1 √ 2π e ±i √ n 2 0 -n 2 x ⊥ √ 2π L 2 (Ω) = 0, |n| ≤ n 0 -1, d = 1, (6.25) G(x 1 , x ⊥ ), e inx 1 √ 2π e ipx ⊥ 2π L 2 (Ω) = 0, p ∈ S 2 √ n 2 0 -n 2 a.e., |n| ≤ n 0 -1, d = 2, (6.26) G(x 1 , x ⊥ ), e ±in 0 x 1 √ 2π L 2 (Ω) = 0, G(x 1 , x ⊥ ), e ±in 0 x 1 √ 2π x ⊥, k L 2 (Ω) = 0, 1 ≤ k ≤ d. (6.27)
Proof. We will use the representation of the function ξ a n (p), n ∈ Z, p ∈ R d as the sum

ξ a n (p)χ {p∈R d , |n|>n 0 } + ξ a n (p)χ {p∈R d , |n|<n 0 } + ξ a n 0 (p)χ {p∈R d , n=n 0 } + ξ a -n 0 (p)χ {p∈R d , n=-n 0 } . (6.28) Obviously |ξ a n (p)χ {p∈R d , |n|>n 0 } | ≤ G n (p) L ∞ n,p <
∞ by means of (6.20). We have trivial estimates on the norms for n ∈ Z

G n (x ⊥ ) L 1 (R d ) ≤ 1 √ 2π 2π 0 dx 1 R d dx ⊥ |G(x 1 , x ⊥ )| < ∞ and x ⊥ G n (x ⊥ ) L 1 (R d ) ≤ 1 √ 2π 2π 0 dx 1 R d dx ⊥ |x ⊥ ||G(x 1 , x ⊥ )| < ∞.
Note that G(x) ∈ L 1 (Ω) and x 2 ⊥ G(x) ∈ L 1 (Ω) by the assumptions of the lemma, which yields x ⊥ G(x) ∈ L 1 (Ω). Thus when dimension d = 1 by means of Lemma A1 ξ a n (p)χ {p∈R d , |n|<n 0 } ∈ L ∞ n,p if and only if orthogonality relations (6.25) hold. For d = 2 the necessary and sufficient conditions for the boundedness of the second term in (6.28) via Lemma A2 are given by (6.26). Lemmas A1 and A2 yield that the third term in (6.28) belongs to L ∞ n,p if and only if conditions (6.27) with the positive sign under the exponents are satisfied. Clearly x 2 ⊥ G n (x ⊥ ) ∈ L 1 (R d ) due to the assumption of the lemma and estimate (6.24). Similarly, we obtain that the necessary and sufficient conditions for the the last term in (6.28) to be contained in L ∞ n,p are given by (6.27) with the negative sign under the exponents. in which the absolute value of the first term has the upper bound (n 2 0 + 1) ξ a n (p) L ∞ n,p < ∞ under the orthogonality conditions of the lemma and of the second one (1+n 2 0 ) G n (p) L ∞ n,p < ∞ via (6.20).

Finally, we study the case when the parameter a is located on an open interval between the squares of two consecutive nonnegative integers.

Lemma A6. Let G(x) ∈ L 1 (Ω), x ⊥ G(x) ∈ L 1 (Ω) and G(0, x ⊥ ) = G(2π, x ⊥ ) for x ⊥ ∈ R d a.e., d = 1, 2 and n 2 0 < a < (n 0 + 1) in which the first and the second terms can be easily bounded above in their absolute values by the quantities finite under the conditions of the lemma, namely 1 + a (n 0 + 1) 2a G n (p) L ∞ n,p and (n 0 + 1) 2 ξ a n (p) L ∞ n,p respectively.

  x→±∞ a(x), b ± = lim x→±∞ b(x), c ± = lim x→±∞ c(x), the limiting operators are:

Lemma A3 .

 A3 Let G(x) ∈ L 1 (I) and G(0) = G(2π). a) If a > 0 and a = n 2 , n ∈ Z then N a < ∞. b) If a = n 2 0 , n 0 ∈ N then N a < ∞ if and only if G(x), e ±in 0 x √ 2π L 2 (I) = 0. (6.17) c) If a = 0 then N 0 < ∞ if and only if (G(x), 1) L 2 (I) = 0. (6.18)

  by means of the Sobolev embedding. Using condition (2.2) we easily arrive at

  which is the only stationary solution of equation (1.2) in H 2 (R d ). Suppose v a (x) vanishes a.e. in R d . This will contradict to the assumption that the Fourier images of G(x) and F (0, x) do not vanish on a set of nonzero Lebesgue measure in R d .

  [START_REF] Agranovich | Elliptic boundary problems[END_REF] of Assumption 1 and Lemma A3 of the Appendix. Hence, for an arbitrary 4.1) and the map τ a : H 2 (I) → H 2 (I) in case I) is well defined. Let us consider any v 1,2 ∈ H 2 (I) with their images under the map mentioned above u 1,2 = τ a v 1,2 ∈ H 2 (I) and arrive easily at the upper bound

v(x) ∈ H 2 (I) there is a unique u(x) ∈ H 2 (I) solving equation (2.3) with its Fourier image given by (

  Then we represent (p 2 + n 2 )ξ a n (p) as the sum (p 2 + n 2 )ξ a n (p)χ {p∈R d , n∈Z, p 2 +n 2 ≤n 2 0 +1} + (p 2 + n 2 )ξ a n (p)χ {p∈R d , n∈Z, p 2 +n 2 >n 2

	0 +1}

  2 , n 0 ∈ Z + = N ∪ {0}. Then M a < ∞ if and only ifG(x 1 , x ⊥ ),Proof. Let us expand ξ a n (p) as the sum of two termsξ a n (p)χ {p∈R d , n∈Z, |n|≥n 0 +1} + ξ a n (p)χ {p∈R d , n∈Z, |n|≤n 0 } ,such that the absolute value of the first one is bounded above byG n (p) L ∞ n,p (n 0 + 1) 2a < ∞and the second one belongs to L ∞ n,p if and only if orthogonality relations (6.29) are satisfied in one dimension by means of Lemma A1 and if and only if conditions (6.30) are fulfilled in two dimensions via Lemma A2. We write (p 2 + n 2 )ξ a n (p) as the sum (p 2 + n 2 )ξ a n (p)χ {p∈R d , n∈Z, p 2 +n 2 ≥(n 0 +1) 2 } + (p 2 + n 2 )ξ a n (p)χ {p∈R d , n∈Z, p 2 +n 2 <(n 0 +1) 2 }

			e inx 1 √ 2π	e ±i	√ √ a-n 2 x ⊥ 2π	L 2 (Ω)	= 0, |n| ≤ n 0 , d = 1,	(6.29)
	G(x 1 , x ⊥ ),	e inx 1 √ 2π	e ipx ⊥ 2π L 2 (Ω)	= 0, p ∈ S 2 √	a-n 2 a.e., |n| ≤ n 0 , d = 2.	(6.30)

Acknowledgement. The first author thanks James Colliander for stimulating discussions.


N Bessonov

V Volpert

Computational Algorithm for Some Problems with Variable Geometrical Structure

Keywords: plants, branching, variable structure, computational algorithm

The work is devoted to the computational algorithm for a problem of plant growth. The plant is represented as a system of connected intervals corresponding to branches. We compute the concentration distributions inside the branches. The originality of the problem is that the geometry of the plant is not a priori given. It evolves in time depending on the concentrations of plant hormones found as a solution of the problem. New branches appear in the process of plant growth. The algorithm is adapted to an arbitrary plant structure and an arbitrary number of branches.

Introduction

Plant modelling attract much attention both from the point of view of understanding of fundamental biological mechanisms and for practical purposes related crop optimization. Various approaches are developed (see [START_REF] Krug | International Symposium on Models for Plant Growth, Environmental Control and Farm Management in Protected Cultivation[END_REF], [START_REF] Godin | 4th International Workshop on Functional-Structural Plant Models[END_REF]). In this paper we describe a numerical algorithm developed in order to model plant growth and structure formation on the basis of realistic biological mechanism suggested in our previous works [START_REF] Bessonov | Dynamic Models of Plant Growth[END_REF], [START_REF] Bessonov | Branching Pattern in Plants[END_REF]. Plants are considered here as interconnected one-dimensional intervals. Each interval corresponds to a branch. There are branches of different generations. There exists only one branch of the first generation. It corresponds to the trunk. One its end is fixed and corresponds to the root, another one moves with some speed. If we introduce a space variable x along this interval, then x = 0 corresponds to its fixed end and x = L(t) to the moving boundary. We do not consider the root system and take it into account only by means of the flux of nutrients passing though the fixed boundary. The moving boundary corresponds to the apical meristem, a narrow layer of cells at the very end of the growing branch. These cells divide and determine the plant growth. The cells along the trunk differentiate and do not divide anymore.

We next introduce the concentration of nutrients C(x, t), which depends on the space variable x and on time t, and the concentration of growth and mitosis factor R(t), which depends only on time. Nutrient are supplied trough the fixed boundary of the interval and consumed at its moving boundary. Growth and mitosis factor is located in the apex. It is produced there and determines the proliferation rate of the cells, that is the rate of branch growth. Thus, we obtain a free boundary problem where the speed of the boundary is determined by the variable R(t) defined at this boundary and by the variable C(x, t) defined in the whole interval. Precise formulation of the problem will be given in the next section. Biological background of this study is presented in the previous works [START_REF] Bessonov | Dynamic Models of Plant Growth[END_REF], [START_REF] Bessonov | Branching Pattern in Plants[END_REF]. Branches of the second generation are also straight intervals with one endpoint inside the first branch (trunk) and another endpoint, which is free. It corresponds to the apex of this new branch and determines its growth. The number of branches of the second generation can be arbitrary, as well as their starting points. These are not fixed parameters of the problem.

In order to describe how branches of the second and of other generations appear, we briefly recall some biological facts. Appearance of new buds is determined by the concentrations of plants hormones auxin and cytokinin. They are produced by the plant itself and redistributed by ascending or descending fluxes. The most recent investigations reveal that auxin efflux carrier PIN1 protein plays the central role in this regulation. The expression of this protein is considered as a key factor for formation of plant organs [START_REF] Heisler | Patterns of Auxin Transport and Gene Expression during Primordium Development Revealed by Live Imaging of the Arabidopsis Inflorescence Meristem[END_REF]- [START_REF] Fleming | Formation of Primordia and Phyllotaxy[END_REF], [START_REF] Reinhardt | Auxin Regulates the Initiation and Radial Position of Plant Lateral Organs[END_REF]- [START_REF] Stieger | The Auxin Influx Carrier is Essential for Correct Leaf Positioning[END_REF]. On the other hand, regulation of the PIN genes is itself under the auxin control (via PLETHORA gene feedback loop [START_REF] Blilou | The PIN Auxin Efflux Facilitator Network Controls Growth and Patterning in Arabidopsis Roots[END_REF]). Therefore, we have an auxin induced auxin efflux which can be modelled with or without the intermediate PIN protein. However, auxin by itself cannot initiate cell proliferation. It can happen only in the presence of cytokinin, another plant hormone (actually a group of hormones), the main role of which is the regulation of cell proliferation. Therefore it can be an interplay of these two hormones that causes formation of a new bud.

Thus, if these concentrations take on some given values at some space point, then a new bud appears there. From the algorithmic point of view, buds represent very short branches with associated values of the concentration of mitosis and growth factor. If this value is low, then the bud remains dormant. If it becomes sufficiently high, then the bud gives a new growing branch. Thus, we introduce two new variables A(x, t) and K(x, t) which describe the distributions of auxin and cytokinin. They are defined inside all branches, as well as the concentration of nutrients C(x, t). The space variable x here is proper to each branch. A more precise description of the model is given below.

We summarize the model of plant growth as follows. Plant is represented as a number of connected and growing straight intervals. Their appearance, location and the rate of growth are determined by concentrations of nutrients, hormones and mitosis and growth factors defined either inside the branches or at their boundaries. The plant architecture, that is the number of intervals and their locations are not a priori given. From the mathematical point of view, we consider several free boundary problems whose solutions influence each other. The number of these problems and their relation to each other evolve in time. This is a new type of free boundary problems which we call one-dimensional problems with branching. The numerical algorithm should be able to describe any possible plant architecture. The main purpose of this work is to present this algorithm (Section 3). We illustrate its application to the study of branching patterns in Section 4.

Model

Without branching

We consider in this work one-dimensional model justified if the length (or height) L of the plant is essentially greater than the diameter of its trunk. Hence we consider the interval 0 ≤ x ≤ L(t) with the length depending on time. The left endpoint x = 0 corresponds to the root. Its role is to provide the flux of nutrients taken into account through the boundary condition. We do not model the root growth here in order to simplify the problem. Therefore the left boundary is fixed. The right endpoint, x = L(t) corresponds to the apex. Its width is much less than that of the plant. We suppose in the model that it is a mathematical point. The value L(t) increases over time. According to the assumption above, the growth rate is determined by the concentration of metabolites at x = L(t), which we denote by R. Thus

dL dt = f (R). (2.1)
The function f (R) will be specified below. We recall that the interval 0 < x < L(t) corresponds to differentiated cells that conduct nutrients from the root to the apex. We suppose that they are in a liquid solution. Denote by C their concentration, which is a function of x and t. Its evolution is described by the diffusion-advection equation

∂C ∂t + u ∂C ∂x = d ∂ 2 C ∂x 2 . (2.2)
Here u is the velocity of the fluid, and d is the diffusion coefficient. Assuming that the fluid is incompressible and fills the xylem uniformly (the part of the plant tissue conducting nutrients from below to above and located inside the cambium layer), we obtain u = dL dt .

We complete equation (2.2) by setting the boundary conditions

x = 0 : C = 1, x = L(t) : d ∂C ∂x = -g(R)C. (2.
3)

The second boundary condition shows that the flux of nutrients from the main body of the plant to the meristem is proportional to the concentration C(L, t). This is a conventional relation for mass exchange at the boundary, Robin boundary conditions. The factor g(R)

shows that this flux can be regulated by proliferating cells. We discuss this assumption as well as the form of the function g(R) below. We now derive the equation describing the evolution of R. At this point we need to return to the model in which the width of the meristem is finite. We denote it by h. Then we have

h dR dt = g(R)C -σR. (2.4)
The first term in the right-hand side of this equation describes production of the GM-factor R in the meristem. The second term corresponds to its consumption. System of equations (2.1)-(2.4) is a generic one-dimensional model of plant growth based on: a) "continuous medium" assumptions of mass conservation (for C + R) and of the proportionality of the flux ∂C/∂x at the boundary to the value of C; and b) a "biological" assumption that there is a chemical species R, the plant growth and mitosis factor, which is produced in the meristem and which determines the plant growth.

We note that the conservation of mass in the case σ = 0 implies that the term g(R)C enters both the boundary condition and equation (2.4). Therefore, the assumption that the rate of the plant growth factor production depends on its concentration R makes the boundary condition depend on it also. We will see below that properties of the function g can be crucial for plant growth. In particular, if it is constant (the production rate is not auto-catalytic), we will not be able to describe the essential difference in plant sizes.

We now specify the form of the functions f and g. We will consider f as a piece-wise constant function equal to 0 if R is less than a critical value R f and equal to some positive constant f 0 if R is greater than R f (Figure 1a). This means that growth begins if the concentration of the plant growth factor exceeds some critical value.

The production of the growth factor R is assumed to be auto-catalytic. To simplify the model, we consider a piece-wise linear function g(R) (Figure 1b). In some cases we also consider smooth functions f and g. These assumptions are consistent with plant morphogenesis. They are discussed in [START_REF] Bessonov | Dynamic Models of Plant Growth[END_REF], [START_REF] Bessonov | Branching Pattern in Plants[END_REF] in more detail.

With branching

The concentrations of nutrients C, and of hormones A and K are described by the diffusion equations with convective terms:

∂C ∂t + V ∂C ∂x = d C ∂ 2 C ∂x 2 , (2.5) ∂K ∂t -V K ∂K ∂x = d K ∂ 2 K ∂x 2 -µK, (2.6)
Figure 1: Functions f and g.

∂A ∂t -V A ∂A ∂x = d A ∂ 2 A ∂x 2 -µA. (2.7)
The convective speed V A in equations (2.6) and (2.7) can be different in comparison with equation (2.5). It corresponds to transport in the phloem in the direction from top (meristem) to bottom (root). The speed V in the first equation is determined as the speed of growth:

dL dt = V, L t=0 = L 0 , h dR dt = Cg(R) -σR, R t=0 = R 0 , V = f (R) (2.8)
Here d C , d K , d A and µ are parameters; the space variable x is defined independently for each branch.

The boundary conditions for C are the same as in the case without branching:

C x=0 = C 0 , d C ∂C ∂x + Cg(R) x=L = 0.
(2.9)

The boundary conditions for K describe its possible production in the root, and its production in the meristem with rate proportional to the rate of growth:

K x=0 = K 0 , d K ∂K ∂x x=L = εV. (2.10)
Finally, the boundary conditions for A are similar, except that the boundary condition at x = 0 takes into account that this horomone can be transported from the stem to the root:

∂A ∂x -βA x=0 = 0, d A ∂A ∂x x=L = εV. (2.11)
We define next the branching conditions. A new branch appears at x = x 0 and t = t 0 if

A(x 0 , t 0 ) = A b , K(x 0 , t 0 ) = K b , (2.12)
where A b and K b are some given values. Appearance of a new branch means that there is an additional interval connected to the previous one at its point x 0 . The variables C n , A n , K n and R n are described at the new interval by the same equations as above. Here the subscript n determines the number of the branch. We should complete the formulation by the initial value of the concentration R n = R n (t 0 ). It cannot be found as a solution of the problem but should be considered as a parameter. There are additional branching conditions for the concentrations

C n (0, t) = C(x 0 , t), A n (0, t) = A(x 0 , t), K n (0, t) = K(x 0 , t),
which means that the concentrations are continuous at branching points. We also need some conditions on the fluxes to provide the conservation of mass. Under the notations shown in Figure 2, we have

SV = S Ṽ + S n V n , S ∂C ∂x = S ∂ C ∂x + S n ∂C n ∂x .
Here S, S, and S n are parameters determined by the cross section area of the corresponding branch. If all branches have the same cross section, then S = S = S n = 1. If the cross section areas are conserved and narrow branches have the same diameter, then S = 1, and S = S n = 1/2. In this work we restrict ourselves to the case where all branches have the same cross section. Otherwise, the diameter of each branch should be considered as a function of time. The diameter would depend on fluxes of metabolites coming through the branch to the root. On the other hand, the fluxes of nutrients going through the branch from the root would depend on its cross section. We obtain a very complex time-dependent problem with many feedbacks. This will essentially complicate the understanding of the mechanism of growth and the interpretation of the results.

As is discussed in the previous section, bud formation is accompanied by production of A and K. The angle of the new branches with respect to the previous one is given as a parameter. It is not related to light and photosynthesis which are not considered in this model.

Computational algorithm

In this section we discuss the numerical algorithm for the one-dimensional problem with branching. The model is based on the reaction-diffusion equations, which can be considered in the following form:

∂ϕ ∂t + u ∂ϕ ∂x = d ∂ 2 ϕ ∂x 2 -µϕ, (3.13)
where ϕ corresponds to the unknown variable (e.g. C, K), x is the space variable defined on each particular branch. Let us introduce a finite-difference mesh along any branch with the nodal coordinates x i (i = 1, . . . , I) (Figure 3a). Let the nodes x i be located uniformly with distance δx between them. We will denote the values of ϕ in the nodes by ϕ i . Consider an explicit upwind approximation for equation (3.13) first in the interior nodes (i = 2, . . . , I -1)

ϕ n i -ϕ n-1 i δt + (|u| + u) ϕ n i+1 -ϕ n i 2δx + (|u| -u) ϕ n i-1 -ϕ n i 2δx = d ϕ n i+1 -2ϕ n i + ϕ n i-1 (δx) 2 -µϕ n i , (3.14)
where δt is the time step and δx the space step. The superscript n corresponds to the time step. To simplify the notations, in what follows we will omit n, that is, we will write ϕ i instead of ϕ n i , and will keep the notation ϕ n-1 i . Let us write equation (3.14) in the form

a i ϕ i-1 -c i ϕ i + b i ϕ i+1 = -f i , i = 2, . . . , I -1, (3.15)
where

a i = |u| -u 2δx + d (δx) 2 , b i = |u| + u 2δx + d (δx) 2 , c i = 1 δt + |u| δx + 2d (δx) 2 + µ, f i = ϕ n-1 i δt . (3.16)
Consider now the approximation of the boundary conditions. If the ends of the branch correspond to the root (x = 0) or to the growing part (Figure 4a), then the approximation of the boundary conditions here yields the expressions:

-c 1 ϕ 1 + b 1 ϕ 2 = -f 1 , (3.17)
and

a I ϕ I-1 -c I ϕ I = -f I . (3.18)
Specific expressions for a 1 , b 1 , . . . , f I in the case of some particular boundary conditions can be given in a similar way.

In the case of the growth of a single branch, the numerical scheme is completely described by system (3.15), (3.17), and (3.18). The matrix of this system is tri-diagonal:

             -c 1 b 1 0 0 0 0 0 . . . 0 a 2 -c 2 b 2 0 0 0 0 . . . 0 . . . 0 . . . 0 a i -c i b i 0 . . . 0 . . . 0 . . . 0 0 0 0 a I-1 -c I-1 b I-1 0 . . . 0 0 0 0 0 a I -c I              . (3.19)
In order to solve a linear system with a tri-diagonal matrix (3.19), we can use the Thomas algorithm.

The situation is more difficult if there are several branches. Moreover, the tree can have an arbitrary structure (see, e.g., Figure 3b). There are nodes from which several branches end or begin. Such nodes will also be called boundary nodes. Each of them will be numbered, and the value of ϕ at the k-th node is denoted by ϕ k . The boundary conditions in the conventional sense are not imposed in such nodes. Instead, we should take into account conservation of species and fluxes.

Let there exist J branches leaving the boundary node k. Then ϕ j,1 = ϕ k at j = 1, . . . , J. The values of ϕ at the second node of each branch will be denoted by ϕ j,2 , at the third ϕ j,3 , and so on. Then the conservation of ϕ at the node k can be written as follows:

J j=1 d ϕ j,2 -ϕ k δx -µ δx 2 ϕ k - δx(ϕ k -ϕ n-1 k ) 2δt - |u j | -u j 2 (ϕ k -ϕ j,2 ) = 0. (3.20)
If there are J branches entering the boundary node k, then the finite difference scheme at this node will be:

J j=1 d ϕ j,I-1 -ϕ k δx -µ δx 2 ϕ k - δx(ϕ k -ϕ n-1 k ) 2δt - |u j | + u j 2 (ϕ k -ϕ j,I-1 ) = 0. (3.21)
If for some of the branches the numeration begins at the node k and for some other it ends there, then the finite difference scheme can be represented as a superposition of (3.20) and (3.21).

To summarize, we recall that, in the internal nodes of the mesh, we use a three-point scheme (3.15). In the boundary nodes, we have expressions (3.20) or (3.21) (see Figure 4b) that include ϕ k and the values ϕ j,I-1 or ϕ j,2 taken at the nodes nearest to the boundary node from the branches connected to this node. Thus we have a completely implicit approximation of our 1D problem.

We describe here the algorithm for the case in which the domain of the computation has several branches arbitrarily connected between each other. Moreover, the number of branches can change during the computation.

The finite difference scheme for the internal nodes of the j-th subsystem has the following form:

               -c 1 ϕ k 1 + b 1 ϕ 2 = -f 1 + other terms, a i ϕ i-1 -c i ϕ i + b i ϕ i+1 = -f i , i = 2, . . . , I -1 a I ϕ k 2 -c I ϕ I-1 = -f I + other terms. (3.22)
The first and the last equations in system (3.22) are obtained from expressions (3.20) and (3.21) corresponding to k 1 -th and k 2 -th boundary nodes. The terms corresponding to the j-th branch are written explicitly; for other branches they are included in "other terms".

Instead of the tri-diagonal matrix (3.19) we have here an "almost tri-diagonal" matrix for the j-th branch:

             -c 1 b 1 0 0 0 0 0 . . . 0 other terms a 2 -c 2 b 2 0 0 0 0 . . . 0 . . . 0 . . . 0 a i -c i b i 0 . . . 0 . . . 0 . . . 0 0 0 0 a I-1 -c I-1 b I-1 0 . . . 0 0 0 0 0 a I -c I other terms              . ( 3.23)
Step 1. Transformation of matrix (3.23) to an "almost two-diagonal" form. Let us write the corresponding system (3.22) in the form

               -c 1 ϕ k 1 + b 1 ϕ 2 = -f 1 + other terms, β i ϕ k 1 -ϕ i + α i ϕ i+1 = -γ i , i = 2, . . . , I -1, a I ϕ k 2 -c I ϕ I-1 = -f I + other terms, (3.24) where α i = b i c i -α i-1 a i , β i = a i β i c i -α i-1 a i , γ i = a i γ i-1 + f i c i -α i-1 a i , i = 3, . . . , I -1 (3.25) and α 2 = b 2 c 2 , β 2 = a 2 c 2 , γ 2 = f 2 c 2 . ( 3.26)
Thus, matrix (3.23) is reduced to the following one:

             -c 1 b 1 0 0 0 0 0 . . . 0 other terms β 2 -1 α 2 0 0 0 0 . . . 0 . . . β i . . . 0 0 -1 α i 0 . . . 0 . . . β I-1 . . . 0 0 0 0 0 -1 α I-1 0 . . . 0 0 0 0 0 a I -c I other terms              . (3.27)
Step 2. Transformation of the matrix of system (3.24) to an almost diagonal form. Let us write the internal part of system (3.24) in the form:

               -c 1 ϕ k 1 + b 1 ϕ 2 = -f 1 + other terms, p i ϕ k 1 -ϕ i + q i ϕ k 2 = -s i , i = 2, . . . , I -1, a I ϕ k 2 -c I ϕ I-1 = -f I + other terms, (3.28)
where

p i = α i p i+1 + γ i , q i = α i q i+1 , s i = α i s i+1 + β i , i = I -3, . . . , 2 (3.29)
and

p I-2 = γ I-2 , q I-2 = α I-2 , s I-2 = β I-2 . ( 3.30)
Now matrix (3.27) is reduced to the following one:

             -c 1 b 1 0 0 0 0 0 . . . 0 other terms p 2 -1 0 0 0 0 0 . . . q 2 . . . p i . . . 0 0 -1 0 0 . . . q i . . . p I-1 . . . 0 0 0 0 0 -1 q I-1 0 . . . 0 0 0 0 0 a I -c I other terms              . ( 3.31)
Step 3. Determination of the new values of ϕ at the boundary nodes. For any branch we substitute the second equation (i = 2) from (3.28) into the first one, and the (I -1)-th equation to the last one. We obtain

(b k 1 p 2 -c k 1 ) ϕ k 1 + b k 1 q 2 ϕ k 2 = -b k 1 s 2 -f k 1 + other terms, (a k 2 -c k 2 q I-1 ) ϕ k 2 -c k 2 p I-1 ϕ k 1 = c k 2 s I-1 -f k 2 + other terms. (3.32)
Let us consider the first and the last equations in (3.32). We see that only the boundary nodes enter these equations. "Other terms" now also contain only the corresponding boundary nodes. Therefore, the boundary nodes form a complete system which allows us to find the unknowns ϕ in all boundary nodes.

Step 4. Determination of new values in the internal nodes. When the values of ϕ are found in all boundary nodes, we can obtain new values of ϕ in the internal nodes for all branches using the corresponding internal equations from system (3.28).

This algorithm is applicable for an arbitrary connection of subsystems. We will finish this section with a short discussion of the sensitivity of the results with respect to numerical discretization. Convergence of the numerical results when we decrease the time and the space steps occurs when we consider finite and not very big time intervals. The situation is much more complex if we consider longtime behavior. Figure 4 shows final (stationary) solutions for different space steps. The absence of convergence of the results can be related to structural instability specific for this model and important for growth of biological organisms (see the discussion in [START_REF] Bessonov | Dynamic Models of Plant Growth[END_REF]).

Branching pattern

A typical example of plant growth in the 1D model with branching is shown in Figure 5. In the beginning of the evolution there is a single branch which grows with an approximately constant speed until the growth period is finished. The apical meristem is located at the upper end of the interval. The plant hormones auxin and cytokinin are produced there and are transported along the whole branch. It can be diffusive or convective transport. Cytokinin can also be produced in the roots. This is taken into account through the boundary condition at x = 0.

If at some point of the branch the concentrations of auxin and of cytokinin take on some prescribed values, then a new bud appears. In the simulation shown in Figure 5, there are five buds that appear one after another at an approximately equal distance. Each of the buds contains its own apical meristem with some value of the GM-factor R. When a new bud appears, the initial value R 0 is prescribed. It can be some given constant or it can depend on some factors (on plant hormones or on the value of R in the apical meristem of the main growing branch). After that, the value of R i in the bud evolves according to the same equation (see (2.4)).

13

When R i becomes larger than R f , proliferation begins and a new branch starts growing from the bud. As it is shown in [START_REF] Bessonov | Dynamic Models of Plant Growth[END_REF], growth of branches can be stepwise, that is a branch has periods of growth and of rest during which the speed of its growth equals zero. When the main branch stops growing, under appropriate conditions a new branch can appear from another bud.

As it is discussed above, we consider diffusive and convective fluxes of nutrients. Diffusive flux acts throughout all branches, whereas convective flux is directed to growing branches. This is related to the continuity equation for the incompressible fluid. The routes of convective flux of nutrients is shown in red in Figure 5.

Thus, initiation of growth of new branches is determined by the interplay between the concentrations of nutrients and of the growth factor. In the example presented in Figure 5 there are five buds formed on the main branch. Three of them give branches of the second generation with three new buds on each of them. Branches of the second generation appear one after another when the main branch stops growing. It is interesting to note that one of the buds on the main branch gives a rudimentary branch which stops growing right after it appears. There are also some branches of the third generation. Similar to the 1D case without branching, the final plant size decreases for larger h. In the case with branching, the final length is determined basically by the number of branches and not by their length (Figure 6). There are four generations of branches for h = 0.0002, three of them for h = 0.0005 and h = 0.001, and only two generations for greater values of h. A possible explanation of the influence of h on the final length from the point of view of nonlinear dynamics is given in [START_REF] Bessonov | Dynamic Models of Plant Growth[END_REF].

Plant evolution in time can be influenced by the initial value of the GM-factor concentration in a new bud. We have described it in the case where R 0 = 0.12 (Figure 5). In fact, it is the same for all values of R 0 between 0.01 and 0.12. Further increase of this parameter changes the plant evolution (Figure 7). When the first bud appears, it does not stay dormant but gives a new branch right away. It grows at the same time with the main branch. The difference between the two cases is determined by the behavior of solutions of equation (2.4). As we have already discussed, when a new bud appears, we prescribe it an initial value R 0 of the GM-factor concentration. We recall that the concentration of nutrients C is a function of space and time. Its value C(x 0 , t 0 ) at the new bud determines the right-hand side of this equation: it equals zero at R = 0 and at two positive values of R, it is negative between first two zeros and positive between the second two. If the initial value R 0 is at the interval of negativity, then the solution rapidly decreases, and the bud remains dormant. When the main branch stops growing, the concentration of nutrients increases. When its value C(x 0 , t) at the bud approaches 1, the right-hand side of equation (2.4) becomes positive for small positive values of R. The concentration of the GM-factor starts growing. After some time it can reach the critical value which is necessary for the new branch to grow. However, it may happen that another branch will start growing before this one. Then the concentration of nutrient can drop down again, and the concentration of the GM-factor may also decrease. Formation of new buds and growth of new branches is determined by a complex interaction of plant hormones, nutrients and mitotic factors. The plant structure depends on the values of the parameters. Figure 8 shows the simulations where the buds are double and give symmetric branches. The symmetry of the plant growth is prescribed by the algorithm.

Figure 2 :

 2 Figure 2: Schematic representation of convective and diffusion fluxes in 1D plant with branching.

Figure 3 :

 3 Figure 3: Numerical mesh with branching

Figure 4 :

 4 Figure 4: The final (stationary) form of the plant for different space steps: for each figure, the step is half of the previous one.

Figure 5 :

 5 Figure 5: Time evolution of the plant. The figure shows the moments in time when new buds or branches appear. Convective flux of nutrients is shown in red.

Figure 6 :

 6 Figure 6: Final plant forms for different values of h.

Figure 7 :

 7 Figure 7: Evolution of the plant structure in time, h = 0.0005, R 0 = 0.13. Increase of R 0 changes completely the plant structure.

Consumption of nutrients and production of plant hormones are different in this case in comparison with the case without symmetry. Hence the number of new buds and the whole plant organization will also be different.

Figure 8: Plant evolution in the case of a double bud and R 0 = 0.12.