16th Conference on Applied Mathematics, Univ. of Central Oklahoma, Electronic Journal of Differential Equations, Conf. 07, 2001, pp. 39–45. ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Approximating parameters in nonlinear reaction diffusion equations ^{∗}

Robert R. Ferdinand

Abstract

We present a model describing population dynamics in an environment. The model is a nonlinear, nonlocal, reaction diffusion equation with Neumann boundary conditions. An inverse method, involving minimization of a least-squares cost functional, is developed to identify unknown model parameters. Finally, numerical results are presented which display estimates of these parameters using computationally generated data.

1 Introduction

In [4] parameter estimation in a nonlinear reaction diffusion equation is discussed and numerical results are presented. In this paper, a similar model is considered with a non-singular, nonlocal diffusion term and Neumann boundary conditions. Model solution is approximated using a Galerkin approximation scheme using finite elements. Certain model parameters are then estimated using an inverse method procedure. Although a similar inverse problem has been considered in [1], the model there has a singular diffusion term, Dirichlet boundary conditions and the solution is estimated using a finite-difference scheme. Also, several inverse problems have also been presented and studied in [2, 3, 4, 5, 6]. In this paper, we consider the following initial boundary value problem which describes population dynamics in an environment:

\begin{matrix}u_{t} − a(l(u(·, t)))u_{xx} = h(u) + f(t) (t, x) ∈ Ω̂_{1} × Ω_{2} \\ u_{x}(t, 0) = 0 = u_{x}(t, x_{max}) t ∈ Ω_{1} \\ u(0, x) = u^{0}(x) x ∈ Ω_{2}. \end{matrix}  (1.1)

Here, l(u(·, t)) = ∫_{0}^{x_{max}}g(x)u(t, x)dx, with Ω_{1} = [0, T_{max}], Ω̂_{1} = (0, T_{max}] and

Ω_{2} = [0, x_{max}]. u(t, x) represents population density of organizm with size x at time t. a(l) is the non-singular, nonlinear, nonlocal diffusion term with kernel g(x) ∈ L^{2}(Ω_{2}). The logistic function h(u) = u(1 − u) and f(t) are the reaction terms. Neumann boundary conditions are given while initial condition u^{0}(x) represents initial population density.

∗Mathematics Subject Classifications: 65N21, 65N30, 65N12, 35K05, 35K55, 35K57.

Key words: Parameter estimation, inverse problem, Galerkin, reaction-diffusion equation.

\begin{matrix}c \\ ◯ \end{matrix}2001 Southwest Texas State University.

Published July 20, 2001.

This paper is organized in the following fashion. In Section 2, a least squares method is developed to identify unknown parameters a and f from observed data z_{i}. Numerical examples are presented in Section 3 which illustrate the estimation of these parameters from computationally generated data.

2 The Inverse Problem

In this section a numerical method is developed to solve the following infinite dimensional least squares problem: Given observations z(t_{i}, x) at times {t_{i}}_{i=1}^{K} with 0 ≤ t_{1} < t_{2} < ⋯ < t_{K} ≤ T_{max} and position x ∈ Ω_{2}, find a parameter q̂ ∈ Q (an infinite dimensional compact set to be specified later) which minimizes the performance index given by the following least-squares cost functional:

J(q) = Φ(u(·; q); z) = ∑_{i=0}^{K} ∫_{0}^{x_{max}} |u(t_{i}, x; q) − z(t_{i}, x)|^{2} dx,  (2.1)

where for each q ∈ Q, the notation u(q) = u(t, x; q) represents the parameter dependent solution of the model equation (1.1). This numerical method involves two levels of numerical approximation. The first approximates u, the solution to (1.1), while the second approximates the parameter space Q. Galerkin approximation is used to approximate the model solution. Similar methods have been used to approximate other heat flow models in [2, 4].

Starting with the first level, equation (1.1) is written in weak form as follows:

\begin{matrix}⟨u_{t}, ϕ⟩ + a(l(u)) ⟨u_{x}, ϕ_{x}⟩ = ⟨[h(u) + f(t)] , ϕ⟩ \\ u(0, x) = u^{0}(x). \end{matrix}  (2.2)

