Open AccessArticle

Inverse Problem Numerical Analysis of Forager Bee Losses in Spatial Environment without Contamination

Atanas Z. Atanasov

Miglena N. Koleva

^2,*

and

Lubin G. Vulkov

Agricultural Machinery Department, Agrarian and Industrial Faculty, University of Ruse, 8 Studentska Str., 7017 Ruse, Bulgaria

Department of Mathematics, Faculty of Natural Sciences and Education, University of Ruse, 8 Studentska Str., 7017 Ruse, Bulgaria

Department of Applied Mathematics and Statistics, Faculty of Natural Sciences and Education, University of Ruse, 8 Studentska Str., 7017 Ruse, Bulgaria

Author to whom correspondence should be addressed.

Symmetry 2023, 15(12), 2099; https://doi.org/10.3390/sym15122099

Submission received: 18 October 2023 / Revised: 9 November 2023 / Accepted: 17 November 2023 / Published: 22 November 2023

(This article belongs to the Special Issue New Trends on the Mathematical Models and Solitons Arising in Real-World Problems, 2nd Edition)

Download

Browse Figures

Versions Notes

Abstract

We consider an inverse problem of recovering the mortality rate in the honey bee difference equation model, that tracks a forage honeybee leaving and entering the hive each day. We concentrate our analysis to the model without pesticide contamination in the symmetric spatial environment. Thus, the mathematical problem is formulated as a symmetric inverse problem for reaction coefficient at final time constraint. We use the overspecified information to transform the inverse coefficient problem to the forward problem with non-local terms in the differential operator and the initial condition. First, we apply semidiscretization in space to the new nonsymmetric differential operator. Then, the resulting non-local nonsymmetric system of ordinary differential equations (ODEs) is discretized by three iterative numerical schemes using different time stepping. Results of numerical experiments which compare the efficiency of the numerical schemes are discussed. Results from numerical tests with synthetic and real data are presented and discussed, as well.

Keywords:

forage bee losses; spatial differential difference equation; final time constraint; inverse problem; ODEs system; non-local forward problem; iterative numerical algorithms

1. Introduction

Bees play a crucial role, not only for human beings but also for all plant species they assist in pollinating. Some of the most globally significant crops rely on bees for their reproduction, as these industrious pollinators visit and fertilize such plants. Furthermore, bees actively contribute to the preservation of natural ecosystems, as their pollination activities aid in the rejuvenation of trees. This, in turn, has a positive impact on the conservation of forest biodiversity and the maintenance of various other ecosystem services. In essence, bees serve as vital regulators of food production, forest equilibrium, and micro-climate dynamics

The loss of honey bee colonies is a prevalent occurrence that has significant economic and ecological consequences. The factors responsible for this phenomenon are currently under active investigation by many researches. Many studies have identified potential contributing factors to this problem, such as pesticides, parasites, climate change, nutritional challenges, etc. [1,2,3,4,5,6,7,8,9,10].

Remarkably, it is recognized that the development of mathematically tractable models is crucial for obtaining valuable insights into the ecological processes and factors that result in colony losses.

The pioneering mathematical model that integrated colony population dynamics influenced by biological and environmental factors was formulated in [11]. Subsequent honeybee mathematical models have been introduced, concerning the activities inside beehives, models investigating the influence of Varroa mites, models that simulate the dynamics of honeybees gathering nectar and pollen to support colony growth, models that gain valuable insights into ecological processes and factors contributing to colony losses, etc., see, e.g., review papers [2,12]. In all these models, the honeybee population dynamics are governed by ordinary differential equations (ODE).

Models involving partial differential equations to study honeybee colony dynamics are proposed in [13,14,15]. In [14,15], the authors model the thermoregulation process in honey bee colonies in winter using the Keller–Segel problem with a sign-changing chemotactic coefficient. The focus of the investigation in [13] is the exposure of forager bees to pesticides within contaminated spatial settings. The model comprises both differential and difference equations governing the spatial patterns of forager bees, distinguishing between those unaffected by contamination and those affected by it. An essential aspect of this model involves the daily return of forager bees to their hive.

Numerical simulations to study honeybee ODE problems are used in many papers. For example, in [16], critical hive sizes were determined numerically across various scenarios, emphasizing the significance of expedited forager recruitment in the depletion of hives during colony collapse. Both adult and immature honeybee populations along with their honey production are examined in [17]. For more details about the existing investigations, see also [2,12].

A second order in space, positivity-preserving numerical method for solving the Keller–Segel problem modeling the thermoregulation in honeybee colonies during the winter is developed in [18].

Over the past few decades, inverse problems have been employed to investigate numerous real-world phenomena and processes. This is primarily due to their robust mathematical formulation, comprehensive theoretical analysis, and effective numerical solutions [19,20,21,22,23,24,25,26,27].

The numerical method for solving the inverse bio-heat transfer problem for identifying space- and time-dependent perfusion coefficient from temperature measurements is constructed in [28]. The parameter identification inverse problem for a system of non-linear ODE modeling honeybee population dynamics taking into account different factors is studied in the literature; see [29,30] and reference therein. In [29], the authors consider a model describing the time evolution of the population count of hive bees, foraging workers, and contaminated foragers. A new problem for interaction of the food stock with the brood, adult bees and produced honey is suggested and analyzed numerically in [30].

Numerical identification approaches for predicting the optimal parameters governing the population dynamics in the beehive in the Keller–Segel problem under temperature measurements are developed in [31,32].

Inverse problems for recovering space-dependent mortality rate of the bees and the rate of contamination of the forager bees by pesticides under final time measurements are posed in [33].

In this work, we consider the inverse problem for recovering space-dependent reaction coefficient and solution, under final time measurements, in the Dirichlet parabolic problem. Such mathematical model is constructed in [13] to describe the spatial distribution of uncontaminated foraging bees, where the unknown solution is the bee density and the reaction coefficient represents the mortality rate of the bees. Determination of the honeybee mortality rate is of great importance for monitoring, studying the dynamics of the bee populations and developing good beekeeping practices.

We transform the inverse problem to a Dirichlet forward problem for non-local parabolic equation and non-local initial condition. Then, we further develop the iteration approach proposed in [34,35,36,37,38,39], where, in contrast to our problem, only initial conditions are non-local.

The remaining part of the paper is organized as follows. In the next section, we present the model problem. The inverse problem is formulated and discussed in Section 3. In Section 4, we construct three discrete schemes—implicit, Crank–Nicolson and averaging Saulyev’s scheme—for solving the 1D inverse problem, and on this base, we develop an iterative algorithm. In Section 5, we extend this iteration method for a 2D inverse problem. Computational results with synthetic and real data are presented in Section 6. Finally, in Section 7, we offer a summary and concluding remarks.

2. Model Problem

In this section, motivated by the results in [13], we formulate the model problem. This is a parabolic equation with Dirichlet boundary conditions that can be used to describe the spatial bee density of forage bees. We briefly discuss the model suggested in [13].

The basic assumption is that the environment is free from pesticide contamination. Forage bees depart from the hive at sunrise and come back to the hive after sunset on a daily basis.

We consider the spatial probability density function of forager bees in the hive

G (t^{*}, x, y)

for a given

t^{*}

\int_{R^{2}} \int G (t^{*}, x, y) d x d y = 1, for t^{*} > 0,

where G is Gaussian function

G (t^{*}, x, y) = \frac{1}{4 π ε t^{*}} exp (- \frac{{(x - x_{0})}^{2} + {(y - y_{0})}^{2}}{4 ε t^{*}})

(1)

with the standard derivation of G for a given value of

ε

and

t^{*}

defined by

σ = \sqrt{2 ε t^{*}},

and

(x_{0}, y_{0})

is the center of the hive.

Let us denote by

u (t, x, y)

the forager bee density at time t and location

(x, y) \in R^{2}

At the beginning of the first day, the initial distribution of forager bees is

u (0, x, y) = u_{0} (x, y) \geq 0, u_{0} \in L_{+}^{1} (R^{2}),

(2)

and the total number of forager bees (TNFB) at time

t = 0

is given by

U_{0} = \int_{R^{2}} \int u_{0} (x, y) d x d y .

Then, the forage bee density u satisfies equation

\begin{matrix} \frac{\partial u}{\partial t} = d ▵ u - μ u, ▵ u = \frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}, t \in [0, t), (x, y) \in R^{2}, \end{matrix}

(3)

where

U (t) = \int_{R^{2}} \int u (t, x, y) d x d y

(4)

is the the total number of forager bees at time

t \in [0, 1)

. Here, d denotes diffusion rate and

μ (x, y) \geq 0

is bee mortality rate. It encompasses homing failure and all other factors contributing to mortality in forager bees, and it is assumed in [13] that

μ (x, y)

is a bounded continuous nonnegative function on

R^{2}

At the conclusion of the initial day, we presume that all forager bees not affected by mortality make their way back to the hive, resulting in the total population at the end of the first day

U (1^{-}) = \int_{R^{2}} \int u (1^{-}, x, y) d x d y .

The forager bees that come back to the hive during the first day actively assist in caring for the hive’s young bees. The Allee function is used to characterize how they contribute to the emergence of forager bees that subsequently depart the hive for foraging:

\frac{β U {(1^{-})}^{2}}{X^{2} + U {(1^{-})}^{2}},

(5)

where

β

is the maximal production parameter and

X

is the sigmoidal Hill function production [40,41].

