Direct and indirect approximations to positive solution for a nonlinear reaction-diffusion problem. Part II. Indirect approximation

7 years ago

Abstract

In this paper, we continue the study initiated in our previous work [3] and design a projection-like algorithm to approximate a hyperbolic unstable “point”. This “point” is in fact the positive solution of the reaction-diffusion problem considered in [3] and the algorithm modifies a finite difference (Euler)–finite elements scheme by incorporating the independence of the length of the domain condition. The numerical results are in good agreement with those obtained by direct methods as well as with those reported in [2], where the problem is solved in a Hamiltonian setting. At the same time we improve our previous results reported in [4]

Authors

C. I. Gheorghiu
Tiberiu Popoviciu Institute of Numerical Analysis,

D. Trif
Babes-Bolyai University, Faculty of Mathematics and Computer Science

Keywords

Nonlinear reaction-diffusion, positive solution, conserved integral,

projection-like method, f.e.m., finite elements–finite differences method, nonlinear stability, energetic method.

Paper coordinates

C.I. Gheorghiu, D. Trif, Direct and indirect approximations to positive solution for a nonlinear reaction-diffusion problem.
Part II. Indirect approximation, Rev. Anal. Numér. Théor. Approx. 31 (2002) 163-170.

PDF

https://ictp.acad.ro/jnaat/journal/article/view/2002-vol31-no2-art6/Gheorghiu-Trif-pdf

About this paper

Journal

Rev. Anal. Numer. Theor. Approx.

Publisher Name

Editura Academiei Romane

Paper on journal website

https://ictp.acad.ro/jnaat/journal/article/view/2002-vol31-no2-art6

Print ISSN

1222-9024

Online ISSN

2457-8126

MR

?

ZBL

?

Google Scholar Profile

google scholar link

[1] Beyn, W.-J., On the numerical approximation of phase portraits near stationary points, SIAM J. Numer. Anal., 24, pp. 1095–1113, 1987.

[2] Gheorghiu, C. I. and Muresan, A., On the significance of integral properties of orbits in some superlinear fixed-period problems, Proc. ICNODEA, Cluj-Napoca, 2001.

[3] Gheorghiu, C. I. and Trif, D., Direct and indirect approximations to positive solutions for a nonlinear reaction-diffusion problem, I. Direct (variational) approximation, this journal, 31 no. 1, pp. 61–70, 2002.

[4] Gheorghiu, C. I. and Trif, D., The numerical approximation to positive solution for some reaction-diffusion problems, Pu.M.A., 11 , pp. 243–253, 2001.

[5] Henry, D., Geometric theory of semilinear parabolic equations, Lecture Notes in Mathematics, 840, Springer-Verlag, 1981.

[6] Iserles, A., A First Course in the Numerical Analysis of Differential Equations, Cambridge Univ. Press, 1996.

[7] Larsson, S. and Sanz-Serna, J. M., The behaviour of finite element solutions of semi-linear parabolic problems near stationary points , SIAM J. Numer. Anal., 31 , pp. 1000– 1018, 1994.

[8] Stuart, A. M. and Humphries, A. R., Dynamical systems and numerical analysis, Cambridge Monographs on Applied and Computational Mathematics, Cambridge Univ. Press., 1996

Paper (preprint) in HTML form

jnaat,+Journal+manager,+2002-2-Gheorghiu-Trif-15-10-30

DIRECT AND INDIRECT APPROXIMATIONS TO POSITIVE SOLUTION FOR A NONLINEAR REACTION-DIFFUSION PROBLEM. II. INDIRECT APPROXIMATION

C. I. GHEORGHIU* and D. TRIF $^{†}$

Abstract

In this paper, we continue the study initiated in our previous work 3 and design a projection-like algorithm to approximate a hyperbolic unstable "point". This "point" is in fact the positive solution of the reaction-diffusion problem considered in 3 and the algorithm modifies a finite difference (Euler)-finite elements scheme by incorporating the independence of the length of the domain condition. The numerical results are in good agreement with those obtained by direct methods as well as with those reported in 2, where the problem is solved in a Hamiltonian setting. At the same time we improve our previous results reported in 4 .