Now setting

u^{N}(t, x) = ∑_{i=0}^{N}w_{i}^{N}(t)ϕ_{i}^{N}(x)

in (2.2), where {ϕ_{j}^{N}}_{j=0}^{N} represent linear B-splines defined on a uniform partition 0 = x_{0} < x_{1} < x_{2} < ⋯ < x_{N} = x_{max} of [0, x_{max}], we arrive at

\begin{matrix}Λ^{N}ẇ^{N} (t) + G^{N} (w^{N}(t); a) = Υ^{N} (t, w^{N}(t)) \\ Λ^{N}w^{N} (0) = (w^{N})^{0} , \end{matrix}  (2.3)

where t ∈ [0, T_{max}] and w^{N}(t) = (w_{0}^{N}, w_{1}^{N}, ⋯, w_{N}^{N}) ∈ ℝ^{N+1}. Here,

• Λ^{N} is an (N + 1) × (N + 1) Gram matrix whose (i, j)^{th} entry is given by Λ_{i,j}^{N} = ⟨ϕ_{i}^{N}, ϕ_{j}^{N}⟩.

• (w^{N})^{0} is an (N + 1)-dimensional vector whose i^{th} element is given by (w^{N})_{i}^{0} = ⟨u^{0}, ϕ_{i}^{N}⟩.

• Furthermore, we have

Υ_{0}^{N} (t, w^{N}(t)) = 0
Υ_{i}^{N} (t, w^{N}(t)) = Δx [h(w_{i}^{N}(t)) + f(t)] for i = 1, ⋯, N

G_{0}^{N}(α; a) = a( ∑_{i=1}^{N} g(x_{i})w_{i}^{N}(t)Δx)(\frac{α_{0} − α_{1}}{Δx} )

G_{N}^{N}(α; a) = a( ∑_{i=1}^{N} g(x_{i})w_{i}^{N}(t)Δx)(\frac{α_{N} − α_{N−1}}{Δx} )

and for i = 1, ⋯, N − 1, α ∈ ℝ^{N+1},

G_{i}^{N}(α; a) = a( ∑_{i=1}^{N} g(x_{i})w_{i}^{N}(t)Δx)(\frac{−α_{i−1} + 2α_{i} − α_{i+1}}{Δx} ) .

Second level of the numerical scheme involves approximating the infinite dimensional parameter space by a sequence {Q^{M}} of finite dimensional spaces. Thus we estimate function a as a one-dimensional function of l, independent of t and x and estimate function f as one-dimensional function of t on these {Q^{M}}, where M = (M_{1}, M_{2}). This results in the following approximations for a and f:

(I_{M_{1}}a)(l) = ∑_{j=0}^{M_{1}}a(l_{a} + j(\frac{l̂_{a} − l_{a}}{M_{1}} ))ψ_{M_{1}}^{j}(l ; l_{a}, l̂_{a}),  (2.4)

where l ∈ ℝ and {ψ_{M_{1}}^{j}(l ; l_{a}, l̂_{a})}_{j=0}^{M_{1}} are linear B-splines defined on uniform

partition of interval [l_{a}, l̂_{a}]. Similarly,

(I_{M_{2}}f) (t) = ∑_{j=0}^{M_{2}}f(\frac{j}{M_{2}} T_{max})λ_{M_{2}}^{j} (t ; T_{max}) ,  (2.5)

where t ∈ [0, T_{max}] and {λ_{M_{2}}^{j}(t ; T_{max})}_{j=0}^{M_{2}} are linear B-splines defined on uniform partition of interval [0, T_{max}].

Thus in the finite space of dimension M = (M_{1}, M_{2}), solve the following initial-value problem:

\begin{matrix}Λ^{N}ẇ^{N} (t) + G^{N} (w^{N}(t); a_{M_{1}}) = Υ^{N} (t, w^{N}(t)) \\ Λ^{N}w^{N} (0) = (w^{N})^{0} . \end{matrix}  (2.6)