Among forager bees, two types of behaviors emerge as they return home at the end of the day and prepare for another day of foraging. A portion of foragers return to the hive, essentially starting anew the next day without any memory of their prior day’s whereabouts. We describe their foraging pattern on the subsequent day using a diffusion model. The second behavior involves foragers that remember productive foraging sites from the prior day. The next morning, these bees head directly to these remembered locations.

Integrating these two behavioral modes and under the assumption that forager bees spread within the hive according to the Gaussian probability density

G (t^{*}, x, y)

, we calculate

u_{1} (x, y)

. This function illustrates the spatial distribution of forager bees within the hive during the morning of the second day.

u_{1} (x, y) = G (t^{*}, x, y) (\frac{β U {(1^{-})}^{2}}{X^{2} + U {(1^{-})}^{2}} + (1 - q) U (1^{-})) + q u (1^{-}, x, y), q \in [0, 1],

where q is the fraction of forager bees that follow the second type of behavior. The portion of forager bees exhibiting this second behavior type consists of individuals with the ability to recall advantageous foraging sites from the prior day. When the next day begins, these bees immediately head towards these locations.

In this work, we consider Problems (2) and (3) defined in domain

{\bar{Q}}_{T} = \bar{Ω} \times [0, T]

\bar{Ω} = [0, X] \times [0, Y]

with Dirichlet boundary conditions

\begin{matrix} \begin{matrix} u (t, 0, y) = u_{l} (t, y), u (t, X, y) = u_{r} (t, y), y \in [0, Y], t \in (0, T], \\ u (t, x, 0) = v_{d} (t, x), u (t, x, Y) = v_{u} (t, x), x \in [0, X], t \in (0, T] . \end{matrix} \end{matrix}

(6)

3. The Inverse Problem

In this section, we discuss the posing of the inverse problem for recovering reaction space-dependent unknown coefficients

μ (x, y)

in parabolic Equation (3) with measurements at the final time

t = T

as it is performed in [33].

We consider Equation (3) at final time constraint

u (T, x, y) = φ (x, y), φ (x, y) > 0, (x, y) \in Ω,

(7)

where at first day

T = 1

, and at second day

T = 2

Frequently, despite the rapid advancements in electronic technology, beekeepers encounter challenges in obtaining adequate observational data to model forager bee behavior in food fields. We now utilize the observations stated in Equation (7) to approximately determine

u (t, x, y)

by solving the inverse problem.

Equation (3), equipped with known coefficient d and

μ

and provided with suitable initial and boundary conditions, is referred to as the direct problem.

The focus of the present work is to find

μ (x, y)

that satisfy Equation (3) and Observation (7). We consider the case of

v (t, x, y) = \frac{\partial u}{\partial t} (t, x, y), (t, x, y) \in {\bar{Q}}_{T} .

(8)

Differentiating Equation (3) with respect to t, we derive

\frac{\partial v}{\partial t} = d ▵ v - μ (x, y) v = 0 .

(9)

From Equation (3) and Observations (7), we find

μ (x, y) = \frac{- \frac{\partial u}{\partial t} (T, x, y) + d ▵ φ (x, y)}{φ (x, y)} = \frac{- v (T, x, y) + d ▵ φ (x, y)}{φ (x, y)} .

(10)

Thus, Problems (2), (3), (6) and (7) are equivalent to the following inverse problem, formulated as a direct problem for unknown solution

v (t, x, y)

(x, y) \in Ω = (0, X) \times (0, Y)

0 < t \leq T

\begin{matrix} \begin{matrix} \frac{\partial v}{\partial t} = d ▵ v + \frac{v (T, x, y) - d ▵ φ (x, y)}{φ (x, y)} v, \\ v (0, x, y) = d ▵ u_{0} (x, y) + \frac{v (T, x, y) - d ▵ φ (x, y)}{φ (x, y)} u_{0} (x, y), \end{matrix} \end{matrix}

(11)

with Dirichlet boundary conditions, derived from (6), applying (8),

\begin{matrix} \begin{matrix} v (t, 0, y) = v_{l} (t, y), v (t, X, y) = v_{r} (t, y), y \in [0, Y], t \in (0, T], \\ v (t, x, 0) = v_{d} (t, x), v (t, x, Y) = v_{u} (t, x), x \in [0, X], t \in (0, T] . \end{matrix} \end{matrix}

(12)

Using the results of [33] (Section 3), the following assertion can easily be proven.

Theorem 1.

Let us assume that

φ (x, y) \in C^{4 + γ} (Ω)

0 < γ < 1

and

φ (x, y) > 0

\bar{Ω}

. Then, there exists at least one solution of Problems (9)–(12).

Corollary 1.

Under the assumptions of Theorem 1, inverse Problems (3) and (7) have at least one solution.

Proof.

Indeed, we saw that if

u (t, x, y)

is a solution of Problems (3) and (7) with zero Dirichlet boundary conditions, then

v (t, x, y)

is a solution to Problems (11) and (12). Conversely, assuming that

u (t, x, y)

is a solution of Problems (3) and (7), it follows that

u (t, x, y) = u_{0} (x, y) + \int_{0}^{t} v (λ, x, y) d λ, μ (x, y) = \frac{- v (T, x, y) + d ▵ φ (x, y)}{φ (x, y)}

is a solution of inverse Problems (3) and (7) with zero Dirichlet boundary conditions. □

4. Iterative Method for Solving a 1D Inverse Problem

In this section, we propose three numerical iterative schemes for solving Problems (11) and (12). For clarity, first, we explain the method for the following one-dimensional prototype problem of (11) and (12):

\begin{matrix} \frac{\partial v}{\partial t} = d \frac{\partial^{2} v}{\partial x^{2}} + \frac{v (T, x) - d \frac{d^{2} φ}{d x^{2}}}{φ (x)} v, (x, t) \in Q_{T} = Ω \times (0, T], \\ v (t, 0) = v_{l} (t), v (t, x) = v_{r} (t), 0 < t \leq T, \\ v (0, x) = d \frac{d^{2} u_{0} (x)}{d x^{2}} + \frac{v (T, x) - d \frac{d^{2} φ}{d x^{2}}}{φ (x)} u_{0} (x), x \in Ω = (0, X) . \end{matrix}

(13)

For the solution of (13), we construct iterative finite difference methods using ideas from papers [34,36,37,38,39], where linear parabolic problems with non-local terms only in the initial condition are solved.

4.1. Discrete Schemes

We introduce uniform partition of the interval

[0, X]

. We let

h > 0

be small and

x_{i} = i h, i = 0, 1, \dots, I

, where

I = X / h

is a positive integer.

The semi-discrete solution

v_{i} (t) = v (x_{i}, t)

of (13) satisfies the following ODE system:

\begin{matrix} \begin{matrix} \frac{\partial v_{i}}{\partial t} = d ▵_{h} v_{i} + \frac{v_{i} (T) - d ▵_{h} φ_{i}}{φ_{i}} v_{i}, i = 1, 2, \dots, I - 1, \\ v_{0} = v_{l} (t), v_{I} = v_{r} (t), \\ v_{i} (0) = d ▵_{h} {(u_{0})}_{i} + \frac{v_{i} (T) - d ▵_{h} φ_{i}}{φ (x_{i})} {(u_{0})}_{i}, i = 1, 2, \dots, I - 1 . \end{matrix} \end{matrix}

(14)

In the time interval

[0, T]

, we define the uniform grid with nodes

t_{k} = k τ

k = 0, 1, \dots, K

τ = T / K

. The values of the solution and mesh functions at grid nodes are denoted by

w_{i}^{k} = w (x_{i}, t_{k})

Further, we apply three different temporal discretizations for the ODE system (14).

The implicit backward Euler scheme: we find ${v_{i}^{k}}$ , $k = 1, 2, \dots, K$ such that

$\begin{matrix} \begin{matrix} \frac{v_{i}^{k} - v_{i}^{k - 1}}{τ} = d ▵_{h} v_{i}^{k} + \frac{v_{i}^{K} - d ▵_{h} φ_{i}}{φ_{i}} v_{i}^{k}, i = 1, 2, \dots, I - 1, \\ v_{0}^{k} = v_{l} (t_{k}), v_{I}^{k} = v_{r} (t_{k}), \\ v_{i}^{0} = d ▵_{h} {(u_{0})}_{i} + \frac{v_{i}^{K} - d ▵_{h} φ_{i}}{φ (x_{i})} {(u_{0})}_{i}, i = 1, 2, \dots, I - 1, \end{matrix} \end{matrix}$

(15)

where $▵_{h} w_{i}^{k} = (w_{i + 1}^{k} - 2 w_{i}^{k} + w_{i - 1}^{k}) / h^{2}$ .
The Crank–Nicolson scheme. Now, ${v_{i}^{k}}$ , $k = 1, 2, \dots, K$ is defined by

$\frac{v_{i}^{k} - v_{i}^{k - 1}}{τ} = d ▵_{h} \frac{v_{i}^{k} + v_{i}^{k - 1}}{2} + \frac{v_{i}^{K} - d ▵_{h} φ_{i}}{φ_{i}} \frac{v_{i}^{k} + v_{i}^{k - 1}}{2}, i = 1, 2, \dots, I - 1,$

(16)

with the same initial and boundary approximate conditions as above.
The Saulyev-type alternating direction explicit (ADE) scheme. For the approximation, we use Saulyev’s first- and second-kind formula in order to obtain an unconditional stable numerical scheme, which can be realized in an explicit manner.

In discrete Scheme (15), we replace

▵_{h} v_{i}^{k}