MSC 2000. 65M20, 65M12, 37M99, 65P40.
Keywords. Nonlinear reaction-diffusion, positive solution, conserved integral, projection-like method, f.e.m., finite elements-finite differences method, nonlinear stability, energetic method.

1. INTRODUCTION

The main aim of this paper is to approximate the positive solution

\bar{u} (x)

to [GT]-(2) as the stationary unstable hyperbolic "point" of continuous dynamical system [GT]-(1) from [GT]

^{1}

. Thus, in the spirit of [8, Ch. 8], we refine the finite elements-finite differences (Euler) algorithm from our previous work [4] and prove its energetic stability. We get a projection-like method, which incorporates the conserved integral [GT]-(3).

Thus, the section 2 is concerned with the finite elements-finite difference (Euler) approximation

{\bar{u}}_{h} (x, t_{n})

of the solution

\bar{u} (x, t)

of [GT]-(1). The numerical approximation

{\bar{u}}_{h} (x, t_{n})

is the solution of the discrete analogous (2) of [GT]-(1). It is of highest importance to observe that whenever the condition [GT]-(3) is enforced at the initial and all subsequent time levels

t_{n}, n = 0, 1, \dots

, the modified discrete solution

{\bar{u}}_{h}^{mod} (x, t_{n})

approaches

{\bar{u}}_{h} (x)

. In any other case

{\bar{u}}_{h} (x, t_{n})

will approach in time the null solution of [GT]-(2) or will be blowingup.

With this in mind, we consider in Section 3 the role of condition [GT]-(3) in the numerical approximation of

\bar{u} (x)

. We will use essentially the complementary works of Larsson and Sanz-Serna [7] and Beyn [1]. They both refer to the approximation of the phase portrait of a dynamical system near a stationary hyperbolic point.

The first two authors show that the local stable and unstable manifolds of the spatial discrete problem converge to their continuous counterparts. Therefore, the spatial discretized dynamical system has the same qualitative behavior near

\bar{u}

as the continuous system. They also obtain error bounds of optimal order in the

L_{2}

and

H^{1}

norms which hold uniformly over arbitrarily long time intervals.

The second work, restricted to ordinary differential equations, shows that the phase portrait of a such dynamical system near a stationary hyperbolic point is reproduced correctly by numerical one-step methods i.e. any continuous trajectory can be approximated by an appropriate discrete trajectory and vice versa. In particular, the stable and unstable manifolds of the discretized system converge to their continuous counterparts.

As is apparent by Liapunov-Schmidt arguments, the stationary point

\bar{u} (x)

and its discrete approximation

{\bar{u}}_{h} (x, t_{n})

are hyperbolic (i.e. 0 is not one of the eigenvalues ) and unstable (i.e. there is at least one positive eigenvalue). Therefore by discretization in space and next in time, and then connecting the above results of Larsson & Sanz-Serna and Beyn, any continuous trajectory is approximated (shadowed) by an appropriate discrete trajectory and vice versa, as long as they remain in a sufficiently small neighborhood of

\bar{u} (x)

.

A discrete analogue of [GT]-(3), thought of as a quadrature formula, is used to "keep" a discrete modified solution

{\bar{u}}_{h}^{m o d} (x, t_{n})

near the stable manifold in a sufficiently small neighborhood of

\bar{u} (x)

. Figure 1 depicts this situation. Therefore, due to the condition [GT]-(3), the discrete modified solution

{\bar{u}}_{h}^{mod} (x, t_{n})

approximates the unstable equilibrium

\bar{u} (x)

as

t_{n} \to \infty

. The stability of the algorithm is proved in Section 4.

2. NUMERICAL APPROXIMATION OF $\bar{U} (X, T)$ BY A PROJECTION-LIKE ALGORITHM

In this section we consider an indirect method to find the positive solution of the problem [GT]-(2). Precisely, we use the finite elements-finite difference (Euler) approximation

{\bar{u}}_{h} (x, t_{n})

of the solution

\bar{u} (x, t)

of the problem [GT]-(1) and we look for the stationary solutions of this problem. The approximation, with the usual finite elements shape functions

φ_{k} (x)

, reads