The definitions of all symbols in (2.6) are same as those in (2.3) with the following exceptions:

Υ_{0}^{N} (t, w^{N}(t)) = 0
Υ_{i}^{N} (t, w^{N}(t)) = Δx [h (w_{i}^{N}(t)) + f_{M_{2}}(t)] for i = 1, ⋯, N.

G_{0}^{N}(α; a_{M_{1}}) = a_{M_{1}}( ∑_{i=1}^{N} g(x_{i})w_{i}^{N}(t)Δx)(\frac{α_{0} − α_{1}}{Δx} )

G_{N}^{N}(α; a_{M_{1}}) = a_{M_{1}}( ∑_{i=1}^{N} g(x_{i})w_{i}^{N}(t)Δx)(\frac{α_{N} − α_{N−1}}{Δx} )

and for i = 1, ⋯, N − 1, α ∈ ℝ^{N+1},

G_{i}^{N}(α; a_{M_{1}}) = a_{M_{1}}( ∑_{i=1}^{N} g(x_{i})w_{i}^{N}(t)Δx) (\frac{−α_{i−1} + 2α_{i} − α_{i+1}}{Δx} ) .

Thus, for the sake of computations, we consider the following approximation to our infinite dimensional minimization problem defined in (2.1)

min_{q∈Q^{M}} J^{N} (q_{M}) = Φ(u^{N} (·; q_{M}) ; z) = ∑_{i=0}^{K} ∫_{0}^{x_{max}} |u^{N}(t_{i}, x; q_{M}) − z(t_{i}, x)|^{2} dx.

In several cases the space of linear splines can be taken to be Q^{M}.

In the next section, numerical examples are presented which shows our numerical scheme actually works.

3 Numerical Results

To numerically test the least-squares method, computational data z(t_{i}, x) is generated. The parameters in model equation (1.1) are chosen as follows:

f(t) = 2 + sin(100t), a(l) = 1/(l + 1)
h(u) = u(1 − u), g(x) = 1 + x^{2}, u^{0}(x) = x(1 − x).

Equation (1.1) is solved using Galerkin approximation described earlier and parameters given above. In these computations, Δt = 10^{−4}, Δx = 0.125 and the constants x_{max} and T_{max} are chosen as 1 and 0.05, respectively. Observations z(t_{i}, x_{j}) are then collected at points t_{i}, i = 0, ⋯, 250, where t_{i} = 0.0002 ∗ i and x_{j}, j = 0, ⋯, 8, where x_{j} = 0.125j.

All parameters are assumed to be known except a(l) and f(t). For an admissible parameter set Q, let D = C_{B} ([0, T] × ℜ), the space of bounded continuous functions on [0, T] × ℜ with the supremum norm. For fixed values of σ, â and l̂, Q is chosen as the D closure of the set A × F where

A = {a ∈ C(0, ∞), such that a′ < σ and a(l) = â for l > l̂}

Figure 1: Exact versus estimated f(t) = 2 + sin(100t) with M_{1} = 5

Figure 2: Exact versus estimated f(t) = 2 + sin(100t) with M_{1} = 11

and

F = {f ∈ C[0, T] such that |f(t_{1}) − f(t_{2})| ≤ L|t_{1} − t_{2}|} ,

where L is a Lipschitz constant for f. It is easily verifiable that Q is a compact subset of D. Further, looking at equation (2.4), lim_{M_{1}→∞} I_{M_{1}} (a) = a in C_{B} ([0, T] × ℝ), uniformly in a, for a ∈ Q. Also, from equation (2.5), lim_{M_{2}→∞} I_{M_{2}} (f) = f in C_{B} ([0, T] × ℜ), uniformly in f, for f ∈ Q (see [7]). Hence, if a_{M_{1}}(l) and f_{M_{2}}(t) are given by

a_{M_{1}}(l) = ∑_{j=0}^{M_{1}} ν_{M_{1}}^{j}ψ_{M_{1}}^{j} (l ; l_{a_{M_{1}}} , l̂_{a_{M_{1}}} )

and