▵_{h}^{1} v_{i}^{k} = (v_{i - 1}^{k} - v_{i}^{k - 1} - v_{i}^{k} + v_{i + 1}^{k - 1}) / h^{2}

and

▵_{h}^{2} v_{i}^{k} = (v_{i + 1}^{k - 1} - v_{i}^{k - 1} - v_{i}^{k} + v_{i - 1}^{k}) / h^{2}

. In order to obtain better accuracy, we combine these approximations and derive the Barakat and Clark scheme.

We find

{v_{i}^{k}}

at each time layer

k = 1, 2, \dots, K

, as an averaging solution of numerical schemes

\begin{matrix} \begin{matrix} \frac{v_{i}^{k} - v_{i}^{k - 1}}{τ} = d ▵_{h}^{1} v_{i}^{k} + \frac{v_{i}^{K} - d ▵_{h} φ_{i}}{φ_{i}} v_{i}^{k}, i = I - 1, I - 2, \dots, 1, \\ \frac{v_{i}^{k} - v_{i}^{k - 1}}{τ} = d ▵_{h}^{2} v_{i}^{k} + \frac{v_{i}^{K} - d ▵_{h} φ_{i}}{φ_{i}} v_{i}^{k}, i = 1, 2, \dots, I - 1, \\ v_{0}^{k} = v_{l} (t_{k}), v_{I}^{k} = v_{r} (t_{k}), \\ v_{i}^{0} = d ▵_{h} {(u_{0})}_{i} + \frac{v_{i}^{K} - d ▵_{h} φ_{i}}{φ_{i}} {(u_{0})}_{i}, 1 \leq i \leq I - 1 . \end{matrix} \end{matrix}

(17)

In general, implicit methods require more computational efforts, but they are unconditionally stable. Moreover, the fully implicit scheme is first-order accurate in time, while the Crank–Nicolson scheme is second-order accurate both in space and time. However, it is known that in some cases, the Crank–Nicolson method may produce oscillations in the numerical solution [42]. On the other hand, Saulyev’s first- and second-kind formulas are unconditionally stable, can be computed in an explicit manner, but they are first-order accurate. To enhance the order of convergence, we need to construct the Barakat and Clark scheme, which couples these two approximations.

There are many other methods that can be used; for example, the unconditionally stable Dufort–Frankel scheme, which is a three-level, second- or fourth-order accurate in space for appropriate choice of the mesh parameters [43]; explicit and unconditionally stable odd–even hopscotch scheme; leapfrog hopscotch scheme and hybrid schemes, combining the above methods. New results in this direction are obtained in [44,45,46,47].

Now, we state our stability and convergence results. We denote by

∥ w ∥ = max_{0 \leq i \leq I} | w_{i} |

Theorem 2.

We let

{v_{i}^{k}}

be the solution of the backward Euler or Crank–Nicolson scheme; then, there exists a constant

C > 0

, independent of τ and h and

τ_{0} > 0

, such that for all

0 < τ \leq τ_{0}

, we have

\begin{matrix} \begin{matrix} max_{0 \leq k \leq K} ∥ v (t_{k}) - V^{k} ∥ \leq C (τ + h^{2}), \\ max_{0 \leq k \leq K} ∥ V^{k} ∥ \leq C \sum_{k = 0}^{K} τ (| v_{l} (t_{k}) | + | v_{r} (t_{k}) |), \end{matrix} \end{matrix}

(18)

where

V^{k} = {(V_{1}^{k}, V_{2}^{k}, \dots, V_{I - 1}^{k})}^{T}

4.2. Iterative Procedure

To solve Problem (15) for identifying

v^{k}

k = 1, 2, \dots, K

, we use Discretizations (15)–(18) and in order to avoid solving large non-linear systems of algebraic equations, we initiate an iteration procedure. As an initial guess, we set

{(v_{i}^{0})}^{(0)} = 0

i = 0, 1, \dots, I

. Then,

{(v_{i}^{0})}^{(l + 1)}

l = 0, 1, \dots

i = 0, 1, \dots, I

is determined by

{(v_{i}^{0})}^{(l + 1)} = d ▵_{h} {(u_{0})}_{i} + \frac{{(v_{i}^{K})}^{(l)} - d ▵_{h} φ_{i}}{φ_{i}} {(u_{0})}_{i}, i = 1, 2, \dots, I - 1,

where

{(v_{i}^{k})}^{(l)}

k = 1, 2, \dots, K

is the solution of the following problems, denoting

v_{i}^{k} = {(v_{i}^{k})}^{(l + 1)}

for simplicity, depending on the scheme used:

(IMIS-1D) Iteration method based on the implicit backward Euler scheme (15):

$\begin{matrix} \begin{matrix} \frac{v_{i}^{k} - v_{i}^{k - 1}}{τ} = d ▵_{h} u_{i}^{k} + \frac{{(v_{i}^{K})}^{(l)} - d ▵_{h} φ_{i}}{φ_{i}} v_{i}^{k}, i = 1, 2, \dots, I - 1, \\ v_{0}^{k} = v_{l} (t_{k}), v_{I}^{k} = v_{r} (t_{k}), \\ v_{i}^{0} = {(v_{i}^{0})}^{(l)} . \end{matrix} \end{matrix}$

(19)
(IMCNS-1D) Iteration method based on the Crank–Nicolson scheme (16):

$\begin{matrix} \begin{matrix} \frac{v_{i}^{k} - v_{i}^{k - 1}}{τ} = d ▵_{h} \frac{v_{i}^{k} + v_{i}^{k - 1}}{2} + \frac{{(v_{i}^{K})}^{(l)} - d ▵_{h} φ_{i}}{φ_{i}} \frac{v_{i}^{k} + v_{i}^{k - 1}}{2}, i = 1, 2, \dots, I - 1, \\ v_{0}^{k} = v_{l} (t_{k}), v_{I}^{k} = v_{r} (t_{k}), \\ v_{i}^{0} = {(v_{i}^{0})}^{(l)} . \end{matrix} \end{matrix}$

(20)
(IMBCS-1D) Iteration method based on the Barakat and Clark scheme obtained by averaging the solution of the first- and second-kind Saulyev’s scheme (17). We suppose that $v^{K}$ in (17) is known. Then, the first-kind approximation is explicit if the computations are performed from the left boundary to the right. Analogically, in a symmetric way, the same is performed for the second-kind approximation. The discretization is realized in an explicit manner if we compute the solution from the right to the left boundary. Therefore, taking $v^{K}$ to be known from previous iteration, we obtain the following iteration process based on the explicit schemes:

$\begin{matrix} \begin{matrix} {(v_{1})}_{i}^{k} = Q_{i}^{- 1} (\frac{v_{i}^{k - 1}}{τ} + d \frac{v_{i + 1}^{k - 1} - v_{i}^{k - 1} + {(v_{1})}_{i - 1}^{k}}{h^{2}}), i = I - 1, I - 2, \dots, 1, \\ {(v_{2})}_{i}^{k} = Q_{i}^{- 1} (\frac{v_{i}^{k - 1}}{τ} + d \frac{{(v_{2})}_{i + 1}^{k} - v_{i}^{k - 1} + v_{i - 1}^{k - 1}}{h^{2}}), i = 1, 2, \dots, I - 1, \\ Q_{i} = \frac{1}{τ} + \frac{d}{h^{2}} + \frac{{(v_{i}^{K})}^{(l)} - d ▵_{h} φ_{i}}{φ_{i}}, \\ {(v_{1})}_{0}^{k} = v_{l} (t_{k}), {(v_{2})}_{I}^{k} = v_{r} (t_{k}), v_{i}^{0} = {(v_{i}^{0})}^{(l)}, \\ v_{i}^{k} = \frac{{(v_{1})}_{i}^{k} + {(v_{2})}_{i}^{k}}{2}, 1 \leq i \leq I - 1, v_{0}^{k} = v_{l} (t_{k}), v_{I}^{k} = v_{r} (t_{k}) . \end{matrix} \end{matrix}$

(21)

The iteration process continues up to reaching the desired accuracy

ϵ

∥ v^{(l + 1)} - v^{(l)} ∥ \leq ϵ,

∥ v^{(l + 1)} - v^{(l)} ∥ \geq ∥ v^{(l)} - v^{(l - 1)} ∥

Solution

u^{k + 1}

is obtained, integrating (8) at each time layer. Applying trapezoidal rule in the right-hand side, we obtain

u^{k + 1} = u^{k} + \frac{τ}{2} (v^{k + 1} + v^{k}), k = 1, 2, \dots, K .

(22)

Thus, the approximation of solution u is obtained by (22).

5. Iterative Method for Solving a 2D Inverse Problem

In this section, we construct the iteration procedure for solving Problems (11) and (12). First, we discretize domain

\bar{Ω} \times [0, T]

\bar{Ω} = [0, X] \times [0, Y]

. The mesh in time is the same as in Section 4. As before, symmetrically in the x- and y-directions, we consider uniform partition of the corresponding intervals

[0, X]

and

[0, Y]

. We let

h_{1} > 0

and

h_{2} > 0

be small numbers. Then,

x_{i} = i h_{1}

i = 0, 1, \dots, I

and

y_{j} = j h_{2}

j = 0, 1, \dots, J

, where

I = X / h_{1}

J = Y / h_{2}

. The numerical solution and the values of the mesh functions at grid nodes are denoted by

w_{i, j}^{k} = w (t_{k}, x_{i}, y_{j})

We construct the same discretizations of Problems (11) and (12) as for the 1D case.

5.1. Discrete Schemes