u_{h} (x, t) = \sum_{k = 1}^{N} c_{k} (t) φ_{k} (x)

and it must be a solution of the semidiscretized system

\int_{0}^{1} \frac{\partial u_{h} (x, t)}{\partial t} φ_{k} (x) d x = \int_{0}^{1} (- \frac{\partial u_{h} (x, t)}{\partial x} \frac{\partial φ_{k} (x)}{\partial x} + u_{h}^{3} (x, t) φ_{k} (x)) d x

for

k = 1, \dots, N

. In the matrix form, the above system is

\begin{matrix} (1) & M C^{'} = F (C) \end{matrix}

where

M = {(m_{i k})}_{i, k = 1, \dots, N}, m_{i i} = 2 h / 3, m_{i k} = h / 6

for

| i - k | = 1, i \neq k

and

m_{i k} = 0

for

| i - k | > 1, C (t) = (c_{1} (t), \dots, c_{N} (t))

and

F

is defined in [GT]-(19), with

μ = 0

.

With an one step method, like Euler method, the system (1) becomes

\begin{matrix} (2) & C ((n + 1) τ) = C (n τ) + τ M^{- 1} F (C (n τ)), n = 0, 1, \dots \end{matrix}

and the starting approximation for

n = 0

is

C (0) = {(u_{0} (x_{k}))}_{k = 1, \dots, N}

from the Cauchy condition of [GT]-(1). The system (2) is the completely discretized form of the problem [GT]-(1). The above approximation, at the time level

t_{n} = n τ, τ > 0

, becomes

\begin{matrix} (3) & u_{h} (x, t_{n}) = \sum_{k = 1}^{N} c_{k} (t_{n}) φ_{k} (x) \end{matrix}

Obviously,

{\bar{u}}_{h} (x_{k}, t_{n}) = c_{k} (t_{n})

for

k = 1, \dots, N, n = 0, 1, \dots

and

t_{n} = n τ

.
With

u_{0} (x) = 6 π / \sqrt{2} x (1 - x)

as an example, for

N = 19

and

τ = h^{2} / 6

,

u_{h} (x, t_{n})

seems to be convergent to

u = 0

when

t_{n} \to \infty

. The same phenomenon occurs with other initial data, even they satisfy the integral condition [GT]-(3). The divergence

({\bar{u}}_{h} \to \infty

, for large values of

t)

of the successive approximations

u_{h} (x, t_{n})

also happens even if

u_{0} (x)

satisfies [GT]-(3). This fact confirms the analytic results of Henry [5, Ch. 3].

We remark that Euler algorithm is only conditionally stable and thus the time step

τ

must be quite small. Another unconditionally stable algorithm, like backward Euler method, will avoid this difficulty but we do not discuss here this issue.

If we modify the algorithm (2) such that at every time level

t_{n}

,

\begin{matrix} (4) & C^{mod} (t_{n}) = α C (t_{n}) \end{matrix}

where

α

is determined from the integral condition [GT]-(3), namely

\int_{0}^{1} \sum_{k = 1}^{N} c_{k}^{m o d} (t_{n}) φ_{k} (x) d x = \frac{π}{\sqrt{2}}

the algorithm becomes convergent to the unique positive solution

{\bar{u}}_{h} (x)

of the problem [GT]-(2). Therefore

lim_{n \to \infty} {\bar{u}}_{h}^{mod} (x_{k}, t_{n}) = lim_{n \to \infty} c_{k}^{mod} (t_{n})

is an approximation of the unique stationary positive solution

\bar{u} (x_{k})

, for

k = 1, \dots, N

. This is in fact our projection-like algorithm.

3. THE ROLE OF THE CONSERVED INTEGRAL CONDITION

In order to explain the role of the integral condition [GT]-(3) on the algorithm (4), we study the stability of the positive stationary solution of the problems [GT]-(2) and [GT]-(18).

Let us consider the linearization of the problem [GT]-(2), around

\bar{u}

,

\begin{aligned} - v^{''} - 3 {\bar{u}}^{2} v = 0, x \in [0, 1] \\ (5) & v (0) = v (1) = 0 \end{aligned}

where

\bar{u}

is the unique positive solution of [GT]-(2). There exists a real and increasing to

+ \infty

sequence of eigenvalues, and this sequence contains but only one negative eigenvalue. Indeed, the smallest eigenvalue is

λ_{1} = min {\int_{0}^{1} (v^{' 2} - 3 {\bar{u}}^{2} v^{2}) d x ∣ v \in H_{0}^{1} (0, 1), \int_{0}^{1} v^{2} d x = 1}

For

v = α \bar{u}

, with

α

so that

\int_{0}^{1} v^{2} d x = 1

, we have

λ_{1} \leq α^{2} \int_{0}^{1} ({\bar{u}}^{' 2} - 3 \bar{u} {\bar{u}}^{3}) d x = α^{2} \int_{0}^{1} ({\bar{u}}^{' 2} + 3 \bar{u} {\bar{u}}^{''}) d x = - 2 α^{2} \int_{0}^{1} {\bar{u}}^{' 2} d x < 0

In our case

λ = 0

is not an eigenvalue. Therefore,

\bar{u}

is a hyperbolic unstable stationary solution. In this case, Larsson and Sanz-Serna [7] show that the discretized in space (by piecewise linear finite elements) dynamical system (1) has the same qualitative behavior near

\bar{u}

as the continuous system. This means that there is an

H^{1}

neighborhood

V

of

\bar{u}

, such that for any exact solution

u (x, t)

contained in

V

for

T_{0} < t < T_{1}

there is a numerical solution

u_{h} (x, t)

which approximates

u

accurately for

T_{0} < t < T_{1}

and vice versa. Moreover, the discrete problem has local stable and unstable manifolds that converge to their continuous counterparts.

Let us consider now the linearization of the semidiscretized problem (1),

\begin{matrix} (6) & Φ^{'} = M^{- 1} d F (C^{*}) Φ \end{matrix}

where

C^{*}

is an approximation of the stationary solution obtained by the modified algorithm, with

t_{n} = 0.5

. The eigenvalues of

M^{- 1} d F (C^{*})

, computed with MATLAB are real and, excepting one, all negative. Table 1 contains the greatest three eigenvalues for some values of

N

and shows that

C^{*}

is an approximation of a hyperbolic unstable stationary solution

\bar{C}

.

Table 1. The greatest eigenvalues of the linearized semidiscretized system.

$N$	$λ_{1}$	$λ_{2}$	$λ_{3}$
99	20.6265	-20.6370	-70.7046
199	20.6256	-20.6284	-70.6563
299	20.6257	-20.6268	-70.6474
399	20.6256	-20.6263	-70.6442

In this case, Beyn 1 shows that the phase portrait of the dynamical system (1) near

\bar{C}

is reproduced correctly by an one-step numerical method and, in particular, the stable and unstable manifolds of the complete discretization

C (t_{n}), n = 1, 2, \dots

converge to their continuous counterparts.

The integral condition [GT]-(3) represents a hyperplane in

L_{2} (0, 1)

which is orthogonal to the constant function 1. The intersection of this hyperplane with the finite-dimensional subspace

H_{N}

of

H_{0}^{1} (0, 1)

spanned by

φ_{1}, \dots, φ_{N}

, is a hyperplane (

H

) in

H_{N}

. In

H_{N}

there exist local stable manifold and local unstable manifold of the system (1), tangent to the eigenspaces

E_{S}

and

E_{U}

of the linearized system (6) in the stationary discrete solution

{\bar{u}}_{h}

.

In our case, the eigenspace

E_{U}

is spanned by the eigenvector

v_{1}

corresponding to the eigenvalue

λ_{1}

in Table 1. The angle between the projection of the constant function 1 on

H_{N}

and

v_{1}

is about

30^{\circ}

and the angle between

v_{1}

and

{\bar{u}}_{h}

is about

5^{\circ}

. Consequently, in a neighborhood of

{\bar{u}}_{h}, (H)

lies near

E_{S}

and then, the modification (4) of the algorithm shifts the predicted value of the iterations

u_{n}^{α}

given by (2) on

u_{n}

, close to the stable manifold (see the heuristic Fig. 1).

Fig. 1. The projection-like algorithm in

H_{N}

.

The importance of the integral condition [GT]-(3) becomes again obvious: it is a way to approach the stable manifold and finally, to approach the unstable hyperbolic equilibrium

\bar{u}

.

4. THE STABILITY OF THE ALGORITHM

As a matter of fact, our problem [GT]-1 reads as follows

{\begin{cases} u_{t} & = u_{x x} + u^{3}, & x \in (0, L), t > 0 \\ u (0, t) & = u (L, t) = & 0 \\ u (x, 0) & = u_{0} (x), & x \in (0, L) \end{cases}

where

u_{0} (x)

satisfies

u_{0} (0) = u_{0} (L) = 0

. A semi-discretization (SD) scheme converts this initial-boundary value problem into a system of coupled ODEs. To proceed, we "stretch" grid points into a long vector, i.e. introduce the notation

u_{Δ x} := {(u_{1}, \dots, u_{N})}^{T}, where u_{i} := u (x_{i}, t), t > 0

and equip the underlying linear space with the Euclidian norm

{‖ u_{Δ x} ‖}_{Δ x} := {(Δ x \sum_{j = 1}^{N} {| u_{j} |}^{2})}^{1 / 2}

(the sum ranges over all space grid points). The (SD) scheme (the method of lines) applied to the above problem leads to the ODEs system

u_{l}^{'} = \frac{1}{(Δ x)^{2}} (u_{l - 1} - 2 u_{l} + u_{l + 1}) + u_{l}^{3}, l = 1, 2, \dots, N

As both, the eigenvalue technique and Fourier analysis, are out of question for such nonlinear system, we recourse to the energy method.

Differentiating

{‖ u_{Δ x} ‖}_{Δ x}^{2}

yields

\begin{aligned} \frac{1}{2} \frac{d}{d t} {‖ u_{Δ x} ‖}_{Δ x}^{2} & = \frac{1}{2} (Δ x) \frac{d}{d t} \sum_{l = 1}^{N} u_{l}^{2} \\ = Δ x \sum_{l = 1}^{N} u_{l} u_{l}^{'} \\ = \sum_{l = 1}^{N} [\frac{1}{Δ x} (u_{l} u_{l - 1} - 2 u_{l}^{2} + u_{l} u_{l + 1}) + (Δ x) u_{l}^{4}] \end{aligned}

At this stage we take into account the condition

[GT] - (3)

and write that as

\sum_{l = 1}^{N + 1} \frac{Δ x}{2} (u_{l - 1} + u_{l}) = \frac{π}{\sqrt{2}}

As we approximate a positive solution, this last equation implies

\begin{matrix} (7) & \sum_{l = 1}^{N} u_{l}^{2} \leq \frac{π^{2}}{2 (Δ x)^{2}} and \sum_{l = 1}^{N} u_{l}^{4} \leq \frac{π^{4}}{4 (Δ x)^{4}} . \end{matrix}

The (SD) scheme lumped together with discrete "isoperimetric" condition represents in fact an extension of our modified algorithm. We want to establish the stability of this algorithm.

First, we observe that

\frac{1}{2} \frac{d}{d t} {‖ u_{Δ x} ‖}_{Δ x}^{2} \leq \sum_{l = 1}^{N} \frac{1}{Δ x} | u_{l} u_{l - 1} - 2 u_{l}^{2} + u_{l} u_{l + 1} | + Δ x \sum_{l = 1}^{N} u_{l}^{4} .

Then, we resort to Cauchy-Schwartz inequality to deduce that

\frac{1}{2} \frac{d}{d t} {‖ u_{Δ x} ‖}_{Δ x}^{2} \leq \frac{4}{Δ x} {(\sum_{l = 1}^{N} u_{l}^{2})}^{1 / 2} {(\sum_{l = 1}^{N} u_{l}^{2})}^{1 / 2} + Δ x \sum_{l = 1}^{N} u_{l}^{4}

and use (7) to have

\frac{d}{d t} {‖ u_{Δ x} ‖}_{Δ x}^{2} \leq \frac{K}{(Δ x)^{3}}, K = 4 π^{2} + \frac{π^{4}}{2}

It follows at once that, for every

T > 0

, there exists a constant

c (T)

such that

{‖ u_{Δ x} ‖}_{Δ x} \leq \frac{1}{(Δ x)^{3 / 2}} c (T), t \in [0, T]

where

c (T) = {({‖ u_{Δ x} (0) ‖}_{Δ x}^{2} \times (Δ x)^{3} + K T)}^{1 / 2}

. For a fixed spatial grid, this is precisely what is required for stability 6 .

5. CONCLUSIONS AND PERSPECTIVES

It is now quite clear that the role of the integral condition [GT]-(3) is far beyond the elementary computation of the average temperature in a reaction diffusion process. In fact, using that, we have considered two coherent strategies, one direct (based on variational calculus) and one indirect (based on a completely discretized scheme) which lead to solutions to [GT]-(2). These solutions conserve the known qualitative behavior (structure) of the underlying differential system. This is the most important benefit.

Extensions of these strategies to multidimensional problems are in progress.

REFERENCES

[1] Beyn, W.-J., On the numerical approximation of phase portraits near stationary points, SIAM J. Numer. Anal., 24, pp. 1095-1113, 1987.
[2] Gheorghiu, C. I. and Mureşan, A., On the significance of integral properties of orbits in some superlinear fixed-period problems, Proc. ICNODEA, Cluj-Napoca, 2001.
[3] Gheorghiu, C. I. and Trif, D., Direct and indirect approximations to positive solutions for a nonlinear reaction-diffusion problem, I. Direct (variational) approximation, this journal, 31 no. 1, pp. 61-70, 2002. ©
[4] Gheorghiu, C. I. and Trif, D., The numerical approximation to positive solution for some reaction-diffusion problems, Pu.M.A., 11, pp. 243-253, 2001.
[5] Henry, D., Geometric theory of semilinear parabolic equations, Lecture Notes in Mathematics, 840, Springer-Verlag, 1981.
[6] Iserles, A., A First Course in the Numerical Analysis of Differential Equations, Cambridge Univ. Press, 1996.
[7] Larsson, S. and Sanz-Serna, J. M., The behaviour of finite element solutions of semilinear parabolic problems near stationary points, SIAM J. Numer. Anal., 31, pp. 10001018, 1994.
[8] Stuart, A. M. and Humphries, A. R., Dynamical systems and numerical analysis, Cambridge Monographs on Applied and Computational Mathematics, Cambridge Univ. Press., 1996.

Received by the editors: June 29, 2000.

*"T. Popoviciu" Institute of Numerical Analysis, P.O. Box 68, 3400 Cluj-Napoca 1, Romania (ghcalin@ictp.acad.ro).
† "Babeş-Bolyai" University, Faculty of Mathematics and Computer Science, st. M. Kogălniceanu 1, 3400 Cluj-Napoca, Romania, e-mail: dtrif@math.ubbcluj.ro.
$^{1}$ We will refer to our paper 3] as [GT] and (1), (2) and (3) are respectively the semilinear parabolic problem, the two-point boundary value problem and the conserved integral condition.

2002

Direct and indirect approximations to positive solution for a nonlinear reaction-diffusion problem. Part II. Indirect approximation

Abstract

Authors

Keywords

Paper coordinates

PDF

About this paper

Journal

Publisher Name

Paper on journal website

Print ISSN

Online ISSN

MR

ZBL

Google Scholar Profile

Paper (preprint) in HTML form

DIRECT AND INDIRECT APPROXIMATIONS TO POSITIVE SOLUTION FOR A NONLINEAR REACTION-DIFFUSION PROBLEM. II. INDIRECT APPROXIMATION

Abstract

1. INTRODUCTION

2. NUMERICAL APPROXIMATION OF U ¯ ( X , T ) U ¯ ( X , T ) bar(U)(X,T)\bar{U}(X, T)U¯(X,T) BY A PROJECTION-LIKE ALGORITHM

3. THE ROLE OF THE CONSERVED INTEGRAL CONDITION

4. THE STABILITY OF THE ALGORITHM

5. CONCLUSIONS AND PERSPECTIVES

REFERENCES

Related Posts

2. NUMERICAL APPROXIMATION OF $\bar{U} (X, T)$ BY A PROJECTION-LIKE ALGORITHM