Solving direct and inverse problems for Fredholm-type integro-differential equations with application to pollution diffusion modeling

Mathematics and Computers in Simulation 223 (2024) 394–404 Available online 22 April 2024 0378-4754/© 2024 The Author(s). Published by Elsevier B.V. on behalf of International Association for Mathematics and Computers in Simulation (IMACS). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Mathematics and Computers in Simulation journal homepage: www.elsevier.com/locate/matcom Original articles Solving direct and inverse problems for Fredholm-type integro-differential equations with application to pollution diffusion modeling M.I. Berenguer a,∗, D. Gámez a, H. Kunze b, D. La Torrec, M. Ruiz Galána aDepartment of Applied Mathematics and Institute of Mathematics (IMAG), University of Granada, Granada, Spain bDepartment of Mathematics and Statistics, University of Guelph, Guelph, Canada cSKEMA Business School, Université Côte d’Azur, France ARTICLE INFO Keywords: Fredholm-type integro-differential equations Direct problem Inverse problem Pollution diffusion modeling ABSTRACT We analyze a particular Fredholm-type partial integro-differential equation. We study the direct problem and prove existence and uniqueness of the solution via a fixed-point argument for generalized contractive maps. This approach also allows us to formulate a collage-type result that can be used to solve inverse problems. We provide numerical examples and we also show how these equations can be used to model pollution diffusion of heavy pollutants and non-volatile substances such as heavy metals, chemical spills, radioactive isotopes, and others. 1. Introduction Integro-differential equations are expressions that involve the derivative of an unknown function, which also appears under an integral sign. They have been used for several applications in numerous real-world problems (see [1,2] and the references therein). Fredholm integro-differential equations are usually defined in such a way that the integration limits are constant. They arise in various areas of science and engineering, including physics, biology, and finance and they are particularly useful in modeling phenomena with non-local interactions where the behavior of the system at a point depends on the behavior of the system over an interval [16,20,32,36]. In this paper we focus on the analysis of the following problem, which belongs to the family of Fredholm-type integro-differential equations: ⎧ ⎪ ⎨ ⎪ ⎩ 𝑢𝑡(𝑥, 𝑡) = 𝑔(𝑥, 𝑡)𝑢(𝑥, 𝑡) + ∫𝑏 𝑎 𝑓(𝑥, 𝑠, 𝑢(𝑠, 𝑡))𝑑𝜇(𝑠),(𝑥, 𝑡)∈[𝑎, 𝑏] × [0, 𝑇 ] 𝑢(𝑥, 0) = 𝑢0(𝑥), 𝑥 ∈ [𝑎, 𝑏] .(1) where 𝑢𝑡(𝑥, 𝑡)denote the partial derivative of 𝑢with respect to 𝑡, the functions 𝑓and 𝑔are given and 𝑢0is the initial condition and we also suppose that 𝜇is a generic probability measure with compact support over an interval [𝑎, 𝑏]. Our interest for this type of equations is motivated by the extensive use of them in pollution diffusion modeling. In that setting, the probability measure 𝜇describes the uncertainty due to unexpected emissions of pollutant or it describes the relative importance and effect of a certain pollutant in a specific location. This family of equations is of particular importance when modeling heavy pollutants or non-volatile substances (heavy metals, non-volatile organic compounds, oil spills, chemical spills, radioactive isotopes ∗Corresponding author. E-mail address: [email protected] (M.I. Berenguer). https://doi.org/10.1016/j.matcom.2024.04.021 Received 24 November 2023; Received in revised form 10 April 2024; Accepted 19 April 2024 Mathematics and Computers in Simulation 223 (2024) 394–404 395 M.I. Berenguer et al. such as uranium, thorium, radium, and radon) that tend to remain localized in their original states or forms, rather than dispersing into the surrounding atmosphere [35,47]. These pollutants can pose significant environmental and health risks, particularly if they contaminate soil, water bodies, or ecosystems, necessitating careful management and remediation efforts. Such pollutants are often associated with industrial processes, and the contamination process is quite slow and it cannot be described by a classical diffusion term. One can refer to [38] to understand various aspects of atmospheric chemistry and physics of heavy pollutants. In [33,45] the authors provide an overview of environmental modeling techniques, including diffusion modeling, for the fate and transport of pollutants in different environmental media while [6] covers the fundamental principles of fluid flow and transport in porous media, which are relevant to modeling the movement of pollutants in groundwater and soil (see also [15]). We will show how the previous integro-differential problem (1) can be stated in terms of a fixed-point equation and, for this abstract problem, we will discuss both the direct and the inverse problems. The direct problem for a certain equation consists in proving the existence and uniqueness of a solution as well as the stability of its solution. It starts with the causes and then calculates the effects. Solving a direct problem involves the following steps: Formulation of the problem, selection of a solution technique, solution process, and verification and interpretation. The inverse problem, instead, is the inverse of a direct problem [5,22,23,34,44,46] and it focuses on the estimation of the values of unknown parameters from a set of observations. That is, an inverse problem’s solution starts from the effects to calculate (or, more likely, to estimate) the causes. This is usually done by estimating the parameters from observed data or numerical approximations of the solution. Solving inverse problems often requires sophisticated mathematical techniques, optimization algorithms, and numerical simulations. Regularization methods and Bayesian approaches are commonly used to address the challenges associated with illposedness and uncertainties in the data. Inverse problems have wide applications in a variety of different fields, ranging from system identification to optics, from acoustics to communication theory, from signal processing to medical imaging, and many others. In order to analyze the direct and inverse problems for problem (1), we will reformulate the equation using a fixed-point framework. This approach has been used in several previous papers in the literature of differential, integral and integro-differential equations for dealing with both direct and inverse problems (see [7,8,19,26–28]). In this context the direct problem is normally studied by means of Banach’s fixed point Theorem (or its extensions), and it is based on a generalized contractivity property of the map and the completeness of the underlying space. The inverse problem, instead, relies on a consequence of Banach’s Theorem known as Collage Theorem. Thanks to this result, the inverse problem boils down to a finite dimensional optimization problem over the space of unknown parameters. This optimization is, most of the time, a non-convex problem and it requires the implementation of deterministic or stochastic global optimization algorithms. The paper proceeds as follows. In Section 2, we recall some well known facts on fixed point theory and we reformulate Problem (1) as a fixed point equation. Under certain hypotheses, an extended fixed-point result allows us to prove existence and uniqueness of the solution. In the same section we also present an algorithm to approximate the solution using the notion of a Schauder basis. We conduct some numerical experiments to test our results and show how the approximation method works in practical contexts. In Section 3, we make use of a generalization of the Collage Theorem to state the inverse problem in terms of an abstract optimization model. By means of practical examples and when specific families of functions are considered, we show how to solve it to estimate the unknown parameters of the model. In Section 4, we present an application of problem (1) to modeling pollution diffusion when heavy pollutants are considered. Finally, in Section 5we present some conclusions and future research avenues. 2. The direct problem: Model analysis and numerical approximation In this section we deal with existence and uniqueness of the solution to problem (1). We also derive a numerical algorithm to approximate it and we illustrate this approach by means of numerical experiments. 2.1. Fixed point formulation In this subsection we provide the theoretical foundations and the fixed-point argument that allow to establish existence and uniqueness of a solution to the above problem (1). The Banach Fixed-Point Theorem is a well known result and it is the basis for proving existence and uniqueness of fixed-point solutions. We are going to use a more general although well-known version known as Caccioppoli-Banach Fixed Point Theorem (see [3, Theorem 2.3]). Theorem 1. Let (𝑋, ‖⋅‖)be a Banach space and let 𝐓∶𝑋⟶𝑋be an operator such that there exits a sequence of nonnegative real numbers {𝜌𝑛}𝑛≥1satisfying ‖𝐓𝑛𝑢1−𝐓𝑛𝑢2‖≤𝜌𝑛‖𝑢1−𝑢2‖∀𝑢1, 𝑢2∈𝑋and the serie ∑𝑛≥1𝜌𝑛is convergent. Then 𝐓has a unique fixed point 𝑢 ∈𝑋and for all 𝑛∈Nand 𝑢∈𝑋, ‖𝐓𝑛𝑢−𝑢‖≤(∞ ∑ 𝑖=𝑛 𝜌𝑖)‖𝐓𝑢−𝑢‖(2) and in particular, 𝑢 = lim 𝑛 𝐓𝑛(𝑢).(3) Mathematics and Computers in Simulation 223 (2024) 394–404 396 M.I. Berenguer et al. We study the nonlinear problem presented below: given 𝑔∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]),𝑓∈𝐶([𝑎, 𝑏]2×R)and 𝑢0∈𝐶([𝑎, 𝑏]), find 𝑢∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]) with 𝑢𝑡∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]) such that: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑢𝑡(𝑥, 𝑡) = 𝑔(𝑥, 𝑡)𝑢(𝑥, 𝑡) + ∫𝑏 𝑎 𝑓(𝑥, 𝑠, 𝑢(𝑠, 𝑡))𝑑𝜇(𝑠) 𝑢(𝑥, 0) = 𝑢0(𝑥) .(4) We consider that 𝑓is Lipschitz in the third variable, and 𝑀is the Lipschitz constant. Regarding this equation, we define the integral operator 𝐓∶𝐶([𝑎, 𝑏] × [0, 𝑇 ]) ⟶𝐶([𝑎, 𝑏] × [0, 𝑇 ]) defined for (𝑥, 𝑡)∈[𝑎, 𝑏] × [0, 𝑇 ]and 𝑢∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]) as (𝐓𝑢)(𝑥, 𝑡) ∶= 𝑢0(𝑥) + ∫𝑡 0 𝑔(𝑥, 𝑟)𝑢(𝑥, 𝑟)𝑑𝑟 +∫𝑡 0∫𝑏 𝑎 𝑓(𝑥, 𝑠, 𝑢(𝑠, 𝑟))𝑑𝜇(𝑠)𝑑𝑟, (5) where ‘‘𝑑𝑟’’ denotes Lebesgue measure. Theorem 2. Let be 𝑔∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]),𝑢0∈𝐶([𝑎, 𝑏]) and 𝑓∈𝐶([𝑎, 𝑏]2×R)satisfying a Lipschitz condition with respect to the last variable, with Lipschitz constant 𝑀. Then, problem (4) admits a unique solution. Proof. Let 𝐓∶𝐶([𝑎, 𝑏] × [0, 𝑇 ]) ⟶𝐶([𝑎, 𝑏] × [0, 𝑇 ]) be the operator defined by (5). We first proof that 𝑢 is solution of problem (4) if, and only if 𝐓𝑢 =𝑢. If 𝑢 is solution of (4) taking into account that 𝑢 is continuous, by the Fundamental Theorem of Calculus we have that 𝑢 is a fixed point of the operator 𝐓. And conversely, if 𝐓𝑢 =𝑢, as 𝑔,𝑢 are continuous, then the Fundamental Theorem of Calculus yields 𝑢𝑡(𝑥, 𝑡) = 𝑔(𝑥, 𝑡)𝑢(𝑥, 𝑡) + ∫𝑏 𝑎 𝑓(𝑥, 𝑠, 𝑢(𝑠, 𝑡))𝑑𝜇(𝑠), and since 𝑢(𝑥, 0) = 𝑢0(𝑥),𝑢 is solution of (4). Therefore, to conclude the proof we will show that 𝐓has a unique fixed point. To this end, if 𝑢1, 𝑢2∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]), (𝑥, 𝑡)∈[𝑎, 𝑏] × [0, 𝑇 ]then, taking 𝑁= max (𝑥,𝑡)∈[𝑎,𝑏]×[0,𝑇 ]|𝑔(𝑥, 𝑡)|, we have by induction that for 𝑛≥1 |𝐓𝑛𝑢1(𝑥, 𝑡) − 𝐓𝑛𝑢2(𝑥, 𝑡)|≤(𝑡(𝑁+𝑀))𝑛 𝑛!‖𝑢1−𝑢2‖. More specifically, for 𝑛= 1, |𝐓𝑢1(𝑥, 𝑡) − 𝐓𝑢2(𝑥, 𝑡)|≤∫𝑡 0|𝑔(𝑥, 𝑟)|⋅|𝑢1(𝑥, 𝑟) − 𝑢2(𝑥, 𝑟)|𝑑𝑟 +∫𝑡 0∫𝑏 𝑎|𝑓(𝑥, 𝑠, 𝑢1(𝑠, 𝑟)) − 𝑓(𝑥, 𝑠, 𝑢2(𝑠, 𝑟))|𝑑𝜇(𝑠)𝑑𝑟 ≤ 𝑁‖𝑢1−𝑢2‖𝑡+𝑀‖𝑢1−𝑢2‖𝜇([𝑎, 𝑏])𝑡=𝑡(𝑁+𝑀) 1! ‖𝑢1−𝑢2‖, and if we suppose that |𝐓𝑛𝑢1(𝑥, 𝑡) − 𝐓𝑛𝑢2(𝑥, 𝑡)|≤𝑡𝑛(𝑁+𝑀)𝑛 𝑛!‖𝑢1−𝑢2‖, then |𝐓𝑛+1𝑢1(𝑥, 𝑡) − 𝐓𝑛+1𝑢2(𝑥, 𝑡)|≤ ∫𝑡 0|𝑔(𝑥, 𝑟)|⋅|𝐓𝑛𝑢1(𝑥, 𝑟) − 𝐓𝑛𝑢2(𝑥, 𝑟)|𝑑𝑟 +∫𝑡 0∫𝑏 𝑎|𝑓(𝑥, 𝑠, 𝐓𝑛𝑢1(𝑠, 𝑟)) − 𝑓(𝑥, 𝑠, 𝐓𝑛𝑢2(𝑠, 𝑟))|𝑑𝜇(𝑠)𝑑𝑟 ≤ ≤𝑁∫𝑡 0 𝑟𝑛(𝑁+𝑀)𝑛 𝑛!‖𝑢1−𝑢2‖𝑑𝑟 +𝑀∫𝑡 0∫𝑏 𝑎 𝑟𝑛(𝑁+𝑀)𝑛 𝑛!‖𝑢1−𝑢2‖𝑑𝜇(𝑠)𝑑𝑟 = 𝑡𝑛+1 (𝑛+ 1)! 𝑁(𝑁+𝑀)𝑛‖𝑢1−𝑢2‖+𝑡𝑛+1 (𝑛+ 1)! 𝑀(𝑁+𝑀)𝑛‖𝑢1−𝑢2‖𝜇([𝑎, 𝑏]) = 𝑡𝑛+1(𝑁+𝑀)𝑛+1 (𝑛+ 1)! ‖𝑢1−𝑢2‖. As consequence, for 𝑛≥1, ‖𝐓𝑛𝑢1−𝐓𝑛𝑢2‖≤𝜌𝑛‖𝑢1−𝑢2‖(6) where 𝜌𝑛=(𝑇(𝑁+𝑀))𝑛 𝑛!.(7) Since the series ∑ 𝑛≥1 𝜌𝑛is convergent, by Theorem 1, there exists a unique 𝑢 ∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]) such that 𝐓(𝑢) = 𝑢, and for each 𝑢∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]) and 𝑛≥1, ‖𝐓𝑛𝑢−𝑢‖≤(∞ ∑ 𝑖=𝑛 𝜌𝑖)‖𝐓𝑢−𝑢‖,(8) and thus lim 𝑛→∞‖𝐓𝑛𝑢−𝑢‖= 0.□ Mathematics and Computers in Simulation 223 (2024) 394–404 397 M.I. Berenguer et al. In view of the above, the solution can be determined as the limit of the iterates of operator 𝐓. Since the calculations of such iterates from a practical point of view are not always feasible, in the next section we propose a method to approximate the solution. 2.2. Numerical method For our purposes, let us recall that a Schauder basis of a Banach space 𝑋is a sequence {𝛾𝑛}𝑛≥1⊂ 𝑋 such that ∀𝛾∈𝑋there is a unique sequence {𝜆𝑛}𝑛≥1of real numbers satisfying 𝛾=∑𝑛≥1𝜆𝑛𝛾𝑛. The sequences of projections 𝑃𝑛∶𝑋→𝑋and coordinate functionals {𝛾∗ 𝑛}𝑛≥1in 𝑋∗, are defined, respectively, as follows: 𝑃𝑛(∑ 𝑛≥1 𝜆𝑛𝛾𝑛)= 𝑛 ∑ 𝑘=1 𝜆𝑘𝛾𝑘and 𝛾∗ 𝑛(∑ 𝑛≥1 𝜆𝑛𝛾𝑛)=𝜆𝑛 and it is verified that lim 𝑛‖𝑃𝑛(𝛾) − 𝛾‖= 0 for all 𝛾∈𝑋. (9) Let {𝐴𝑛}𝑛≥1be a Schauder basis in 𝐶([𝑎, 𝑏] × [0, 𝑇 ]) with {𝑄𝑛}𝑛≥1as their associated projections and {𝐵𝑛}𝑛≥1a Schauder basis in 𝐶([𝑎, 𝑏]2× [0, 𝑇 ]), with its associated projections {𝑅𝑛}𝑛≥1. For 𝑝∈N, we define 𝐆𝑝∶𝐶([𝑎, 𝑏] × [0, 𝑇 ]) ⟶𝐶([𝑎, 𝑏] × [0, 𝑇 ]): (𝐆𝑝𝑣)(𝑥, 𝑡) ∶= 𝑢0(𝑥) + ∫𝑡 0 𝑄𝑛𝑝(𝑔⋅𝑣)(𝑥, 𝑟)𝑑𝑟 +∫𝑡 0∫𝑏 𝑎 𝑅𝑛∗ 𝑝(𝐿0(𝑣))(𝑥, 𝑠, 𝑟)𝑑𝜇(𝑠)𝑑𝑟 (10) where 𝑛𝑝∈N, 𝑛∗ 𝑝∈N,𝑣∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]) and 𝐿0∶𝐶([𝑎, 𝑏] × [0, 𝑇 ]) →𝐶([𝑎, 𝑏]2× [0, 𝑇 ]), 𝐿0(𝑣)(𝑥, 𝑠, 𝑡) ∶= 𝑓(𝑥, 𝑠, 𝑣(𝑠, 𝑡)),(𝑥, 𝑠, 𝑡)∈[𝑎, 𝑏]2× [0, 𝑇 ].(11) Let 𝑢 ∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]), and we define: 𝑧0(𝑥, 𝑡) ∶= 𝑢(𝑥, 𝑡)(12) and for 𝑚∈N, 𝑧𝑚(𝑥, 𝑡) ∶= 𝐆𝑚◦...◦𝐆1(𝑢)(𝑥, 𝑡).(13) The following result allows us to choose 𝑛1,…, 𝑛𝑚, 𝑛∗ 1,…, 𝑛∗ 𝑚∈Nso that 𝑧𝑚approximates the unique solution of (4). Theorem 3. Let 𝑢 ∈𝐶([𝑎, 𝑏] × [0, 𝑇 ]) and let 𝑢 be the unique solution of problem (1), then for each 𝜀 > 0there exist 𝑚≥1and 𝑛1,…, 𝑛𝑚, 𝑛∗ 1,…, 𝑛∗ 𝑚≥1such that ‖𝑢 −𝑧𝑚‖< 𝜀 with 𝑧0and 𝑧𝑚defined by (12) and (13). Proof. If we consider 𝜌0= 1 and for all 𝑛≥1,𝜌𝑛is given by (7), then we have that for all 𝑚≥0, ‖𝐓𝑚𝑢 −𝑧𝑚‖≤ ‖𝐓𝑚𝑢 −𝐓𝑚−1𝑧1‖+‖𝐓𝑚−1𝑧1−𝐓𝑚−2𝑧2‖+⋯+‖𝐓2𝑧𝑚−2 −𝐓𝑧𝑚−1‖+‖𝐓𝑧𝑚−1 −𝑧𝑚‖≤ (𝑚−1 ∑ 𝑝=1 ‖𝐓𝑚−𝑝+1𝑧𝑝−1 −𝐓𝑚−𝑝𝑧𝑝‖)+‖𝐓𝑧𝑚−1 −𝑧𝑚‖≤ (𝑚−1 ∑ 𝑝=1 ‖𝐓𝑚−𝑝𝐓𝑧𝑝−1 −𝐓𝑚−𝑝𝑧𝑝‖)+𝜌0‖𝐓𝑧𝑚−1 −𝑧𝑚‖≤ (𝑚−1 ∑ 𝑝=1 𝜇𝑚−𝑝‖𝐓𝑧𝑝−1 −𝑧𝑝‖)+𝜌0‖𝐓𝑧𝑚−1 −𝑧𝑚‖= 𝑚 ∑ 𝑝=1 𝜌𝑚−𝑝‖𝐓𝑧𝑝−1 −𝑧𝑝‖. Therefore, ‖𝐓𝑚𝑢 −𝑧𝑚‖≤𝑚 ∑ 𝑝=1 𝜌𝑚−𝑝‖𝐓𝑧𝑝−1 −𝑧𝑝‖. For 𝑝∈ {1,…, 𝑚}, we define 𝐿𝑝∶𝐶([𝑎, 𝑏] × [0, 𝑇 ]) →𝐶([𝑎, 𝑏]2× [0, 𝑇 ]), 𝐿𝑝(𝑣)(𝑥, 𝑠, 𝑡) ∶= 𝑓(𝑥, 𝑠, 𝐆𝑝◦...◦𝐆1(𝑣)(𝑥, 𝑡)),(𝑥, 𝑠, 𝑡)∈[𝑎, 𝑏]2× [0, 𝑇 ].(14) Then |𝐓𝑧𝑝−1(𝑥, 𝑡) − 𝑧𝑝(𝑥, 𝑡)|≤ ∫𝑡 0|||(𝑔⋅𝑧𝑝−1)(𝑥, 𝑟) − 𝑄𝑛𝑝(𝑔⋅𝑧𝑝−1)(𝑥, 𝑟)|||𝑑𝑟 +∫𝑡 0∫𝑏 𝑎|||𝑓(𝑥, 𝑠, 𝑧𝑝−1(𝑠, 𝑟)) − 𝑅𝑛∗ 𝑝(𝐿𝑝−1(𝑢))(𝑥, 𝑠, 𝑟)|||𝑑𝜇(𝑠)𝑑𝑟 Mathematics and Computers in Simulation 223 (2024) 394–404 398 M.I. Berenguer et al. Table 1 Absolute error |𝑧7(𝑥, 𝑡) − 𝑢(𝑥, 𝑡)|for Example 2.1 with 𝑛𝑖= 92and 𝑛∗ 𝑖= 93. 𝑡,𝑥0.2 0.4 0.6 0.8 1 0.2 1.03155 × 10-3 2.52247 × 10-3 4.62988 × 10-3 7.65072 × 10-3 1.18063 × 10-2 0.4 1.67067 × 10-3 4.06346 × 10-3 7.42057 × 10-3 1.21464 × 10-2 1.86538 × 10-2 0.6 1.92418 × 10-3 4.61077 × 10-3 8.2941 × 10-3 1.33447 × 10-2 2.01397 × 10-2 0.8 1.77083 × 10-3 4.12469 × 10-3 7.21951 × 10-3 1.13068 × 10-2 1.66163 × 10-2 11.32592 × 10-3 2.97533 × 10-3 5.03546 × 10-3 7.65843 × 10-3 1.09648 × 10-2 Table 2 Absolute error |𝑧7(𝑥, 𝑡) − 𝑢(𝑥, 𝑡)|for Example 2.1 with 𝑛𝑖= 172and 𝑛∗ 𝑖= 173. 𝑡,𝑥0.2 0.4 0.6 0.8 1 0.2 4.72292 × 10-4 1.14505 × 10-3 2.08206 × 10-3 3.36593 × 10-3 5.10239 × 10-3 0.4 3.68087 × 10-4 8.85775 × 10-4 1.59959 × 10-3 2.57172 × 10-3 3.88733 × 10-3 0.6 4.47135 × 10-4 1.07321 × 10-3 1.93185 × 10-3 3.09301 × 10-3 4.64846 × 10-3 0.8 1.3224 × 10-3 3.0642 × 10-3 5.33105 × 10-3 8.25494 × 10-3 1.2 × 10-2 15.04882 × 10-4 1.14554 × 10-3 1.95615 × 10-3 2.9809 × 10-3 4.27573 × 10-3 ≤𝑇‖𝑔⋅𝑧𝑝−1 −𝑄𝑛𝑝(𝑔⋅𝑧𝑝−1)‖+𝑇 𝜇([𝑎, 𝑏])‖𝐿𝑝−1(𝑢) − 𝑅𝑛∗ 𝑝(𝐿𝑝−1(𝑢))‖, therefore, ‖𝐓𝑚𝑢 −𝑧𝑚‖≤𝑚 ∑ 𝑝=1 𝑇 𝜌𝑚−𝑝(‖𝑔⋅𝑧𝑝−1 −𝑄𝑛𝑝(𝑔⋅𝑧𝑝−1)‖+‖𝐿𝑝−1(𝑢) − 𝑅𝑛∗ 𝑝(𝐿𝑝−1(𝑢))‖).(15) For all 𝜀 > 0,by(3), there exits 𝑚∈Nverifying that ‖𝑢 −𝐓𝑚𝑢‖< 𝜀∕2.By(9) and (15), there exist 𝑛1,…, 𝑛𝑚, 𝑛∗ 1,…, 𝑛∗ 𝑚≥1such that ‖𝐓𝑚𝑢 −𝑧𝑚‖< 𝜀∕2. Then, by triangle inequality, ‖𝑢 −𝑧𝑚‖≤‖𝐓𝑚𝑢 −𝑢‖+‖𝐓𝑚𝑢 −𝑧𝑚‖<𝜀 2+𝜀 2< 𝜀. □ 2.3. Numerical experiments The presented numerical method has been tested on two numerical examples. Given the dense sequence {𝑠𝑛}𝑛≥1(respectively {𝑡𝑛}𝑛≥1) taken on the diadic partition of the interval [𝑎, 𝑏]([0, 𝑇 ], respectively) with 𝑠1=𝑎and 𝑠2=𝑏(𝑡1= 0 and 𝑡2=𝑇respectively), we consider the Faber Schauder basis {𝜑𝑛}𝑛≥1in 𝐶([𝑎, 𝑏]) and {𝜓𝑛}𝑛≥1in 𝐶([0, 𝑇 ]), and we fix the concrete Schauder basis {𝐴𝑛}𝑛≥1in 𝐶([𝑎, 𝑏]×[0, 𝑇 ]) and {𝐵𝑛}𝑛≥1in 𝐶([𝑎, 𝑏]2×[0, 𝑇 ]) as the tensor product of bases {𝜑𝑛}𝑛≥1and {𝜓𝑛}𝑛≥1[39]. The advantage of introducing these bases is that, via the following change-of-variable, the calculations are feasible (see [25, Section 5]): Let us suppose that 𝜇is a measure defined on the Borel subsets of Rand suppose that supp(𝜇) = [𝑎, 𝑏]and that is non-atomic. Let 𝐹∶𝐾= supp(𝜇)→[0,1] be the cumulative of 𝜇and 𝐹−1 ∶ [0,1] →𝐾be its inverse. Then, given a function 𝑓∶𝐾→R, the following change of variable rule holds: ∫𝐾 𝑓(𝑥)𝑑𝜇(𝑥) = ∫1 0 𝑓(𝐹−1(𝑥))𝑑𝑥 (16) where 𝑑𝑥 is the Lebesgue measure. The computations to implement the algorithm for the approximate solution have been done using Mathematica. Example 2.1. Consider the following problem ⎧ ⎪ ⎨ ⎪ ⎩ 𝑢𝑡(𝑥, 𝑡) = 2 + 𝑒𝑥−2𝑥2−𝑒𝑥𝑥2 2𝑢(𝑥, 𝑡) + ∫1 0 𝑥2𝑒−𝑠𝑢(𝑠, 𝑡)𝑑𝑠 𝑢(𝑥, 0) = 𝑒−𝑥(𝑥, 𝑡) ∈ [0,1] × [0,1] .(17) where the exact solution is given by 𝑢(𝑥, 𝑡) = 𝑒𝑡−𝑥. To compare numerical and exact solutions see Tables 1 and 2where it is shown the absolute errors |𝑧7(𝑥, 𝑡) − 𝑢(𝑥, 𝑡)|taking, on one hand, 𝑛𝑖= 92and 𝑛∗ 𝑖= 93, and in other hand 𝑛𝑖= 172and 𝑛∗ 𝑖= 173for 𝑖= 1,…,7. The sup norm of the error ‖𝑧7−𝑢‖is 2.27689 × 10−2 in the first case and 1.66485 × 10−2 in the second one. Mathematics and Computers in Simulation 223 (2024) 394–404 399 M.I. Berenguer et al. Table 3 Absolute error |𝑧3(𝑥, 𝑡) − 𝑢(𝑥, 𝑡)|for Example 2.2 with 𝑛𝑖= 92and 𝑛∗ 𝑖= 93. 𝑡,𝑥0.2 0.4 0.6 0.8 1 0.2 1.51415 × 10−4 3.01377 × 10−4 4.43192 × 10−4 6.06209 × 10−4 1.61912 × 10−3 0.4 3.04166 × 10−4 6.16909 × 10−4 9.16421 × 10−4 1.04999 × 10−3 7.86755 × 10−4 0.6 4.5769 × 10−4 9.33775 × 10−4 1.41162 × 10−3 1.70325 × 10−3 1.22657 × 10−3 0.8 6.4003 × 1041.30108 × 10−3 1.98926 × 10−3 2.58337 × 10−3 2.4922 × 10−3 17.50114 × 10−4 1.55333 × 10−3 2.51532 × 10−3 3.74581 × 10−3 4.89973 × 10−3 Example 2.2. Let us suppose that the density 𝜉of the measure 𝜇over the interval [0,1] is equal to 𝜉(𝑥) = 3 4+𝑥 2and 0otherwise and consider the problem ⎧ ⎪ ⎨ ⎪ ⎩ 𝑢𝑡(𝑥, 𝑡) = 𝑔(𝑥, 𝑡)𝑢(𝑥, 𝑡) + ∫1 0 𝑥𝑢(𝑠, 𝑡)𝑑𝜇(𝑠) 𝑢(𝑥, 0) = cos(𝑥) (𝑥, 𝑡) ∈ [0,1] × [0,1] .(18) where 𝑔(𝑥, 𝑡)is such that 𝑢(𝑥, 𝑡) = cos(𝑥−𝑡)is the solution. The obtained results are presented for 𝑚= 3 with 𝑛𝑖= 92and 𝑛∗ 𝑖= 93for 𝑖= 1,2,3can be viewed in Table 3. The sup norm of the error ‖𝑧3−𝑢‖is 4.89973 × 10−3. 3. The inverse problem The inverse problem consists in determining an estimation of the unknown parameters of the model starting from data that have been gathered by samples or experiments. In the case of our model Eq. (1), the data could be values of 𝑢at particular observation points (𝑥𝑖, 𝑡𝑖), and the parameters we wish to estimate might be the coefficients in the proposed functional forms of 𝑔and 𝑓. This problem can be solved by using a Collage-type Theorem. The classical ‘‘Collage Theorem’’ for the standard version of Banach’s fixed point theorem has been known for many years, but it was first used in the context of solving an inverse problem in [5] for the analysis of inverse problems in fractal image compression. This classical Collage Theorem states an approximation result between any target object 𝑢and the solution 𝑢 of a fixed point equation 𝑇 𝑢 =𝑢. In the fractal imaging context, the action of the fractal transform operator 𝑇on an image 𝑢produces a new image 𝑣that is an assemblage of adjusted subsets of 𝑢; it is this action that gives rise to name ‘‘Collage’’. The Collage Theorem is a crucial result for solving inverse problems for fixed point equations and it has been used in a variety of different formulations and frameworks (see [4,9,10,12–14,17,18,37,40,41]). In the current paper, the next result Theorem 4 allows us to formulate a computable inverse problem approach to problem (1). The proof of this result is similar to that in [24], which is based on the application of the standard Collage Theorem to a contraction map 𝐓𝑛0. Theorem 4 ([24]).Let (𝑋, ‖⋅‖)be a Banach space and let 𝐓∶𝑋⟶𝑋be an operator such that there exits a sequence of nonnegative real numbers {𝜌𝑛}𝑛≥1satisfying ‖𝐓𝑛𝑢1−𝐓𝑛𝑢2‖≤𝜌𝑛‖𝑢1−𝑢2‖∀𝑢1, 𝑢2∈𝑋and the series ∑𝑛≥1𝜌𝑛is convergent. Let 𝑢 be the unique fixed point of 𝐓, as guaranteed by Theorem 1, and 𝑢be a chosen element of 𝑋. Then there exists 𝑛0such that 𝜌𝑛0<1and this implies ‖𝑢−𝑢‖≤∑𝑛0−1 𝑖=0 𝜌𝑖 1 − 𝜌𝑛0‖𝐓𝑢−𝑢‖.(19) When we label the observational data or their interpolant as 𝑢, we can identify the right-hand side of (19) as the error in approximating 𝑢by the fixed point 𝑢 of 𝐓. If 𝜌𝑛0≪1, then the magnitude of the error can be controlled by minimizing ‖𝐓𝑢−𝑢‖. If, as usually happens in this context, the operator 𝐓belongs to a fixed family of functions and it is parametrized by a vector 𝜆∈𝛬 ⊂ R𝑝 (in the inverse problem framework the target element 𝑢∈𝑋is fixed and known), the minimization of the term 𝛥(𝜆) ∶= ‖𝐓𝜆𝑢−𝑢‖ boils down to solving the optimization problem: min 𝜆∈𝛬‖𝐓𝜆𝑢−𝑢‖(20) which has at least one solution whenever 𝛥is continuous and 𝛬is a compact subset of R𝑝. Now let us apply the previous abstract inverse problem approach to our specific model. We suppose we know or have an estimate of the target solution 𝑢to the model and we want to determine an estimation of 𝑔and 𝑓. Using the Collage Theorem, this can be stated in terms of the minimization of the Collage distance: min 𝑓,𝑔 ‖𝑢−𝐓𝑢‖=‖‖‖‖‖𝑢(𝑥, 𝑡) − 𝑢0(𝑥) − ∫𝑡 0 𝑔(𝑥, 𝑟)𝑢(𝑥, 𝑟)𝑑𝑟 −∫𝑡 0∫𝑏 𝑎 𝑓(𝑥, 𝑠, 𝑢(𝑠, 𝑟))𝑑𝜇(𝑠)𝑑𝑟‖‖‖‖‖,(21) where the minimization is done over all Lipschitz functions 𝑓and 𝑔. The complexity of this optimization problem depends on the assumptions of 𝑓and 𝑔. A regularization term [42,43], expressed in terms of the zero-norm of 𝑓and 𝑔, can also be added to reduce the dimensionality. Mathematics and Computers in Simulation 223 (2024) 394–404 400 M.I. Berenguer et al. Table 4 Results from Example 3.1. 𝑁𝑥𝑁𝑡𝐷 𝜀 𝛥 True values 1.0000000 −0.4323324 −0.4323324 −0.2161662 𝑔0𝑔2𝑔3𝑔4 20 10 6 0.00 0.0000000 0.9996606 −0.5236164 −0.1846185 −0.4777624 20 10 10 0.00 0.0000000 0.9996503 −0.5218805 −0.1879165 −0.4759510 100 50 6 0.00 0.0000000 0.9996606 −0.5236164 −0.1846185 −0.4777624 100 50 10 0.00 0.0000000 0.9996503 −0.5218805 −0.1879165 −0.4759510 20 10 6 0.01 0.0000037 1.0009614 −0.5334301 −0.3618969 −0.3097484 20 10 10 0.01 0.0000078 1.0037317 −0.2392890 −0.6389388 −0.2441403 100 50 6 0.01 0.0000037 1.0009614 −0.5334301 −0.3618969 −0.3097484 100 50 10 0.01 0.0000078 1.0037317 −0.2392890 −0.6389388 −0.2441403 20 10 6 0.05 0.0000619 0.9953297 −1.3535893 0.8862019 −0.7843862 20 10 10 0.05 0.0001412 1.0012706 −0.5801000 −1.1627976 0.3364531 100 50 6 0.05 0.0000619 0.9953297 −1.3535893 0.8862019 −0.7843862 100 50 10 0.05 0.0001411 1.0012706 −0.5801000 −1.1627976 0.3364531 As mentioned above, in a practical inverse problem, both functions 𝑓and 𝑔may be assumed to take specific forms. Therefore, the essence of the inverse problem lies in estimating the unknown coefficients within these predefined forms, namely 𝑔𝑖,𝑖= 1,…, 𝑛𝑔, and 𝑓𝑗,𝑗= 1,…, 𝑛𝑓. The inverse problem boils down to the following form: min 𝑔𝑖,𝑓𝑗‖𝑢−𝐓𝑢‖.(22) The following examples show a practical implementations of the Collage Theorem. Example 3.1. We return to (17) considered in Example 2.1. We sample the solution 𝑒𝑡−𝑥on [0,1]2on a uniform grid with 𝑁𝑥(𝑁𝑡) nodes in the 𝑥- (𝑡-) direction, add Gaussian noise with amplitude 𝜀to the sample values, and then fit a polynomial 𝑢(𝑥, 𝑡)of degree 𝐷to the noised data. We seek an equation of the form 𝑢𝑡(𝑥, 𝑡) = 𝑔(𝑥)𝑢(𝑥, 𝑡) + ∫1 0 𝑥2𝑒−𝑠𝑢(𝑠, 𝑡)𝑑𝑠, (23) supposing that 𝑔has the form 𝑔(𝑥) = 𝑔0+𝑔2𝑥2+𝑔3𝑥3+𝑔4𝑥4. The resulting collage distance ‖𝑢−𝐓𝑢‖, with 𝐓the associated integral operator and 𝑢our fitted target, is a function of the 𝑔𝑖. The results of minimizing the collage distance in various scenarios are presented in Table 4. Writing the first terms of the Taylor expansion of the true function 𝑔𝑡𝑟𝑢𝑒(𝑥)about 𝑥= 0, we note that near 𝑥= 0, 𝑔𝑡𝑟𝑢𝑒(𝑥)≈1.0000000 − 0.4323324𝑥2− 0.4323324𝑥3− 0.2161662𝑥4, with coefficients to 7 decimal places. We see that the quality of the results is alright when there is no noise or very low noise, but the worsens as the noise level grows to 5%. Example 3.2. Let us suppose that [𝑎, 𝑏] = [−1,1],𝑔(𝑥, 𝑡)is linear and takes the form 𝑔(𝑥, 𝑡) = 𝑔0+𝑔1𝑥+𝑔2𝑡(24) and the function 𝑓, instead, can be written as 𝑓(𝑥, 𝑠, 𝑢(𝑥, 𝑡)) = 𝜙(𝑥, 𝑠)𝑢(𝑠, 𝑡)(25) where 𝜙is a normal kernel taking the form 𝜙(𝑥, 𝑠) =  𝜙(𝑥)𝑒−1 1−𝑠2.(26) Let us suppose that the density 𝜉of the measure 𝜇over the interval [−1,1] is equal to 𝜉(𝑥) = 𝑥+1 2and 0otherwise. Then the cumulative 𝐹(𝑥) = 𝜇([−1, 𝑥]) is equal to 𝑥2 4+𝑥 2+1 4over [−1,1],0when 𝑥≤−1 and 1when 𝑥≥1. The inverse of the cumulative 𝑥=𝐹−1(𝑦) = √4𝑦− 1 and then the operator 𝐓𝑢can be written as (𝐓𝑢)(𝑥, 𝑡) = 𝑢0(𝑥) + ∫𝑡 0 𝑔(𝑥, 𝑟)𝑢(𝑥, 𝑟)𝑑𝑟 +∫𝑡 0∫1 −1 𝜙(𝑥, 𝑠)𝑢(𝑠, 𝑡)𝑑𝜇(𝑠)𝑑𝑟 =𝑢0(𝑥) + ∫𝑡 0 (𝑔0+𝑔1𝑥+𝑔2𝑟)𝑢(𝑥, 𝑟)𝑑𝑟 + 𝜙(𝑥)∫𝑡 0∫1 −1 𝑒−1 1−𝑠2𝑢(𝑠, 𝑟)𝑑𝜇(𝑠)𝑑𝑟 Mathematics and Computers in Simulation 223 (2024) 394–404 401 M.I. Berenguer et al. Fig. 1. Solution for Example 3.2, (a) profiles and (b) solution surface for 0≤𝑡≤1. Table 5 Results for Example 3.2. True values are (𝑔0, 𝑔1, 𝑔2) = (0.5,2,−1). 𝑁𝑥𝑁𝑡𝜀 𝛥 𝑔0𝑔1𝑔2 20 10 0.00 0.0001917 0.4950936 2.0116270 −1.0357785 50 25 0.00 0.0001518 0.4948047 2.0129139 −1.0356112 100 50 0.00 0.0001435 0.4946804 2.0133074 −1.0354633 20 10 0.01 0.0002193 0.4938846 2.0039143 −1.0298406 50 25 0.01 0.0001650 0.4869434 2.0167247 −1.0247163 100 50 0.01 0.0001622 0.5015758 2.0094539 −1.0446763 20 10 0.05 0.0009047 0.5205707 1.9106357 −1.0320081 50 25 0.05 0.0002995 0.5208807 2.0356040 −1.0885161 100 50 0.05 0.0002057 0.4713763 2.0317943 −1.0070384 20 10 0.10 0.0021695 0.4923433 1.9702788 −1.0065235 50 25 0.10 0.0009982 0.4551648 2.0330032 −0.9799383 100 50 0.10 0.0002851 0.4853762 2.0299821 −1.0296401 =𝑢0(𝑥) + ∫𝑡 0 (𝑔0+𝑔1𝑥+𝑔2𝑟)𝑢(𝑥, 𝑟)𝑑𝑟 + 𝜙(𝑥)∫𝑡 0∫1 0 𝑒−1 1−(√4𝑦−1)2𝑢(√4𝑦− 1, 𝑟)𝑑𝑦𝑑𝑟 (27) Taking  𝜙(𝑥) = 𝑒−𝑥2,𝑢0(𝑥) = 1 − 𝑥2, and, temporarily, (𝑔0, 𝑔1, 𝑔2) = (0.5,2,−1), we solve (1) numerically and the solution surface is illustrated in Fig. 1. We display the solution on a uniform grid with 𝑁𝑥(𝑁𝑡) nodes in the 𝑥- (𝑡-) direction, add Gaussian noise with amplitude 𝜀to the sample values, and then fit a polynomial 𝑢(𝑥, 𝑡)of degree 6 to the noised data. This 𝑢plays the role of the target in the Collage Theorem, as described above, with (𝐓𝑢)(𝑥, 𝑡)given in (27). The inverse problem we seek to solve is min 𝑔0,𝑔1,𝑔2 𝛥∶= ‖‖‖𝑢− (𝐓𝑔0,𝑔1,𝑔2𝑢)‖‖‖ Results are presented in Table 5. We see that the recovered values are reasonable; that improving the target by increasing 𝑁𝑥and 𝑁𝑡improves the results; and that the results are negatively, but not severely, impacted by the addition of noise. 4. An application to pollution diffusion modeling The control of pollution emissions is probably one of the most challenging issue to guarantee the long-run sustainability and growth. There is a broad research interest in modeling pollution diffusion and determine which actions can be put in place to control the emissions and implement abatement activities and optimal policies. In general, in a one-dimensional environment [𝑎, 𝑏], the level of pollution 𝑃(𝑥, 𝑡)at the time 𝑡∈ [0, 𝑇 ]and location 𝑥can be model through a reaction–diffusion equation taking the form (𝑥, 𝑡)∈[𝑎, 𝑏] × [0, 𝑇 ]: ⎧ ⎪ ⎨ ⎪ ⎩ 𝑃𝑡(𝑥, 𝑡) = 𝑐𝑃𝑥𝑥(𝑥, 𝑡) + 𝑆(𝑥, 𝑡)𝑃(𝑥, 𝑡) − 𝛿𝑃𝑃(𝑥, 𝑡) + ∫𝑏 𝑎 𝜙(𝑠, 𝑥)𝑃(𝑠, 𝑡)𝑑𝜇(𝑠) 𝑃(𝑥, 0) = 𝑃0(𝑥) (28) Mathematics and Computers in Simulation 223 (2024) 394–404 402 M.I. Berenguer et al. Table 6 Numerical results 𝑧8(𝑥, 𝑡)obtained in Example 4.1. 𝑡,𝑥0.2 0.4 0.6 0.8 1 0.2 0.256062 0.293058 0.30049 0.276818 0.229786 0.4 0.56453 0.712074 0.803802 0.813097 0.738161 0.6 0.932108 1.28944 1.58932 1.74599 1.71097 0.8 1.36636 2.0656 2.76622 3.28263 3.45529 1 1.8741 3.08456 4.46908 5.70699 6.43501 subject to an initial condition at 𝑡= 0 where 𝑃𝑡(𝑥, 𝑡)denote the partial derivative of 𝑃with respect to 𝑡. Here 𝑃𝑥𝑥 is the second-order spatial derivative of 𝑃(𝑥, 𝑡),𝑐is the diffusion coefficient, 𝑆(𝑥, 𝑡)is the growth rate of pollution and it depends on the level of emissions determined by the production output, −𝛿𝑃𝑃takes into account the natural environmental cleaning, and the integral term models the contribution of the level of pollution at each location 𝑠to the accumulation of pollution at 𝑥and the probability measure 𝜇 describes the relative importance and effect of certain pollutant in a specific location (see [11,21,29–31]). The kernel 𝜙weights the contribution of each location in terms of total emissions. When heavy pollutants or non-volatile substances such as heavy metals, non-volatile organic compounds, oil spills, chemical spills, radioactive isotopes such as uranium, thorium, radium, and radon are considered, the diffusion coefficient can be neglected and the equation boils down to ⎧ ⎪ ⎨ ⎪ ⎩ 𝑃𝑡(𝑥, 𝑡) = 𝑆(𝑥, 𝑡)𝑃(𝑥, 𝑡) − 𝛿𝑃𝑃(𝑥, 𝑡) + ∫𝑏 𝑎 𝜙(𝑠, 𝑥)𝑃(𝑠, 𝑡)𝑑𝜇(𝑠) 𝑃(𝑥, 0) = 𝑃0(𝑥) .(29) In fact, heavy pollutants tend to remain localized in their original states or forms, rather than dispersing into the surrounding atmosphere. In this case there is no flow of material associated with pollution diffusion but, instead, the level of pollution increases as a consequence of pollution stocks and water or soil contamination. We apply the numerical method described in Section 2for the following pollution model in [0,1] × [0,1]. Example 4.1. Let us suppose that the density 𝜉of the measure 𝜇over the interval [0,1] is equal to 𝜉(𝑥) = 2𝑥and 0otherwise and consider the problem with 𝑆(𝑥, 𝑡) = 2 + 2𝑥− 3𝑡,𝛿𝑃= 1,𝜙(𝑠, 𝑥) = 𝑥2cos(𝑠)and 𝑃0(𝑥) = 𝑥. The obtained results are presented for 𝑚= 8 with 𝑛𝑖= 92and 𝑛∗ 𝑖= 93for 𝑖= 1,…,8can be viewed in Table 6. That 𝑚has been chosen in such away that ‖𝑧𝑚−𝑧𝑚−1‖≤10−3. An interesting inverse problem for the above problem (29) is the following: Given gathered data of 𝑃and suppose we know both 𝛿𝑃and 𝜙, let us determine an estimation of the growth rate 𝑆. The following example shows a numerical implementation of the inverse problem with simulated data. Example 4.2. As an example inverse problem, we proceed similar to Example 3.2, choosing the same measure 𝜇, and again taking  𝜙(𝑥) = 𝑒−𝑥2and 𝜙(𝑥, 𝑠) =  𝜙(𝑥)𝑒−1 1−𝑠2. Here we set, 𝑢0(𝑥)=2𝑥6+𝑥5− 4𝑥4− 2𝑥3+3 2𝑥2+𝑥+1 2, with its graph identifying two peak pollution levels on [−1,1]. We suppose that 𝑆(𝑥, 𝑡) = 𝑆0+𝑆1𝑥+𝑆2𝑡, using the values (𝑆0, 𝑆1, 𝑆2) = (1,−2,0.8) to generate a target 𝑃, in the manner of Example 3.2. For each (𝑁𝑥, 𝑁𝑡, 𝜀)triplet, we construct two targets, one of degree 8 and the other of degree 10. The numerical solution (with no noise added) is presented in Fig. 2 We now seek to solve the corresponding inverse problem for (𝑆0, 𝑆1, 𝑆2). The results are presented in Table 7. The recovered values are close to the true values, and they hold up when a low level of noise is added. It is helpful to look again at the solution surface in Fig. 2 to recognize the challenges of using a polynomial basis for the target function. In the thirdand second-last rows of the table, for example, we see that increasing the degree of the fitted polynomial reduces the accuracy of the recovered parameter values. 5. Conclusion In this paper we have analyzed a particular Fredholm-type integro-differential equation. In the first part we have considered the direct problem and proved existence and uniqueness of the solution via a fixed-point argument. In the second part, instead, we have formulated a collage-type result that can be used to solve inverse problems. Numerical examples as well as a potential application to modeling of pollution diffusion conclude the paper. Future avenues include the analysis of more general models of macroeconomic geography in which demography and pollution dynamics affect each other. These problems are modeled with systems of two partial differential equations, one describing the evolution of pollution over time and space and the other describing the evolution of the human population.