The implicit backward Euler scheme: we find ${v_{i, j}^{k}}$ , $k = 1, 2, \dots, K$ , such that

$\begin{matrix} \begin{matrix} \frac{v_{i, j}^{k} - v_{i, j}^{k - 1}}{τ} = d ▵_{h} v_{i, j}^{k} + \frac{v_{i, j}^{K} - d ▵_{h} φ_{i}}{φ_{i, j}} v_{i, j}^{k}, i = 1, 2, \dots, I - 1, j = 1, 2, \dots, J - 1, \\ v_{0, j}^{k} = v_{l} (t_{k}, y_{j}), v_{I, j}^{k} = v_{r} (t_{k}, y_{j}), j = 1, 2, \dots, J, \\ v_{i, 0}^{k} = v_{d} (t_{k}, x_{i}), v_{i, J}^{k} = v_{u} (t_{k}, x_{i}), i = 1, 2, \dots, I, \\ v_{i, j}^{0} = d ▵_{h} {(u_{0})}_{i, j} + \frac{v_{i, j}^{K} - d ▵_{h} φ_{i, j}}{φ_{i, j}} {(u_{0})}_{i, j}, i = 1, 2, \dots, I - 1, j = 1, 2, \dots, J - 1, \end{matrix} \end{matrix}$

(23)

where $▵_{h} w_{i, j}^{k} = (w_{i + 1, j}^{k} - 2 w_{i, j}^{k} + w_{i - 1, j}^{k}) / h_{1}^{2} + (w_{i, j + 1}^{k} - 2 w_{i, j}^{k} + w_{i, j - 1}^{k}) / h_{2}^{2}$ .
The Crank–Nicolson scheme. We find ${v_{i, j}^{k}}$ , $k = 1, 2, \dots, K$ , solving the system

$\frac{v_{i, j}^{k} - v_{i, j}^{k - 1}}{τ} = d ▵_{h} \frac{v_{i, j}^{k} + v_{i, j}^{k - 1}}{2} + \frac{v_{i, j}^{K} - d ▵_{h} φ_{i, j}}{φ_{i, j}} \frac{v_{i, j}^{k} + v_{i, j}^{k - 1}}{2},$

(24)

$i = 1, 2, \dots, I - 1$ , $j = 1, 2, \dots, J - 1$ , with the same initial and boundary discrete conditions as in (23).
The ADE approximation, using the averaging of Saulyev’s first- and second-kind formulae. In discrete scheme (23), we replace $▵_{h} v_{i, j}^{k}$ by $▵_{h}^{1} v_{i, j}^{k}$ and $▵_{h}^{2} v_{i, j}^{k}$ , where

$\begin{matrix} ▵_{h}^{1} v_{i, j}^{k} = \frac{v_{i + 1, j}^{k - 1} - v_{i, j}^{k - 1} - v_{i, j}^{k} + v_{i - 1}^{k}}{h_{1}^{2}} + \frac{v_{i, j + 1}^{k - 1} - v_{i, j}^{k - 1} - v_{i, j}^{k} + v_{i, j - 1}^{k}}{h_{2}^{2}}, \\ ▵_{h}^{2} v_{i, j}^{k} = \frac{v_{i + 1, j}^{k} - v_{i, j}^{k - 1} - v_{i, j}^{k} + v_{i - 1, j}^{k - 1}}{h_{1}^{2}} + \frac{v_{i, j + 1}^{k} - v_{i, j}^{k - 1} - v_{i, j}^{k} + v_{i, j - 1}^{k - 1}}{h_{2}^{2}} . \end{matrix}$

We combine these approximations in order to construct the Barakat and Clark scheme. Namely, we find

{v_{i, j}^{k}}

k = 1, 2, \dots, K

, averaging the solution of the numerical schemes at each time level

\begin{matrix} \begin{matrix} \frac{v_{i, j}^{k} - v_{i, j}^{k - 1}}{τ} = d ▵_{h}^{1} v_{i, j}^{k} + \frac{v_{i, j}^{K} - d ▵_{h} φ_{i}}{φ_{i, j}} v_{i, j}^{k}, i = I - 1, I - 2, \dots, 1, j = J - 1, J - 2, \dots, 1, \\ \frac{v_{i, j}^{k} - v_{i, j}^{k - 1}}{τ} = d ▵_{h}^{2} v_{i, j}^{k} + \frac{v_{i, j}^{K} - d ▵_{h} φ_{i, j}}{φ_{i, j}} v_{i, j}^{k}, i = 1, 2, \dots, I - 1, j = 1, 2, \dots, J - 1, \\ v_{0, j}^{k} = v_{l} (t_{k}, y_{j}), v_{I, j}^{k} = v_{r} (t_{k}, y_{j}), j = 1, 2, \dots, J, \\ v_{i, 0}^{k} = v_{d} (t_{k}, x_{i}), v_{i, J}^{k} = v_{u} (t_{k}, x_{i}), i = 1, 2, \dots, I, \\ v_{i, j}^{0} = d ▵_{h} {(u_{0})}_{i, j} + \frac{v_{i, j}^{K} - d ▵_{h} φ_{i, j}}{φ_{i, j}} {(u_{0})}_{i, j}, 1 \leq i \leq I - 1, 1 \leq j \leq J - 1 . \end{matrix} \end{matrix}

(25)

5.2. Iterative Method

Problems (11) and (12) for finding

v^{k}

k = 1, 2, \dots, K

are solved by the iterative method, based on Discretizations (23)–(25). The initial guess is

{(v_{i, j}^{0})}^{(0)} = 0

i = 0, 1, \dots, I

j = 0, 1, \dots, J

. Then,

{(v_{i, j}^{0})}^{(l + 1)}

l = 0, 1, \dots

i = 0, 1, \dots, I

j = 0, 1, \dots, J

is determined from

{(v_{i, j}^{0})}^{(l + 1)} = d ▵_{h} {(u_{0})}_{i, j} + \frac{{(v_{i, j}^{K})}^{(l)} - d ▵_{h} φ_{i, j}}{φ_{i, j}} {(u_{0})}_{i, j}, i = 1, 2, \dots, I - 1, j = 1, 2, \dots, J - 1,

where

{(v_{i, j}^{k})}^{(l)}

k = 1, 2, \dots, K

is the solution of discrete Problems (23)–(25), where, for simplicity, we denote

v_{i, j}^{k} = {(v_{i, j}^{k})}^{(l + 1)}

(IMIS-2D) The iteration method based on the implicit backward Euler scheme (23):

$\begin{matrix} \begin{matrix} \frac{v_{i, j}^{k} - v_{i, j}^{k - 1}}{τ} = d ▵_{h} u_{i, j}^{k} + \frac{{(v_{i, j}^{K})}^{(l)} - d ▵_{h} φ_{i}}{φ_{i, j}} v_{i, j}^{k}, i = 1, 2, \dots, I - 1, j = 1, 2, \dots, J, \\ v_{0, j}^{k} = v_{l} (t_{k}, y_{j}), v_{I, j}^{k} = v_{r} (t_{k}, y_{j}), j = 1, 2, \dots, J, \\ v_{i, 0}^{k} = v_{d} (t_{k}, x_{i}), v_{i, J}^{k} = v_{u} (t_{k}, x_{i}), i = 1, 2, \dots, I, \\ v_{i, j}^{0} = {(v_{i, j}^{0})}^{(l)} . \end{matrix} \end{matrix}$

(26)
(IMCNS-2D) The iteration method based on the Crank–Nicolson scheme (24):

$\begin{matrix} \begin{matrix} \frac{v_{i, j}^{k} - v_{i, j}^{k - 1}}{τ} = d ▵_{h} \frac{v_{i, j}^{k} + v_{i, j}^{k - 1}}{2} \\ + \frac{{(v_{i, j}^{K})}^{(l)} - d ▵_{h} φ_{i, j}}{φ_{i, j}} \frac{v_{i, j}^{k} + v_{i, j}^{k - 1}}{2}, i = 1, 2, \dots, I - 1, j = 1, 2, \dots, J, \\ v_{0, j}^{k} = v_{l} (t_{k}, y_{j}), v_{I, j}^{k} = v_{r} (t_{k}, y_{j}), j = 1, 2, \dots, J, \\ v_{i, 0}^{k} = v_{d} (t_{k}, x_{i}), v_{i, J}^{k} = v_{u} (t_{k}, x_{i}), i = 1, 2, \dots, I, \\ v_{i, j}^{0} = {(v_{i, j}^{0})}^{(l)} . \end{matrix} \end{matrix}$

(27)
(IMBCS-2D) The iteration method obtained by averaging the solution of the first- and the second-kind Saulyev’s scheme (25). Taking $v^{K}$ to be known from previous iteration, the first approximation is explicit if the computations are performed passing the nodes from the bottom boundary to the top and from the left boundary to the right. Similarly, the second-kind approximation is performed in an explicit manner if we compute the solution from the upper to the bottom boundary and from the right to the left boundary. Thus, we obtain the following iteration process based on the explicit schemes:

$\begin{matrix} \begin{matrix} {(v_{1})}_{i, j}^{k} = Q_{i, j}^{- 1} (\frac{v_{i, j}^{k - 1}}{τ} + d \frac{v_{i + 1, j}^{k - 1} - v_{i, j}^{k - 1} + {(v_{1})}_{i - 1, j}^{k}}{h_{1}^{2}} + d \frac{v_{i, j + 1}^{k - 1} - v_{i, j}^{k - 1} + {(v_{1})}_{i, j - 1}^{k}}{h_{1}^{2}}), \\ i = I - 1, I - 2, \dots, 1, j = J - 1, J - 2, \dots, 1, \\ {(v_{2})}_{i, j}^{k} = Q_{i, j}^{- 1} (\frac{v_{i, j}^{k - 1}}{τ} + d \frac{{(v_{2})}_{i + 1, j}^{k} - v_{i, j}^{k - 1} + v_{i - 1, j}^{k - 1}}{h_{1}^{2}} + d \frac{{(v_{2})}_{i, j + 1}^{k} - v_{i, j}^{k - 1} + v_{i, j - 1}^{k - 1}}{h_{2}^{2}}), \\ i = 1, 2, \dots, I - 1, j = 1, 2, \dots, J - 1, \\ Q_{i, j} = \frac{1}{τ} + \frac{2 d}{h^{2}} + \frac{{(v_{i, j}^{K})}^{(l)} - d ▵_{h} φ_{i, j}}{φ_{i}}, \\ {(v_{1})}_{0, j}^{k} = v_{l} (t_{k}, y_{j}), {(v_{2})}_{I, j}^{k} = v_{r} (t_{k}, y_{j}), j = 1, 2, \dots, J, \\ {(v_{1})}_{i, 0}^{k} = v_{d} (t_{k}, x_{i}), {(v_{2})}_{i, J}^{k} = v_{u} (t_{k}, x_{i}), i = 1, 2, \dots, I, \\ v_{i, j}^{0} = {(v_{i, j}^{0})}^{(l)}, \\ v_{i, j}^{k} = \frac{{(v_{1})}_{i, j}^{k} + {(v_{2})}_{i, j}^{k}}{2}, 1 \leq i \leq I - 1, 1 \leq j \leq J - 1, \\ v_{0, j}^{k} = v_{l} (t_{k}, y_{j}), v_{I, j}^{k} = v_{r} (t_{k}, y_{j}), j = 1, 2, \dots, J, \\ v_{i, 0}^{k} = v_{d} (t_{k}, x_{i}), v_{i, J}^{k} = v_{u} (t_{k}, x_{i}), i = 1, 2, \dots, I . \end{matrix} \end{matrix}$

(28)

The iteration process continues up to reaching the desired accuracy

ϵ

∥ v^{(l + 1)} - v^{(l)} ∥ \leq ϵ, where ∥ w_{i, j} ∥ = max_{0 \leq i \leq I} max_{0 \leq j \leq J} | w_{i, j} |,

(29)

∥ v^{(l + 1)} - v^{(l)} ∥ \geq ∥ v^{(l)} - v^{(l - 1)} ∥

Solution u is obtained from (22).

6. Numerical Tests

In this section, we illustrate the accuracy and efficiency of the proposed iterative procedure for solving 1D and 2D inverse Problems (11)–(13).

Example 1

(1D problem: convergence test). In this test example, we verify the order of convergence of iteration methods IMIS-1D, IMCNS-1D, IMBCS-1D for solving 1D Problem (13) for given functions

u_{0} (x)

and

φ (x)

. The accuracy of these approaches is computed for modified Problem (13),

T = 1

X = 1

with exact solution

v (t, x) = - e^{- t} sin (π x)

. To this aim, we add residual functions in the right-hand side of the equation in (13) and set

φ (x) = e^{- T} sin (π x)

u_{0} (x) = sin (π x)

We offer the absolute errors and orders of convergence in maximal and

L_{2}

norms at a final time, as well as the relative errors in the maximal discrete norm and the corresponding orders of convergence

\begin{matrix} E_{\infty} = E_{\infty}^{I} = ∥ v^{K} - v (x, T) ∥, C R_{\infty} = {log}_{2} \frac{E^{I}}{E^{2 I}}, \\ E_{2} = E_{2}^{I} = ∥ v^{K} {- v (x, T) ∥}_{2}, C R_{2} = {log}_{2} \frac{E_{2}^{I}}{E_{2}^{2 I}}, {∥ w ∥}_{2} = \sqrt{\sum_{i = 0}^{I} h w_{i}^{2}}, \\ E_{R} = E_{R}^{I} = ∥\frac{v^{K} - v (x, T)}{v (x, T)}∥, C R_{R} = {log}_{2} \frac{E_{R}^{I}}{E_{R}^{2 I}} . \end{matrix}

Since the expected order of convergence of IMIS-1D and IMBCS-1D is

O (τ + h^{2})

, the computations are performed for

τ \approx h^{2}

. Similarly, since the expected order of convergence of IMCNS-1D is

O (τ^{2} + h^{2})

, for all runs, we take

τ \approx h

. The accuracy is

ϵ = 1 \times 10^{- 7}

. The results from computations are given in Table 1, Table 2 and Table 3. The number of iterations is denoted by

i t e r

. We observe that only a few iterations are needed to reach the desired accuracy, and because of the fixed ratio between mesh step sizes, the order of convergence is as expected.

Example 2

(2D problem: convergence test). Now, we verify the order of convergence of the iteration methods for solving the 2D inverse problem for identifying

u (t, x, y)

and

μ (x, y)

in Equation (3) associated with zero Dirichlet boundary conditions and given observations (7). To this aim, we consider the corresponding modified problem, setting

μ (x, y) = 2 x^{2} + y + 1

and adding residual function

f (t, x, y)

in the right-hand side of Equation (3), such that the exact solution is

u = e^{- t} sin (π x) sin (π y)

and

φ (x, y) = e^{- T} sin (π x) sin (π y)

We solve the inverse problem by iterative methods IMIS-2D, IMCNS-2D, IMBCS-2D for

I = J

, i.e.,

h_{1} = h_{2} = h

. As in Example 1, the mesh step sizes are fixed according to the expected order of convergence of the corresponding iteration method, namely

τ \approx h^{2}

for IMIS-2D and IMBCS-2D, and

τ \approx h

for IMCNS-2D. For the stopping criteria, we take

ϵ = 1 \times 10^{- 7}

Errors and order of convergence of the restored solution u and

μ

, in maximal and

L_{2}

norms at the final time are defined as follows:

\begin{matrix} E_{\infty} = E_{\infty}^{I} = ∥ u^{K} - u (x, y, T) ∥, C R_{\infty} = {log}_{2} \frac{E^{I}}{E^{2 I}}, \\ E_{2} = E_{2}^{I} = ∥ u^{K} {- u (x, y, T) ∥}_{2}, C R_{2} = {log}_{2} \frac{E_{2}^{I}}{E_{2}^{2 I}}, {∥ w ∥}_{2} = \sqrt{\sum_{i = 0}^{I} \sum_{j = 0}^{J} h_{1} h_{2} w_{i, j}^{2}}, \\ E_{\infty}^{μ} = {(E^{μ})}_{\infty}^{I} = ∥ μ - μ (x, y) ∥, C R_{\infty} = {log}_{2} \frac{{(E^{μ})}^{I}}{{(E^{μ})}^{2 I}}, \\ E_{2}^{μ} = {(E^{μ})}_{2}^{I} = {∥ μ - μ (x, y) ∥}_{2}, C R_{2} = {log}_{2} \frac{{(E^{μ})}_{2}^{I}}{{(E^{μ})}_{2}^{2 I}} . \end{matrix}

Since the modified problem for Equation (3) has a right-hand side

f (t, x, y)

, we also modify Equation (11) in the corresponding to (11) and (12) problems as follows:

\begin{matrix} \begin{matrix} \frac{\partial v}{\partial t} = d ▵ v + \frac{v (T, x, y) - d ▵ φ (x, y) - f (t, x, y)}{φ (x, y)} v + \frac{\partial f (t, x, y)}{\partial t}, \\ v (0, x, y) = d ▵ u_{0} (x, y) + \frac{v (T, x, y) - d ▵ φ (x, y) - f (t, x, y)}{φ (x, y)} u_{0} (x, y), \end{matrix} \end{matrix}

(30)

and

μ (x, t) = [- v (T, x, y) + d ▵ φ (x, y) + f (t, x, y)] / φ (x, y)

. Similarly, iterative methods IMIS-2D, IMCNS-2D, IMBCS-2D are modified according to (30).

Computational results presented in Table 4, Table 5 and Table 6 confirm that the iteration process for recovering both reaction coefficient

μ

and solution u has the same order of convergence as the underlying numerical scheme.

Example 3

(2D problem: noisy data). The test problem is as in Example 2, but with perturbed measurements (7) defined as in [48]:

φ_{ρ} (x, y) = φ (x, y) + ρ_{0} (ϱ (x, y) - 0.5),

(31)

where

ϱ (x, y)

is a random function, uniformly distributed in domain [0, 2] × [0, 2] and

ρ_{0}

is the level of noise.

We compare the precision of the three methods for solving the inverse problem for one and the same temporal mesh. In Table 7, we present computational results for different noise levels and time step sizes. For all runs, we set

N_{1} = 40

and

ϵ = 1 \times 10^{- 5}

. The results show that the better precision is obtained b IMIS-2D. Moreover, the bigger the noise level, the finer the time mesh should be.

Example 4

(2D problem: real data). For the simulations, we take the data from [13], in which the authors collected the results from investigations and observations reported in many other papers. Just as in [13], we assume that