f_{M_{2}} (t) = ∑_{j=0}^{M_{2}} β_{M_{2}}^{j}λ_{M_{2}}^{j} (t ; T_{max}) ,

Figure 3: Exact versus estimated a(l) = 1/(l + 1) with M_{2} = 5

Figure 4: Exact versus estimated a(l) = 1/(l + 1) with M_{2} = 11

respectively, then the least squares problem involves the identification of the (M_{1}+3) coefficients {ν_{M_{1}}^{j}}_{j=0}^{M_{1}}, l_{a_{M_{1}}} and l̂_{a_{M_{1}}} and (M_{2}+1) coefficients {β_{M_{2}}^{j}}_{j=0}^{M_{2}}, from a compact subset of ℜ^{M_{1}+M_{2}+4}. Initial guesses are as follows: ν^{j} = 1, j = 0, ⋯, M_{1}, β^{j} = 2, j = 0, ⋯, M_{2}, l_{a_{M_{1}}} = 0 and l̂_{a_{M_{1}}} = 1. The subroutine LMDIF1, obtained from NETLIB, is used in the computations. This FORTRAN software is an application of the Levenberg-Marquardt algorithm. Computations were executed on a SCO Unix 5.0.5 machine at East Central University, consisting of two 550 mhz Xeon processors in parallel. The first two figures which immediately follow, show a comparison between exact and estimated function f(t) for M_{1} = 5 and 11, respectively with dots and solid line representing estimated and exact function, respectively. It can be seen from these figures that as the value of M_{1} and subsequently the number of approximating elements increases, the estimated function gets closer to the exact function, thereby demonstrating a convergence to the minimizer of the original least squares cost functional. The last two figures show the same for function a(l) for M_{2} = 5 and 11, respectively . Values of the least squares cost functional at the end of the computer program range from between 10^{−13} to 10^{−8} in these experiments with execution time being approximately 2 minutes in each case, showing the numerical scheme to be working effectively.

Acknowledgment: The author wishes to thank his Ph.D. adviser Dr. Azmy Ackleh for his continued guidance and inspiration through the years.

References

[1] A. S. Ackleh (1999) “Parameter Estimation Problems for a Nonlinear Parabolic Equation with a Singular Nonlocal Diffusion Term”, Proceedings of the 7th IEEE Mediterranean Conference on Control and Automation, Haifa, Israel, 2191-2199.

[2] A. S. Ackleh, S. Aizicovici, R. R. Ferdinand and S. Reich (1999) “Parameter Identification in a Nonautonomous Nonlinear Volterra Integral Equation”, Proceedings of the 7th IEEE Mediterranean Conference on Control and Automation, Haifa, Israel, 2200-2206.

[3] A. S. Ackleh, S. Aizicovici and S. Reich (1998) “Parameter Identification in Nonlocal Nonlinear Evolution Equations”, Numerical Functional Analysis and Optimization, 21, 553-570.

[4] A. S. Ackleh, R. R. Ferdinand and S. Reich (1998) “Numerical Studies of Parameter Estimation Techniques in Nonlinear Evolution Equations”, Kybernetika, 34, 693-712.

[5] A. S. Ackleh and B. G. Fitzpatrick (1996) “Estimation of Discontinuous Parameters in General Nonautonomous Parabolic Systems”, Kybernetika, 32, 543-556.

[6] A. S. Ackleh and L. Ke (2000) “Existence-Uniqueness and Long Time Behavior for a Class of Nonlocal Nonlinear Parabolic Evolution Equations”, Proceedings of the American Mathematical Society, 128, 3483-3492.

[7] Banks H. T. and Kunisch K. (1989) “Estimation Techniques for Distributed Parameter Systems” (Birkhäuser, Boston).

[8] Schultz M. H. (1973) “Spline Analysis” (Prentice-Hall, Englewood Cliffs, New Jersey).

[9] Smoller J. (1994) “Shock Waves and Reaction Diffusion Equations” (Springer-Verlag, New York).

Robert R. Ferdinand Department of Mathematics, East Central University, Ada, OK 74820-6899 USA e-mail: robert.ferdinand@ecok.edu