-: Approximately 25% of the bees in the colony are forager;
-: The colony’s bee population ranges from 20,000 to 60,000 individuals;
-: The foraging domain with rectangular symmetry measures 2 km by 2 km and is represented as [0, 2] × [0, 2];
-: The beehive is positioned in the center (1.0 km, 1.0 km) of the terrain;
-: The diffusion rate is 0.1 km²/day;
-: Forager bees tend to stay away from the boundaries of this area;
-: The initial state at the beehive is $u_{0} (x, y) = N_{0} G (x, y)$ , where $N_{0}$ = 10,000 and $G (x, y)$ represents a symmetric Gaussian density function (1), centered at $(1, 1)$ with standard deviation 1.5812 × 10⁻⁴;
-: Mortality rate $μ = 1 / 6.5 {day}^{- 1}$ , maximal production parameter $β =$ 2421.13, sigmoidal Hill production parameter $χ =$ 7173.56, stable and unstable uncontaminated equilibrium population size ${\bar{U}}_{+} = 10, 000$ , ${\bar{U}}_{-} =$ 7000, respectively;
-: In the case $q = 0$ , $μ = c o n s t$ , at the sunrise of the days $n =$ 1,2,…, the density and the total number of the forager bees are represented by [13]

$\begin{matrix} u (n + 1, x, y) = G (x, y) U (n + 1), \end{matrix}$

(32)

$\begin{matrix} U (n + 1) = \frac{β U {(n)}^{2}}{{\tilde{χ}}^{2} + U {(n)}^{2}} + e^{- μ} U (n), n = 0, 1, \dots ., \\ {\tilde{χ}}^{2} = \frac{{({\bar{U}}_{+} + {\bar{U}}_{-})}^{2} - {({\bar{U}}_{+} - {\bar{U}}_{-})}^{2}}{4} = χ^{2} e^{2 μ} . \end{matrix}$

(33)

For the simulations, we take

I = J = 100

τ = 0.005

and

ϵ = 1 . e^{- 3}

. The computations are performed day by day, starting with initial Condition (32), namely

u^{0} : = u (n + 1, x, y)

. In order to avoid an extremely large initial condition for the inverse problem, for this practical example, we obtain

v^{0}

, integrating (8) from

t_{0} : = t_{n}

to T, depending on which day the simulations are performed. Thus, applying the exact and rectangle method integration, we derive

v^{0} = (u^{K} - u^{0}) / (T - t_{0})

. Moreover, we rescale both the inverse and direct problems, setting

u = p u

p = 2 π σ^{2}

. As a result, we obtain the same direct and inverse problems (2), (3), (11) and (12), with initial condition

u_{0} : = p u_{0}

. To generate measurements

φ

, we solve the direct problem.

We compute the solution of the inverse problem by IMIS-2D in the first two days, since in that period the bees may remain uncontaminated [13], even in regions with pesticide.

In Figure 1, we depict the recovered mortality rate

μ = max (0, μ)

in the first day after 6 h and on the second day after 8 h. Although the recovered rate is not a constant, we observe that the values of

μ

are close to the exact one, namely

μ = 0.1538 / 1 d a y^{- 1}

The center-symmetric pictures in Figure 2, Figure 3 and Figure 4 show the spatial density of uncontaminated bees obtained by the inverse problem. Over the first two days, uncontaminated bees are progressively moving away from the hive. This behavior and location pattern matches the one presented in [13]. The regions with different density are smoother, since we study only the contamination-free phenomenon, and as a result we consider a simpler model in comparison with [13].

In Figure 5, the decreasing value of the norm (29) at each iteration is illustrated.

7. Conclusions

In this paper, we presented a numerical analysis of forager bee looses in a spatial environment without contamination. The approach was based on the inverse problem for recovering the mortality rate in a honeybee difference equation model that describes forages leaving and entering the hive each day. The mathematical formulation of this process lead to the study of the reaction coefficient identification on the base of final time measurements of the forage bee concentration. After introducing a new solution variable, we reduced the inverse problem to a forward one with a non-local difference operator and a non-local initial condition. For the numerical solution of the problem, different time-stepping techniques were applied to the ODEs system, resulting from second-order difference approximation.

We illustrated, by numerical tests, that for exact measurements, the order of convergence of the solution of the inverse problem for all three methods is the same as in the underlying methods. In the case of perturbed measurements, the most efficient was the iterative approach based on implicit time-stepping. We also observed that in order to achieve good accuracy of the recovered solution, the larger the deviation, the smaller the time step size should be.

The simulations with real data illustrate that the solution of the inverse problem successfully recovers the bee density and the mortality rate of the forage bees and produces relevant results.

In our forthcoming work, we plan to propose a numerical method for solving the inverse problem for a spatial model describing the collapse of honey bee colonies caused by the contamination of foraging bees with pesticides.

Author Contributions

Conceptualization, L.G.V.; methodology, M.N.K. and L.G.V.; investigation, M.N.K. and L.G.V.; resources, A.Z.A., M.N.K. and L.G.V.; writing—original draft preparation, M.N.K. and L.G.V.; writing—review and editing, A.Z.A., M.N.K. and L.G.V.; validation, M.N.K. All authors have read and agreed to the published version of the manuscript.

Funding

The first author is supported by the Bulgarian National Science Fund under the Project KP-06-PN 46-7 “Design and research of fundamental technologies and methods for precision apiculture”. The second and the third authors are supported by the Bulgarian National Science Fund under the Project KP-06-N 62/3 “Numerical methods for inverse problems in evolutionary differential equations with applications to mathematical finance, heat mass transfer, honeybee population and environmental pollution”, 2022.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Acknowledgments

The authors are very grateful to the anonymous reviewers whose valuable comments and suggestions improved the quality of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

Bagheri, S.; Mirzaie, M. A mathematical model of honey bee colony dynamics to predict the effect of pollen on colony failure. PLoS ONE 2019, 14, e0225632. [Google Scholar] [CrossRef]
Chen, J.; DeGrandi-Hoffman, G.; Ratti, V.; Kang, Y. Review on mathematical modeling of honeybee population dynamics. J. Math. Biosci. Eng. 2021, 18, 9606–9650. [Google Scholar] [CrossRef]
Van Dooremalen, C.; Gerritsen, L.; Cornelissen, B.; van der Steen, J.J.; van Langevelde, F.; Blacquiere, T. Winter survival of individual honey bees and honey bee colonies depends on level of varroa destructor infestation. PLoS ONE 2012, 7, e36285. [Google Scholar] [CrossRef]
Fisher, A., II; DeGrandi-Hoffman, G.; Smith, B.H.; Johnson, M.; Kaftanoglu, O.; Cogley, T.; Fewell, J.H.; Harrison, J.F. Colony field test reveals dramatically higher toxicity of a widely-used mito-toxic fungicide on honey bees (apis mellifera). Environ. Pollut. 2021, 269, 115964. [Google Scholar] [CrossRef]
Genersch, E. American foulbrood in honeybees and its causative agent, paenibacillus larvae. J. Invertebr. Pathol. 2010, 103, S10–S19. [Google Scholar] [CrossRef]
Laomettachit, T.; Liangruksa, M.; Termsaithong, T.; Tangthanawatsakul, A.; Duangphakdee, O. A model of infection in honeybee colonies with social immunity. PLoS ONE 2021, 16, e0247294. [Google Scholar] [CrossRef]
Mayack, C.; Naug, D. Energetic stress in the honeybee apis mellifera from nosema ceranae infection. J. Invertebr. Pathol. 2009, 100, 185–188. [Google Scholar] [CrossRef]
Paris, L.; El Alaoui, H.; Delbac, F.; Diogon, M. Effects of the gut parasite nosema ceranae on honey bee physiology and behavior. Curr. Opin. Insect Sci. 2018, 26, 149–154. [Google Scholar] [CrossRef]
Smith, K.M.; Loh, E.H.; Rostal, M.K.; Zambrana-Torrelio, C.M.; Mendiola, L.; Daszak, P. Pathogens, pests, and economics: Drivers of honey bee colony declines and losses. Ecohealth 2013, 10, 434–445. [Google Scholar] [CrossRef]
Williams, B.A. Unique physiology of host-parasite interactions in microsporidia infections. Cell. Microbiol. 2009, 11, 1551–1560. [Google Scholar] [CrossRef]
DeGrandi-Hoffman, G.; Roth, S.A.; Loper, G.; Erickson, E., Jr. Beepop: A honeybee population dynamics simulation model. Ecol. Modell. 1989, 45, 133–150. [Google Scholar] [CrossRef]
Becher, M.A.; Osborne, J.L.; Thorbek, P.; Kennedy, P.J.; Grimm, V. Towards a systems approach for understanding honeybee decline: A stocktaking and synthesis of existing models. J. Appl. Ecol. 2013, 50, 868–880. [Google Scholar] [CrossRef]
Magal, P.; Webb, G.F. A spatial model of honeybee colony collapse due to pesticide contamination of foraging bees. J. Math. Biol. 2020, 80, 2363–2393. [Google Scholar] [CrossRef]
Bastaansen, R.; Doelman, A.; van Langevede, F.; Rottschafer, V. Modeling honey bee colonies in winter using a Keller-Segel model with a sign-changing chemotectic coefficient. SIAM J. Appl. Math. 2020, 80, 839–863. [Google Scholar] [CrossRef]
Watmough, J.; Camazine, S. Self-organized thermoregulation of honeybee clusters. J. Theor. Biol. 1995, 176, 391–402. [Google Scholar] [CrossRef]
Kribs-Zaleta, C.M.; Mitchell, C. Modeling colony collapse disorder in honeybees as a contagion. Math. Biosci. Engn. 2014, 11, 1275–1294. [Google Scholar] [CrossRef]
Romero-Leiton, J.P.; Gutierrez, A.; Benavides, I.F.; Molina, O.E.; Pulgarin, A. An approach to the modeling of honey bee colonies. Web Ecol. 2022, 22, 7–19. [Google Scholar] [CrossRef]
Atanasov, A.Z.; Koleva, M.N.; Vulkov, L.G. Numerical analysis of thermoregulation in honey bee colonies in winter based on sign-changing chemotactic coefficient model. In International Conference on New Trends in the Applications of Differential Equations in Sciences; Springer Proceedings in Mathematics & Statistics; Springer: Cham, Switzerland, 2023; Volume 4122, pp. 69–279. [Google Scholar]
Chavent, G. Nonlinear Least Squares Problems: Theoretical Foundation and Step-by Guide for Applications; Springer: Cham, Switzerland, 2009. [Google Scholar]
Hasanov, A.H.; Romanov, V.G. Introduction to Inverse Problems for Differential Equations, 1st ed.; Springer: Cham, Switzerland, 2017; 261p. [Google Scholar]
Isakov, V. Inverse Problems for Partial Differential Equations, 3rd ed.; Springer: Cham, Switzerland, 2017; p. 406. [Google Scholar]
Ivanov, V.K.; Vasin, V.V.; Tanana, V.P. Theory of Linear Ill-Posed Problems and Its Approximations; Walter de Gruyter: Nauka, Moscow, 1978. (In Russian) [Google Scholar]
Kabanikhin, S.I. Inverse and Ill-Posed Problems; DeGruyer: Berlin, Germany, 2011. [Google Scholar]
Lesnic, D. Inverse Problems with Applications in Science and Engineering; CRC Press: Abingdon, UK, 2021; p. 349. [Google Scholar]
Prilepko, A.I.; Orlovsky, D.G.; Vasin, I.A. Methods for Solving Inverse Problems in Mathematical Physics; Marcel Dekker: New York, NY, USA, 2000. [Google Scholar]
Prilepko, A.I.; Kostin, A.B.; Solov’ev, V.V. Inverse source and coefficient problems for elliptic and parabolic equations in Hölder and Sobolev spaces. J. Math. Sci. 2019, 237, 576–594. [Google Scholar] [CrossRef]
Romanov, V.G.; Hasanov, A.H. Uniquiness and stability analysis of final data inverse source problems for evolution equations. J. Inverse Ill-Posed Probl. 2022, 30, 425–446. [Google Scholar] [CrossRef]
Cao, K.; Lesnic, D. Reconstruction of the perfusion coefficient from temperature measurements using the conjugate gradient method. Int. J. Comput. Math. 2018, 95, 797–814. [Google Scholar] [CrossRef]
Atanasov, A.; Georgiev, S.; Vulkov, L. Reconstruction analysis of honeybee colony collapse disorder modeling. Optim. Eng. 2021, 22, 2481–2503. [Google Scholar] [CrossRef]
Atanasov, A.Z.; Georgiev, S.G.; Vulkov, L.G. Parameter Estimation Analysis in a Model of Honey Production. Axioms 2023, 12, 214. [Google Scholar] [CrossRef]
Atanasov, A.Z.; Koleva, M.N.; Vulkov, L.G. Numerical optimization identification of a Keller-Segel model for thermoregulation in honey bee colonies in winter. Model. Dev. Intell. Syst. Commun. Comput. Inf. Sci. 2023, 1761, 279–293. [Google Scholar]
Atanasov, A.Z.; Koleva, M.N.; Vulkov, L.G. Parameter estimation inspired by temperature measurements for a chemotactic model of honeybee thermoregulation. In International Conference on Numerical Methods and Applications; Lecture Notes in Computer Science; Springer: Cham, Switzerland, 2023; Volume 13858, pp. 36–47. [Google Scholar]
Koleva, M.; Vulkov, L. Reconstruction coefficient analysis of honeybee collapse due to pesticide contamination. J. Phys. Conf. Proc. 2023. accepted. [Google Scholar]
Martin-Vaquero, J.; Sajvicius, S. The two level finite difference scheme for the heat equation with nonlocal initial condition. Appl. Math. Comput. 2019, 342, 160–177. [Google Scholar] [CrossRef]
Cong, C.; Li, D.; Li, L.; Zhao, D. Crank-Nicolson compact difference scheme for a class of nonlocal nonlinear parabolic problems. Comput. Math. Appl. 2023, 132, 1–17. [Google Scholar]
Dehghan, M. Numerical schemes for one-dimensional parabolic equations with nonstandard initial condition. Appl. Math. Comput. 2004, 147, 321–331. [Google Scholar] [CrossRef]
Dehghan, M. Implicit collocation technique for heat equation with non-classic initial condition. Int. J. Nonlinear Sci. Numer. Simul. 2006, 7, 461–466. [Google Scholar] [CrossRef]
Lin, Y. Analytical and numerical solutions for a class nonlocal nonlinear parabolic differential equations. SIAM J. Math. Anal. 1994, 25, 1577–1594. [Google Scholar] [CrossRef]
Lin, Y. Finite-difference solutions for parabolic equations with time weighting initial conditions. Appl. Math. Comput. 1994, 65, 49–61. [Google Scholar] [CrossRef]
Magal, P.; Webb, G.F.; Wu, Y. Environmental model of honey bee colony collapse due to pesticide contamination. Bull. Math. Biol. 2019, 81, 4908–4931. [Google Scholar] [CrossRef] [PubMed]
Ratti, V.; Kevan, P.G.; Eberl, H.J. A mathematical model of forager loss in honeybee colonies infested with Varroa destructor and accute bee paralysis virus. Bull. Math. Biol. 2017, 79, 1218–1253. [Google Scholar] [CrossRef] [PubMed]
Il’in, A.M. Differencing Scheme for a differential equation with a small parameter affecting the highest derivative. Mat. Zametki 1969, 6, 237–248. [Google Scholar] [CrossRef]
Dehghan, M. Three-level techniques for one-dimensional parabolic equation with nonlinear initial condition. Appl. Math. Comput. 2004, 151, 567–579. [Google Scholar] [CrossRef]
Askar, A.H.; Omle, I.; Kovács, E.; Majár, J. Testing some different implementations of heat convection and radiation in the Leapfrog-Hopscotch algorithm. Algorithms 2022, 15, 400. [Google Scholar] [CrossRef]
Jalghaf, H.K.; Kovács, E.; Majár, J.; Nagy, Á.; Askar, A.H. Explicit stable finite difference methods for diffusion-reaction type equations. Mathematics 2021, 9, 3308. [Google Scholar] [CrossRef]
Nagy, Á.; Omle, I.; Kareem, H.; Kovacs, E.; Barna, I.F.; Bognar, G. Stable, Explicit, Leapfrog-Hopscotch Algorithms for the Diffusion Equation. Computation 2021, 9, 92. [Google Scholar] [CrossRef]
Saleh, M.; Kovács, E.; Barna, I.F.; Mátyás, L. New analytical results and comparison of 14 numerical schemes for the diffusion equation with space-dependent diffusion coefficient. Mathematics 2022, 10, 2813. [Google Scholar] [CrossRef]
Samarskii, A.A.; Vabishchevich, P.N. Numerical Methods for Solving Inverse Problems in Mathematical Physics; de Gruyter: Berlin, Germany, 2007; 438p. [Google Scholar]

Figure 1. Recovered

μ

at Day 1, after 6 h (left), and at Day 2 after 8 h (right), Example 4.

Figure 1. Recovered

μ

at Day 1, after 6 h (left), and at Day 2 after 8 h (right), Example 4.

Figure 2. Recovered bee density at Day 1 after 2 h (left), Day 2 after 2 h (right), Example 4.

Figure 3. Recovered bee density at Day 1 after 6 h (left), Day 2 after 6 h (right), Example 4.

Figure 4. Recovered bee density at Day 1 after 8 h (left), Day 2 after 8 h (right), Example 4.

Figure 5. Values of the norm (29) at each iteration at Day 1 after 6 h (left), Day 2 after 6 h (right), Example 4.

Table 1. Errors and convergence rate of the solution of IMIS-1D,

τ = h^{2}

, Example 1.

Table 1. Errors and convergence rate of the solution of IMIS-1D,

τ = h^{2}

, Example 1.

I	$τ$	$E_{\infty}$	${CR}_{\infty}$	$E_{2}$	${CR}_{2}$	$E_{R}$	${CR}_{R}$	$iter$
20	2.500000 × 10⁻³	6.8620 × 10⁻⁴		4.8594 × 10⁻⁴		1.9222 × 10⁻³		4
40	6.250000 × 10⁻⁴	1.7141 × 10⁻⁴	2.0012	1.2139 × 10⁻⁴	2.0012	4.8026 × 10⁻⁴	2.0009	4
80	1.562500 × 10⁻⁴	4.2850 × 10⁻⁵	2.0001	3.0340 × 10⁻⁵	2.0002	1.2006 × 10⁻⁴	2.0002	4
160	3.906250 × 10⁻⁵	1.0713 × 10⁻⁵	1.9999	7.5846 × 10⁻⁶	2.0001	3.0010 × 10⁻⁵	2.0001	4
320	9.765625 × 10⁻⁶	2.6782 × 10⁻⁶	2.0001	1.8961 × 10⁻⁶	2.0001	7.5030 × 10⁻⁶	2.0000	4

Table 2. Errors and convergence rate of the solution of IMCNS-1D,

τ = h

, Example 1.

Table 2. Errors and convergence rate of the solution of IMCNS-1D,

τ = h

, Example 1.

I	$E_{\infty}$	${CR}_{\infty}$	$E_{2}$	${CR}_{2}$	$E_{R}$	${CR}_{R}$	$iter$
20	5.4756 × 10⁻⁴		3.8788 × 10⁻⁴		1.5432 × 10⁻³		6
40	1.3681 × 10⁻⁴	2.0009	9.6909 × 10⁻⁵	2.0009	3.8567 × 10⁻⁴	2.0005	5
80	3.4207 × 10⁻⁵	1.9997	2.4223 × 10⁻⁵	2.0002	9.6408 × 10⁻⁵	2.0001	5
160	8.5516 × 10⁻⁶	2.0000	6.0555 × 10⁻⁶	2.0001	2.4100 × 10⁻⁵	2.0001	4
320	2.1378 × 10⁻⁶	2.0000	1.5138 × 10⁻⁶	2.0001	6.0229 × 10⁻⁶	2.0005	4

Table 3. Errors and convergence rate of the solution of IMBCS-1D,

τ = h^{2}

, Example 1.

Table 3. Errors and convergence rate of the solution of IMBCS-1D,

τ = h^{2}

, Example 1.

I	$E_{\infty}$	${CR}_{\infty}$	$E_{2}$	${CR}_{2}$	$E_{R}$	${CR}_{R}$	$iter$
20	1.5768 × 10⁻³		1.0677 × 10⁻³		4.2919 × 10⁻³		4
40	4.3687 × 10⁻⁴	1.8517	3.0230 × 10⁻⁴	1.8204	1.1912 × 10⁻³	1.8493	4
80	1.1413 × 10⁻⁴	1.9366	7.98838 × 10⁻⁵	1.9200	3.1212 × 10⁻⁴	1.9322	4
160	2.9114 × 10⁻⁵	1.9708	2.0494 × 10⁻⁵	1.9627	7.9964 × 10⁻⁵	1.9647	4
320	7.3485 × 10⁻⁶	1.9862	5.1876 × 10⁻⁶	1.9821	2.0279 × 10⁻⁵	1.9794	4

Table 4. Errors and convergence rate of the solution of IMIS-2D,

τ = h^{2}

, Example 2.

Table 4. Errors and convergence rate of the solution of IMIS-2D,

τ = h^{2}

, Example 2.

I	$E_{\infty}$	${CR}_{\infty}$	$E_{2}$	${CR}_{2}$	$E_{\infty}^{μ}$	${CR}_{\infty}^{μ}$	$E_{2}^{μ}$	${CR}_{2}^{μ}$	$iter$
20	5.9575 × 10⁻⁴		2.9777 × 10⁻⁴		9.9173 × 10⁻⁴		9.3245 × 10⁻⁴		5
40	1.4881 × 10⁻⁴	2.0012	7.4380 × 10⁻⁵	2.0012	2.4773 × 10⁻⁴	2.0011	2.3904 × 10⁻⁴	1.9638	5
80	3.7195 × 10⁻⁵	2.0003	1.8591 × 10⁻⁵	2.0003	6.1904 × 10⁻⁵	2.0007	6.0498 × 10⁻⁵	1.9823	5
160	9.2975 × 10⁻⁶	2.0002	4.6404 × 10⁻⁶	2.0002	1.5472 × 10⁻⁵	2.0004	1.5192 × 10⁻⁵	1.9936	5
320	2.3242 × 10⁻⁶	2.0001	1.1601 × 10⁻⁶	2.0000	3.8674 × 10⁻⁶	2.0002	3.8038 × 10⁻⁶	1.9978	5

Table 5. Errors and convergence rate of the solution of IMCNS-2D,

τ = h

, Example 2.

Table 5. Errors and convergence rate of the solution of IMCNS-2D,

τ = h

, Example 2.

I	$E_{\infty}$	${CR}_{\infty}$	$E_{2}$	${CR}_{2}$	$E_{\infty}^{μ}$	${CR}_{\infty}^{μ}$	$E_{2}^{μ}$	${CR}_{2}^{μ}$	$iter$
20	5.6981 × 10⁻⁴		2.8480 × 10⁻⁴		9.4954 × 10⁻⁴		8.9176 × 10⁻⁴		5
40	1.4234 × 10⁻⁴	2.0012	7.1144 × 10⁻⁵	2.0012	2.3712 × 10⁻⁴	2.0016	2.2863 × 10⁻⁴	1.9636	5
80	3.5577 × 10⁻⁵	2.0003	1.7782 × 10⁻⁵	2.0003	5.9261 × 10⁻⁵	2.0005	5.7880 × 10⁻⁵	1.9819	5
160	8.8939 × 10⁻⁶	2.0001	4.4454 × 10⁻⁶	2.0001	1.4814 × 10⁻⁵	2.0002	1.4561 × 10⁻⁵	1.9910	5
320	2.2236 × 10⁻⁶	2.0000	1.1114 × 10⁻⁶	2.0000	3.7034 × 10⁻⁶	2.0000	3.6511 × 10⁻⁶	1.9957	5

Table 6. Errors and convergence rate of the solution of IMBCS-2D,

τ = h^{2}

, Example 2.

Table 6. Errors and convergence rate of the solution of IMBCS-2D,

τ = h^{2}

, Example 2.

I	$E_{\infty}$	${CR}_{\infty}$	$E_{2}$	${CR}_{2}$	$E_{\infty}^{μ}$	${CR}_{\infty}^{μ}$	$E_{2}^{μ}$	${CR}_{2}^{μ}$	$iter$
20	1.1899 × 10⁻³		5.7278 × 10⁻⁴		2.6630 × 10⁻³		1.6889 × 10⁻³		5
40	3.1606 × 10⁻⁴	1.9126	1.5524 × 10⁻⁴	1.8835	8.9879 × 10⁻⁴	1.5670	4.8498 × 10⁻⁴	1.8001	5
80	8.0990 × 10⁻⁵	1.9644	4.0237 × 10⁻⁵	1.9479	2.8785 × 10⁻⁴	1.6426	1.3132 × 10⁻⁴	1.8849	5
160	2.0468 × 10⁻⁵	1.9844	1.0228 × 10⁻⁵	1.9760	7.9171 × 10⁻⁵	1.8623	3.4373 × 10⁻⁵	1.9337	5
320	5.1263 × 10⁻⁶	1.9974	2.5777 × 10⁻⁶	1.9884	2.0558 × 10⁻⁵	1.9453	8.6548 × 10⁻⁶	1.9897	5

Table 7. Errors and convergence rate of the solution of inverse methods,

N_{1} = 40

, Example 3.

Table 7. Errors and convergence rate of the solution of inverse methods,

N_{1} = 40

, Example 3.

$ρ_{0}$	$τ$	Method	$E_{\infty}$	$E_{2}$	$iter$	Reach Accuracy
						$ϵ = 1 \times 10^{- 5}$
0.002	0.01	IMIS-2D	6.3931 × 10⁻⁴	2.9009 × 10⁻⁴	4	yes
		IMCNS-2D	fails	fails		no
		IMBCS-2D	3.0327 × 10⁻¹	1.5163 × 10⁻¹	2	no
0.002	0.005	IMIS-2D	6.3919 × 10⁻⁴	2.9014 × 10⁻⁴	4	yes
		IMCNS-2D	2.6857 × 10⁻²	9.7986 × 10⁻³	4	yes
		IMBCS-2D	3.0327 × 10⁻¹	1.5473 × 10⁻¹	2	no
0.002	0.001	IMIS-2D	6.3909 × 10⁻⁴	2.9020 × 10⁻⁴	4	yes
		IMCNS-2D	6.3907 × 10⁻⁴	2.9022 × 10⁻⁴	4	yes
		IMBCS-2D	8.7514 × 10⁻⁴	3.4776 × 10⁻⁴	4	yes
0.01	0.001	IMIS-2D	3.8821 × 10⁻³	1.4497 × 10⁻³	4	yes
		IMCNS-2D	fails	fails		no
		IMBCS-2D	fails	fails		no

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Atanasov, A.Z.; Koleva, M.N.; Vulkov, L.G. Inverse Problem Numerical Analysis of Forager Bee Losses in Spatial Environment without Contamination. Symmetry 2023, 15, 2099. https://doi.org/10.3390/sym15122099

AMA Style

Atanasov AZ, Koleva MN, Vulkov LG. Inverse Problem Numerical Analysis of Forager Bee Losses in Spatial Environment without Contamination. Symmetry. 2023; 15(12):2099. https://doi.org/10.3390/sym15122099

Chicago/Turabian Style

Atanasov, Atanas Z., Miglena N. Koleva, and Lubin G. Vulkov. 2023. "Inverse Problem Numerical Analysis of Forager Bee Losses in Spatial Environment without Contamination" Symmetry 15, no. 12: 2099. https://doi.org/10.3390/sym15122099

APA Style

Atanasov, A. Z., Koleva, M. N., & Vulkov, L. G. (2023). Inverse Problem Numerical Analysis of Forager Bee Losses in Spatial Environment without Contamination. Symmetry, 15(12), 2099. https://doi.org/10.3390/sym15122099

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Inverse Problem Numerical Analysis of Forager Bee Losses in Spatial Environment without Contamination

Abstract

1. Introduction

2. Model Problem

3. The Inverse Problem

4. Iterative Method for Solving a 1D Inverse Problem

4.1. Discrete Schemes

4.2. Iterative Procedure

5. Iterative Method for Solving a 2D Inverse Problem

5.1. Discrete Schemes

5.2. Iterative Method

6. Numerical Tests

7. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI