Aliou Sonko¹, Maurice Leno¹, Aladji Babacar Niang^1,2, Ibrahima Sory Mamikouny Camara^1
1Faculty of Sciences and Techniques, University of N’Zérékoré, N’Zérékoré, Guinea.
²Imhotep International Mathematical Center, Gaston Berger University, Saint-Louis, Senegal.
DOI: 10.4236/jamp.2026.146104   PDF    HTML   XML   2 Downloads   25 Views   Citations

Abstract

Understanding how knowledge spreads within a learner population and how pedagogical methods shape this propagation is a fundamental challenge in educational science. In this paper, we propose a reaction-diffusion partial differential equation (PDE) framework to model the spatiotemporal dynamics of knowledge acquisition in a structured population. The state variable $u\left(x,t\right)\in \left[0,1\right]$ represents the local density of mastered competence at position $x$ in an abstract learner space and at time $t\ge 0$ . The diffusion term models the social contagion of knowledge through peer interactions, modulated by a set of pedagogical methods $\left(M_i\right)$ whose individual diffusion efficiencies are captured by coefficients $\left({\delta }_i\right)$ . The reaction term combines a logistic growth law, a Hill-type threshold function accounting for the triggering effect in learning, a social facilitation integral operator, and a linear forgetting term consistent with the Ebbinghaus forgetting curve. The learning rate $r\left(x,t\right)$ incorporates both the positive contributions of pedagogical methods through coefficients $\left({\alpha }_i\right)$ and the inhibitory effect of cognitive stress through a parameter $\beta$ . We establish the mathematical well-posedness of the model, derive sufficient conditions for the existence of nontrivial steady states, and perform a local stability analysis. A fictitious numerical simulation using MATLAB is conducted to validate the qualitative behavior of the model and to illustrate the comparative impact of five pedagogical methods on the long-term knowledge distribution. The extended two-compartment system coupling cognitive and non-cognitive competences is also discussed.

Keywords

Reaction-Diffusion Equation, Knowledge Propagation, Pedagogical Methods, Logistic Growth, Hill Function, Social Facilitation, Ebbinghaus Forgetting, Stability Analysis, MATLAB Simulation

Share and Cite:

1. Introduction

The question of how knowledge spreads through a population of learners has attracted growing interest at the intersection of applied mathematics, cognitive science, and educational research. While classical models of epidemiology [1] and population dynamics [2] have long exploited reaction-diffusion frameworks to describe the propagation of biological phenomena, their systematic transposition to the realm of learning and knowledge acquisition remains largely underexplored.

The analogy between the spread of an infectious disease and the diffusion of knowledge within a learning community is both natural and instructive. In both cases, a “carrier”—an individual who has acquired a competence—can transmit it to neighbors through interaction, while decay mechanisms (forgetting, disengagement) continuously erode the acquired level. This structural similarity motivates the use of partial differential equations (PDEs) of reaction-diffusion type as a theoretical backbone for modeling knowledge propagation.

Early contributions to the mathematical modeling of learning date back to the work of [3] and [4], who proposed stochastic and linear models of individual learning dynamics. More recently, agent-based models [5] and network diffusion models [6] have been used to study collective knowledge acquisition in social systems. However, these approaches often lack the analytical tractability that continuous PDE models offer, particularly with respect to steady-state analysis and stability theory.

From a pedagogical standpoint, the literature consistently emphasizes the heterogeneous effectiveness of teaching methods. Active learning strategies, collaborative work, and personalized tutoring have been shown to outperform purely lecture-based instruction in terms of long-term retention and competence transfer [7] [8]. Yet, no unified mathematical framework exists that simultaneously accounts for the spatial (social) propagation of knowledge, the temporal dynamics of forgetting, the nonlinear threshold effects in learning, and the modulating role of pedagogical methods.

The present paper aims to fill this gap by proposing a reaction-diffusion PDE model in which:

(i) The diffusion operator captures the social contagion of knowledge, with a diffusion coefficient that is linearly modulated by $n$ pedagogical methods $M_i$ through efficiency parameters ${\delta }_i$ ;

(ii) The reaction term combines a logistic saturation law $u\left(1-u\right)$ , a Hill-type threshold function $h\left(u\right)=u^p/\left({\theta }^p+u^p\right)$ encoding the existence of a triggering threshold $\theta$ below which learning is self-inhibiting, a nonlocal social facilitation integral weighted by a Gaussian kernel $w\left(x,x^{\prime }\right)$ , and a linear forgetting term $-\lambda u$ consistent with the Ebbinghaus forgetting curve [9];

(iii) The local learning rate $r\left(x,t\right)$ integrates the pedagogical efficiency coefficients ${\alpha }_i$ of each method and the inhibitory effect of cognitive stress through a parameter $\beta$ , in the spirit of the Yerkes-Dodson law [10].

The model is formulated on a bounded domain $\Omega \subset {ℝ}^d$ with homogeneous Neumann boundary conditions, ensuring a closed-population interpretation. We establish the mathematical well-posedness of the problem, analyze the existence and local stability of equilibria, and conduct a fictitious numerical simulation under MATLAB to validate the qualitative predictions of the model. An extended two-compartment formulation coupling cognitive competences $u$ and non-cognitive competences $v$ (motivation, self-efficacy) is also introduced and briefly analyzed.

Comparison with existing models. While agent-based [5] and network diffusion [6] models have successfully captured qualitative aspects of knowledge spread, they typically lack the analytical tractability of continuous PDEs for equilibrium and stability analysis. Classical reactiondiffusion models of innovation diffusion (e.g., Bass-type equations) do not incorporate pedagogical modulation, threshold effects, or forgetting. The present framework is, to the best of our knowledge, the first PDE model that simultaneously integrates (i) Pedagogically modulated diffusion, (ii) A Hill-type threshold, (iii) Nonlocal social facilitation, (iv) Ebbinghaus forgetting, and (v) Coupling between cognitive and non-cognitive competences. This combination distinguishes it from previous approaches and provides a mathematically rigorous foundation for educational modeling.

The remainder of the paper is organized as follows. Section 2 presents the full model formulation and the interpretation of each term. Section 3 is devoted to the mathematical analysis: well-posedness, equilibria, and stability. Section 4 describes the fictitious numerical simulation and discusses the results. Section 5 introduces the two-compartment extended model. Section 6 concludes the paper and outlines perspectives for calibration on real educational data.

2. Model Formulation

2.1. General Framework

Let $\Omega \subset {ℝ}^d$ ($d\ge 1$ ) be a bounded, connected, open domain with smooth boundary $\partial \Omega$ , and let $T>0$ be a finite time horizon. We denote by $Q_T=\Omega \times \left(0,T\right)$ the space-time cylinder. The state variable $u:{\overline{Q}}_T\to \left[0,1\right]$ represents the local density of mastered competence at position $x\in \Omega$ and time $t\in \left[0,T\right]$ , where $u=0$ corresponds to total absence of mastery and $u=1$ to complete mastery of the considered competence. The spatial coordinate $x$ is an abstract index that orders learners according to their initial knowledge level (e.g., from lowest to highest baseline competence). This interpretation is used consistently throughout the paper: in particular, it allows the initial condition $u_0\left(x\right)$ to represent a heterogeneous population where learners with higher $x$ are not necessarily spatially contiguous but are simply labeled by their initial proficiency.

Definition 2.1. The function $u\left(x,t\right)$ is called the knowledge density field. For a given population of learners indexed by $x\in \Omega$ , the quantity

$\overline{u}\left(t\right)=\frac{1}{\left|\Omega \right|}{\displaystyle {\int }_{\Omega }u\left(x,t\right)\text{d}x}$(1)

defines the mean competence level of the population at time $t$ .

2.2. The Reaction-Diffusion Equation

The spatiotemporal evolution of $u$ is governed by the following nonlinear reaction-diffusion equation:

$\boxed{\frac{\partial u}{\partial t}=\nabla \cdot \left[D\left(x,t\right)\nabla u\right]+f\left(u,x,t\right),\left(x,t\right)\in Q_T,}$ (2)

supplemented by homogeneous Neumann boundary conditions:

$\frac{\partial u}{\partial n}=0,\left(x,t\right)\in \partial \Omega \times \left(0,T\right),$ (3)

and an initial condition:

$u\left(x,0\right)=u_0\left(x\right),x\in \Omega ,$ (4)

where $u_0\in L^{\infty }\left(\Omega \right)$ with $0\le u_0\left(x\right)\le 1$ almost everywhere.

Remark 2.2. The Neumann boundary condition (3) encodes the assumption of a closed population: there is no flux of knowledge across the boundary of Ω. This is appropriate when Ω represents a self-contained learning community (e.g., a university cohort) with negligible external exchanges.

2.3. The Diffusion Term

The diffusion operator $\nabla \cdot \left[D\left(x,t\right)\nabla u\right]$ models the social contagion of knowledge: a learner who has mastered a competence can transmit it to neighboring learners through peer interaction, collaborative work, or informal exchanges. The position-dependent diffusion coefficient is defined as:

$D\left(x,t\right)=D_0+{\displaystyle \sum }_{i=1}^n{\delta }_i{M}_i\left(x,t\right),$ (5)

where:

• $D_0\ge 0$ is the baseline diffusion coefficient, representing the natural propagation of knowledge in the absence of any structured pedagogical intervention (e.g., informal peer exchanges);

• $M_i\left(x,t\right)\in \left[0,1\right]$ is the normalized intensity of pedagogical method $i$ at position $x$ and time $t$ . Five methods are considered in this paper:

$M_1$ : Lecture-based teaching (cours magistral),

$M_2$ : Active learning/practical work (TD/TP),

$M_3$ : Collaborative learning (group work),

$M_4$ : Digitalre sources (online platforms, videos),

$M_5$ : Personalized tutoring;

• ${\delta }_i\ge 0$ is the diffusion efficiency of method $i$ , quantifying its contribution to the spread of knowledge through social interaction. Methods that promote exchange and peer interaction (e.g., $M_3$ , $M_5$ ) are expected to yield higher values of ${\delta }_i$ .

Why pedagogical methods affect both diffusion and reaction. The coefficients ${\delta }_i$ and ${\alpha }_i$ capture two distinct but complementary roles of a pedagogical method $M_i$ . The diffusion efficiency ${\delta }_i$ quantifies the method’s capacity to promote social contagion of knowledge: methods that encourage peer interaction, discussion, and collaborative problem-solving (e.g., $M_3$ , $M_5$ ) increase the effective diffusivity $D\left(x,t\right)$ , accelerating the spatial spread of competence across the learner population. The pedagogical efficiency ${\alpha }_i$ measures the method’s direct impact on an individual’s local learning rate, independent of peer influence (e.g., through clarity of explanation, structured exercises, or feedback quality). We assume that the effects of different methods are additive and independent (no interaction terms such as ${\delta }_{ij}{M}_i{M}_j$ ), which is a simplifying firstorder approximation. Interaction effects (e.g., synergy between active learning and digital resources) are intentionally excluded from the present model but could be incorporated in future extensions.

Assumption 2.3. There exist constants $D_{\min },D_{\max }>0$ such that

$D_{\text{min}}\le D\left(x,t\right)\le D_{\text{max}},\forall \left(x,t\right)\in {\overline{Q}}_T.$(6)

This is satisfied whenever $D_0>0$ and the $M_i$ are bounded, which we assume throughout.

2.4. The Reaction Term

The reaction term $f\left(u,x,t\right)$ collects all local processes governing knowledge acquisition, social facilitation, and forgetting:

$f\left(u,x,t\right)=r\left(x,t\right)u\left(1-u\right)g\left(u,x,t\right)-\lambda u,$ (7)

where each component is described in the following subsections.

2.4.1. Local Learning Rate $r\left(x,t\right)$

The local learning rate integrates the contributions of pedagogical methods and the inhibitory effect of cognitive stress:

$r\left(x,t\right)=r_0+{\displaystyle \sum }_{i=1}^n{\alpha }_i{M}_i\left(x,t\right)-\beta S\left(x,t\right),$ (8)

where:

• $r_0\ge 0$ is the intrinsic learning rate, modeling spontaneous knowledge acquisition without any pedagogical intervention;

• ${\alpha }_i\ge 0$ is the pedagogical efficiency of method $i$ , measuring its direct contribution to accelerating learning. These parameters are the primary control variables of the model and can be calibrated from standardized assessment data;

• $S\left(x,t\right)\ge 0$ denotes the cognitive stress level perceived by learners at position $x$ and time $t$ . It may encompass time pressure, perceived difficulty, or environmental factors;

• $\beta \ge 0$ is the stress attenuation coefficient. A high value of $\beta$ indicates that stress strongly inhibits learning efficiency, consistently with the Yerkes-Dodson law [10]: moderate stress may be stimulating, but excessive stress is detrimental.

Remark 2.4. We assume $r\left(x,t\right)>0$ for all $\left(x,t\right)\in Q_T$ , which requires that the net contribution of pedagogical methods dominates the stress effect: $r_0+{\displaystyle {\sum }_i{\alpha }_i{M}_i}>\beta S$ . If this condition fails locally, the model predicts a region of learning inhibition, which is pedagogically meaningful.

2.4.2. Logistic Growth Term $u\left(1-u\right)$

The factor $u\left(1-u\right)$ is the classical logistic saturation term:

• When $u\approx 0$ , growth is approximately linear in $u$ (initial exponential phase);

• As $u\to 1$ , growth decelerates and vanishes, reflecting the saturation of competence: one cannot exceed complete mastery.

This term ensures that $u=0$ and $u=1$ are invariant states of the reaction dynamics in the absence of forgetting.

2.4.3. Threshold and Social Facilitation Function $g\left(u,x,t\right)$

The function $g\left(u,x,t\right)$ combines two fundamental mechanisms:

$g\left(u,x,t\right)=h\left(u\right)\cdot F\left(u,x,t\right),$ (9)

where $h\left(u\right)$ is a Hill-type threshold function and $F\left(u,x,t\right)$ is a social facilitation operator.

Hill-type threshold function.

$h\left(u\right)=\frac{u^p}{{\theta }^p+u^p},$ (10)

with:

• $\theta \in \left(0,1\right)$ the critical mastery threshold. Below $\theta$ , learning is self-inhibiting (dropout effect); above $\theta$ , learning becomes self-sustaining (snowball effect);

• $p\ge 1$ the Hill coefficient, controlling the sharpness of the transition. For $p=1$ , the transition is smooth (sigmoid-like); for $p\gg 1$ , it approaches a step function.

Proposition 2.5. The Hill function $h:\left[0,1\right]\to \left[0,1\right]$ defined in (10) satisfies:

(i) $h\left(0\right)=0$ , ${\lim }_{u\to \infty }h\left(u\right)=1$ ;

(ii) $h$ is strictly increasing on $\left[0,1\right]$ ;

(iii) $h\left(u\right)<\frac{1}{2}$ for $u<\theta$ , and $h\left(u\right)>\frac{1}{2}$ for $u>\theta$ ;

(iv) $h^{\prime }\left(u\right)=\frac{p{\theta }^p{u}^{p-1}}{{\left({\theta }^p+u^p\right)}^2}>0$ for all $u>0$ .

Proof. Properties (i)-(iii) follow directly from the definition (10) by elementary algebraic manipulations. Property (iv) is obtained by differentiating (10) with respect to $u$ using the quotient rule. □

Social facilitation operator.

$F\left(u,x,t\right)=1+\gamma {\displaystyle {\int }_{\Omega }w\left(x,x^{\prime }\right)u\left(x^{\prime },t\right)\text{d}{x}^{\prime }},$ (11)

where:

• $\gamma \ge 0$ is the social facilitation intensity. A higher $\gamma$ means that a learner’s progress is more strongly influenced by the competence level of their peers;

• $w\left(x,x^{\prime }\right)\ge 0$ is the interaction kernel, weighting the influence of learners at $x\text{'}$ on those at $x$ . We adopt the standard Gaussian kernel:

$w\left(x,x^{\prime }\right)=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{d/2}}\text{exp}\left(-\frac{{\Vert x-x^{\prime }\Vert }^2}{2{\sigma }^2}\right),$ (12)

where $\sigma >0$ is the social interaction range.

Remark 2.6. The operator $\mathcal{W}\left[u\right]\left(x,t\right)={\displaystyle {\int }_{\Omega }w\left(x,x^{\prime }\right)u\left(x^{\prime },t\right)\text{d}{x}^{\prime }}$ is a convolution-type operator. Its boundedness on $L^2\left(\Omega \right)$ follows from Young’s convolution inequality, since $w\in L^1\left({ℝ}^d\right)$ by virtue of (12).

2.4.4. Forgetting Term $-\lambda u$

The term $-\lambda u$ , with $\lambda \ge 0$ , models the natural decay of knowledge when competences are not regularly reactivated. This linear decay is consistent with the Ebbinghaus forgetting curve [9], which posits that memory retention decays approximately exponentially over time in the absence of reinforcement.

2.5. Complete Model

Assembling all terms, the complete model reads:

$\boxed{\begin{array}{l}\frac{\partial u}{\partial t}=\nabla \cdot \left[\left(D_0+{\displaystyle \sum }_{i=1}^n{\delta }_i{M}_i\right)\nabla u\right]+\left(r_0+{\displaystyle \sum }_{i=1}^n{\alpha }_i{M}_i-\beta S\right)u\left(1-u\right)\\ \cdot \frac{u^p}{{\theta }^p+u^p}\left(1+\gamma {\displaystyle {\int }_{\Omega }w\left(x,x^{\prime }\right)u\left(x^{\prime },t\right)\text{d}{x}^{\prime }}\right)-\lambda u,\end{array}}$ (13)

for $\left(x,t\right)\in Q_T$ , with boundary condition (3) and initial condition (4).

2.6. Parameter Summary

Table 1 summarizes all parameters of the model, their mathematical role, and their pedagogical interpretation.

Table 1. Summary of model parameters and their pedagogical interpretation.

Symbol	Mathematical role	Pedagogical interpretation	Sign
$D_0$	Baseline diffusion coefficient	Natural knowledge propagation (informal exchanges)	≥0
${\delta }_i$	Diffusion efficiency of method $i$	Contribution of method $i$ to social knowledge spread	≥0
$r_0$	Intrinsic learning rate	Spontaneous learning without intervention	≥0
${\alpha }_i$	Pedagogical efficiency of method $i$	Direct impact of method $i$ on learning speed	≥0
$\beta$	Stress attenuation coefficient	Inhibitory effect of cognitive stress (Yerkes-Dodson)	≥0
$\theta$	Critical mastery threshold	Triggering threshold for self-sustained learning	(0,1)
$p$	Hill coefficient	Sharpness of the learning threshold transition	≥1
$\gamma$	Social facilitation intensity	Peer influence on individual learning	≥0
$\sigma$	Social interaction range	Spatial reach of peer influence (Gaussian kernel)	>0
$\lambda$	Forgetting rate	Speed of knowledge decay (Ebbinghaus curve)	≥0

2.7. Extended Two-Compartment Model

To account for the interplay between cognitive and non-cognitive competences, we introduce a second state variable $v\left(x,t\right)\in \left[0,1\right]$ representing non-cognitive skills (motivation, self-efficacy, metacognition). The extended system reads:

$\{\begin{array}{l}\frac{\partial u}{\partial t}=D_u\Delta u+r_u{u}^a{v}^b\left(1-u\right)-{\lambda }_{u}u,\hfill \\ \frac{\partial v}{\partial t}=D_v\Delta v+r_v{u}^c{v}^d\left(1-v\right)-{\lambda }_{v}v,\hfill \end{array}\left(x,t\right)\in Q_T,$ (14)

supplemented by homogeneous Neumann boundary conditions for both $u$ and $v$ . The cross-coupling terms $u^a{v}^b$ and $u^c{v}^d$ model the bidirectional interaction: cognitive competences reinforce motivation, and motivation facilitates knowledge acquisition. The exponents $a,b,c,d\ge 0$ control the strength of each coupling.

Remark 2.7. When $b=d=0$ , the two equations in (14) decouple, and each reduces to an independent reaction-diffusion equation of Fisher-KPP type [11] [12]. The full coupled system (14) is therefore a generalization of the classical Fisher-KPP model to the educational context.

3. Mathematical Analysis

This section establishes the theoretical foundations of the model (13). We successively address the well-posedness of the problem (existence, uniqueness, and positivity of solutions), the existence of equilibrium states, and their local stability.

3.1. Functional Setting and Assumptions

We work in the following functional framework. Let $q\in \left(1,\infty \right)$ and define the Sobolev space $W^{1,q}\left(\Omega \right)$ endowed with its standard norm. For the parabolic problem, we use the Bochner space $L^2\left(0,T;H^1\left(\Omega \right)\right)$ , where $H^1\left(\Omega \right)=W^{1,2}\left(\Omega \right)$ .

We collect the standing assumptions used throughout this section.

Assumption 3.1. The following conditions hold:

(A1) $\Omega \subset {ℝ}^d$ is a bounded, connected, open set with boundary $\partial \Omega$ of class $C^{1,1}$ .

(A2) The diffusion coefficient satisfies (6): $0<D_{\text{min}}\le D\left(x,t\right)\le D_{\text{max}}<\infty$ for all $\left(x,t\right)\in {\overline{Q}}_T$ .

(A3) The pedagogical method intensities satisfy $M_i\in L^{\infty }\left(Q_T\right)$ , $0\le M_i\le 1$ , for all $i=1,\cdots ,n$ .

(A4) The stress field satisfies $S\in L^{\infty }\left(Q_T\right)$ , $S\ge 0$ .

(A5) The initial datum satisfies $u_0\in L^{\infty }\left(\Omega \right)$ , $0\le u_0\left(x\right)\le 1$ almost everywhere in Ω.

(A6) The interaction kernel $w\in L^1\left({ℝ}^d\right)\cap L^{\infty }\left({ℝ}^d\right)$ is nonnegative and symmetric: $w\left(x,x^{\prime }\right)=w\left(x^{\prime },x\right)\ge 0$ .

3.2. Well-Posedness

Theorem 3.2. (Existence and Uniqueness) Under Assumption 3.1, there exists a unique weak solution

$u\in L^2\left(0,T;H^1\left(\Omega \right)\right)\cap C\left(\left[0,T\right];L^2\left(\Omega \right)\right)\cap L^{\infty }\left(Q_T\right)$ (15)

to problem (13)-(3)-(4), in the sense that for all $\phi \in H^1\left(\Omega \right)$ and a.e. $t\in \left(0,T\right)$ :

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{\Omega }u\phi \text{d}x}+{\displaystyle {\int }_{\Omega }D\left(x,t\right)\nabla u\cdot \nabla \phi \text{d}x}={\displaystyle {\int }_{\Omega }f\left(u,x,t\right)\phi \text{d}x},$ (16)

with $u\left(\cdot ,0\right)=u_0$ in $L^2\left(\Omega \right)$ .

Proof. The proof proceeds in three steps.

Step 1: Galerkin approximation. Let ${\left\{e_k\right\}}_{k\ge 1}$ be the orthonormal basis of $L^2\left(\Omega \right)$ formed by the eigenfunctions of $-\text{Δ}$ with homogeneous Neumann boundary conditions. For each $m\ge 1$ , we seek an approximate solution of the form $u_m\left(x,t\right)={\displaystyle {\sum }_{k=1}^m{c}_k^m\left(t\right)e_k\left(x\right)}$ , where the coefficients $c_k^m$ satisfy the finite-dimensional ODE system obtained by projecting (16) onto $\text{span}\left\{e_1,\cdots ,e_m\right\}$ . By the Cauchy-Lipschitz theorem, this system admits a unique local solution on some interval $\left[0,T_m\right)$ .

Step 2: A priori estimates. Testing (16) against $u_m$ and using Assumption 3.1(A2), we obtain:

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert u_m\Vert }_{L^2\left(\Omega \right)}^2+D_{\min }{\Vert \nabla u_m\Vert }_{L^2\left(\Omega \right)}^2\le {\displaystyle {\int }_{\Omega }f\left(u_m,x,t\right)u_m\text{d}x}.$ (17)

Since $f$ is locally Lipschitz and $\left|f\left(u,x,t\right)\right|\le C\left(1+\left|u\right|\right)$ for a constant $C>0$ independent of $m$ , Gronwall’s lemma yields uniform bounds:

$\underset{t\in \left[0,T\right]}{\sup }{\Vert u_m\left(t\right)\Vert }_{L^2\left(\Omega \right)}^2+{\displaystyle {\int }_0^T{\Vert \nabla u_m\Vert }_{L^2\left(\Omega \right)}^2\text{d}t}\le C_T,$ (18)

for a constant $C_T$ independent of $m$ . This extends the local solution to the full interval $\left[0,T\right]$ .

Step 3: Passage to the limit. The uniform bounds (18) allow extraction of a subsequence (not relabeled) such that $u_m\rightharpoonup u$ weakly in $L^2\left(0,T;H^1\left(\Omega \right)\right)$ and strongly in $L^2\left(Q_T\right)$ by the Aubin-Lions compactness lemma. The nonlinear terms converge by dominated convergence. The uniqueness follows from a standard energy estimate applied to the difference of two solutions. □

Theorem 3.3 (Positivity and $L^{\infty }$ Bound). Under Assumption 3.1, if $0\le u_0\left(x\right)\le 1$ a.e., then the solution $u$ of (13)-(4) satisfies:

$0\le u\left(x,t\right)\le 1,a.e.in{Q}_T.$ (19)

Proof. Lower bound. Testing (16) against $u^{-}=\text{max}\left(-u,0\right)\ge 0$ (the negative part of $u$ ) and using $f\left(0,x,t\right)=0$ , a standard comparison argument shows that $u^{-}\left(x,t\right)=0$ a.e., hence $u\ge 0$ .

Upper bound. Let $w=u-1$ and test the equation for $w$ against $w^{+}=\text{max}\left(w,0\right)$ . Since $f\left(1,x,t\right)=-\lambda \le 0$ and the reaction term satisfies $f\left(u,x,t\right)\le 0$ for $u\ge 1$ (the logistic factor $u\left(1-u\right)$ is nonpositive for $u\ge 1$ ), the same comparison argument yields $u\le 1$ a.e. □

Remark 3.4. Theorem 3.3 confirms that the closed interval $\left[0,1\right]$ is a positively invariant region for the dynamics (13), which is essential for the biological (pedagogical) consistency of the model: the knowledge density remains bounded between total ignorance and complete mastery at all times.

3.3. Equilibrium States

We now analyze the spatially homogeneous steady states of (13), obtained by setting $\partial u/\partial t=0$ and $\nabla u=0$ (uniform solutions). In this case, the diffusion term vanishes and the equilibria $u^{\text{*}}\in \left[0,1\right]$ satisfy:

$F\left(u^{\text{*}}\right):=r^{\text{*}}{u}^{\text{*}}\left(1-u^{\text{*}}\right)\frac{{\left(u^{\text{*}}\right)}^p}{{\theta }^p+{\left(u^{\text{*}}\right)}^p}\left(1+\gamma \overline{w}{u}^{\text{*}}\right)-\lambda u^{\text{*}}=0,$(20)

where $r^{\text{*}}=r_0+{\displaystyle {\sum }_i{\alpha }_i{\overline{M}}_i}-\beta \overline{S}$ is the spatially averaged learning rate, ${\overline{M}}_i$ and $\overline{S}$ are the spatial means of $M_i$ and $S$ , and $\overline{w}={\displaystyle {\int }_{\Omega }w\left(x,x^{\prime }\right)\text{d}{x}^{\prime }}$ (constant by normalization of the Gaussian kernel).

Proposition 3.5 (Trivial Equilibrium). $u^{*}=0$ is always an equilibrium of (20).

Proof. Direct substitution of $u^{\text{*}}=0$ into (20) gives $F\left(0\right)=0$ . □

Proposition 3.6 (Nontrivial Equilibria). Assume $\theta \in \left(0,1\right)$ and that the net learning rate satisfies

$r^{\text{*}}>\frac{2\lambda }{\left(1-\theta \right)\left(1+\gamma \overline{w}\theta \right)}.$(21)

Then there exists at least one nontrivial equilibrium $u^{\text{**}}\in \left(\theta ,1\right)$ satisfying (20).

Proof. Define $\Phi \left(u\right)=F\left(u\right)/u$ for $u\in \left(0,1\right]$ . Then:

$\Phi \left(u\right)=r^{*}\left(1-u\right)\frac{u^p}{{\theta }^p+u^p}\left(1+\gamma \overline{w}u\right)-\lambda .$

At $u={\theta }^{+}$ : since $h\left(\theta \right)=1/2$ , we have $\Phi \left(\theta \right)=r^{*}\left(1-\theta \right)/2\cdot \left(1+\gamma \overline{w}\theta \right)-\lambda$ . Condition (21) guarantees $\Phi \left(\theta \right)>0$ . At $u=1^{-}$ : $\Phi \left(1\right)=-\lambda <0$ . By the intermediate value theorem, there exists $u^{\text{**}}\in \left(\theta ,1\right)$ such that $\Phi \left(u^{**}\right)=0$ , i.e., $F\left(u^{\text{**}}\right)=0$ . □

Remark 3.7. The condition $r^{*}>\lambda$ has a clear pedagogical meaning: the net learning rate (boosted by pedagogical methods) must exceed the forgetting rate for nontrivial knowledge equilibria to exist. If $r^{*}\le \lambda$ , the only equilibrium is $u^{*}=0$ (total knowledge loss), illustrating the critical role of pedagogical intervention in sustaining long-term competence.

3.4. Local Stability Analysis

We analyze the local stability of the equilibria identified in Section 3.3 via linearization.

3.4.1. Linearization around $u^{\ast }=0$

Linearizing (13) around $u^{\text{*}}=0$ (with spatially homogeneous perturbation $u=\epsilon \varphi \left(x\right){\text{e}}^{\mu t}$ , $\epsilon \ll 1$ ), the leading-order eigenvalue problem is:

$\mu \varphi =D^{\text{*}}\Delta \varphi -\lambda \varphi ,x\in \Omega ,\frac{\partial \varphi }{\partial n}=0\text{on}\partial \Omega ,$ (22)

where $D^{\text{*}}=D_0+{\displaystyle {\sum }_i{\delta }_i{\overline{M}}_i}$. The eigenvalues are:

${\mu }_k=-D^{\text{*}}{\nu }_k-\lambda ,k=0,1,2,\cdots ,$ (23)

where $0={\nu }_0<{\nu }_1\le {\nu }_2\le \cdots$ are the Neumann eigenvalues of $-\text{Δ}$ on Ω.

Proposition 3.8 (Stability of the Trivial Equilibrium). The trivial equilibrium $u^{*}=0$ is always locally asymptotically stable with respect to the linearized dynamics (22).

Proof. Since ${\nu }_k\ge 0$ and $\lambda \ge 0$ , we have ${\mu }_k=-D^{\text{*}}{\nu }_k-\lambda \le 0$ for all $k\ge 0$ . Moreover, the reaction term $f\left(u,x,t\right)$ vanishes to higher order at $u=0$ due to the Hill factor $h\left(u\right)=u^p/\left({\theta }^p+u^p\right)~u^p/{\theta }^p\to 0$ as $u\to 0$ (for $p\ge 1$ ). Hence, the linearization at $u^{\text{*}}=0$ gives ${\mu }_k=-\lambda <0$ for $k=0$ , confirming local asymptotic stability. □

Remark 3.9. The stability of $u^{*}=0$ reflects the dropout effect: a population with very low initial knowledge density (below $\theta$ ) tends to converge back to zero mastery in the absence of sufficient pedagogical stimulation. This is consistent with the Hill threshold mechanism.

3.4.2. Linearization around the Nontrivial Equilibrium $u^{\ast \ast }$

Let $u^{\text{**}}\in \left(\theta ,1\right)$ be a nontrivial equilibrium given by Proposition 3.6. Setting $u=u^{\text{**}}+\epsilon v$ and linearizing (13), the perturbation $v$ satisfies:

$\frac{\partial v}{\partial t}=D^{*}\Delta v+a^{*}v+b^{*}{\displaystyle {\int }_{\Omega }w\left(x,x^{\prime }\right)v\left(x^{\prime },t\right)\text{d}{x}^{\prime }},$ (24)

where:

$\begin{array}{c}{a}^{\text{*}}={\frac{\partial f}{\partial u}|}_{u=u^{\text{**}}}\\ =r^{\text{*}}\left(1-2u^{\text{**}}\right)h\left(u^{\text{**}}\right)\left(1+\gamma \overline{w}{u}^{\text{**}}\right)\\ +r^{\text{*}}\left(u^{\text{**}}\right)\left(1-u^{\text{**}}\right)h^{\prime }\left(u^{\text{**}}\right)\left(1+\gamma \overline{w}{u}^{\text{**}}\right)-\lambda ,\end{array}$(25)

$b^{\text{*}}=r^{\text{*}}{u}^{\text{**}}\left(1-u^{\text{**}}\right)h\left(u^{\text{**}}\right)\gamma .$ (26)

Theorem 3.10 (Stability Criterion for $u^{**}$ ). The nontrivial equilibrium $u^{**}$ is locally asymptotically stable if and only if:

$a^{*}+b^{*}\overline{w}\left|\Omega \right|<D^{*}{\nu }_1,$ (27)

where ${\nu }_1>0$ is the first nonzero Neumann eigenvalue of $-\text{Δ}$ on Ω.

Proof. Seeking solutions of (24) of the form $v\left(x,t\right)=\varphi \left(x\right){\text{e}}^{\mu t}$ , we obtain for each Neumann eigenmode $\varphi =e_k$ :

${\mu }_k=-D^{\text{*}}{\nu }_k+a^{\text{*}}+b^{\text{*}}{\widehat{w}}_k,$(28)

where ${\widehat{w}}_k={\displaystyle {\int }_{\Omega }w\left(x,x^{\prime }\right)e_k\left(x^{\prime }\right)\text{d}{x}^{\prime }}$ is the action of the kernel on the $k$ -th eigenfunction. For stability, we require ${\mu }_k<0$ for all $k\ge 0$ . The most restrictive condition corresponds to $k=0$ (spatially uniform mode, ${\nu }_0=0$ , ${\widehat{w}}_0=\overline{w}\left|\Omega \right|$ ), giving ${\mu }_0=a^{*}+b^{*}\overline{w}\left|\Omega \right|<0$ , which is (27). □

Corollary 3.11. If $u^{**}$ is close to 1 (near-complete mastery), then $1-2u^{**}<0$ and the term $a^{*}$ is dominated by the negative contribution $r^{*}\left(1-2u^{**}\right)h\left(u^{**}\right)(\cdots )<0$ , making the stability condition (27) more easily satisfied. Pedagogically, this means that a population that has collectively reached high mastery is self-stabilizing.

3.5. Summary of Theoretical Results

Table 2 summarizes the main theoretical results established in this section.

Table 2. Summary of theoretical results for model (13).

Result	Condition	Pedagogical meaning
Existence & uniqueness	Assumptions 3.1	Model is mathematically consistent
Positivity & $L^{\infty }$ bound	$0\le u_0\le 1$	Knowledge density stays in $\left[0,1\right]$
Trivial equilibrium $u^{\text{*}}=0$	Always exists	Total knowledge loss is always possible
Nontrivial equilibrium $u^{\text{**}}$	$r^{\text{*}}>\lambda$	Sustained learning requires $r^{\text{*}}>\lambda$
Stability of $u^{\text{*}}=0$	Always (Hill threshold)	Dropout is stable below threshold $\theta$
Stability of $u^{\text{**}}$	(27)	High-mastery state is self-stabilizing

4. Numerical Simulation

4.1. Objectives and Strategy

The purpose of this section is twofold: (i) to illustrate the qualitative behavior of model (13) in a set of illustrative scenarios, and (ii) to compare the impact of the five pedagogical methods under a plausible but fictive parameter configuration. Because no empirical dataset is used, this simulation should be viewed as a proofofconcept and a demonstration of the model’s diagnostic capabilities, not as a validation against real data. All parameter values (Table 3) are chosen to be pedagogically plausible and consistent with the qualitative findings of the educational literature [7] [8]. Importantly, the resulting ranking of methods (Section 4.5.4) is conditional on the chosen parameter set; a different calibration (e.g., different ${\alpha }_i$ , ${\delta }_i$ , or $\theta$ ) could modify the ordering. The simulation thus serves as a flexible tool for “whatif” analyses rather than a definitive assessment.

The simulation is conducted on the one-dimensional spatial domain $\Omega =\left[0,1\right]$ , discretized with $N_x=100$ uniform grid points, over the time interval $\left[0,T\right]$ with $T=50$ (arbitrary time units representing weeks). The spatial variable $x\in \left[0,1\right]$ indexes learners ordered by their initial knowledge level. All computations are performed in MATLAB using finite differences for space discretization and the implicit-explicit (IMEX) scheme for time integration, which treats the diffusion term implicitly (for stability) and the reaction term explicitly (for computational efficiency).

4.2. Parameter Values

The parameter values used in the simulation are summarized in Table 3. Five scenarios are considered, each activating a single pedagogical method at full intensity ($M_i\equiv 1$ , all other $M_j\equiv 0$ , $j\ne i$ ) to isolate the individual contribution of each method.

Table 3. Parameter values used in the numerical simulation.

Parameter	Description	Value
$D_0$	Baseline diffusion coefficient	0.02
${\delta }_1$	Diffusion efficiency-Lecture	0.01
${\delta }_2$	Diffusion efficiency-Active learning	0.04
${\delta }_3$	Diffusion efficiency-Collaborative	0.08
${\delta }_4$	Diffusion efficiency-Digital	0.03
${\delta }_5$	Diffusion efficiency-Tutoring	0.06
$r_0$	Intrinsic learning rate	0.10
${\alpha }_1$	Pedagogical efficiency-Lecture	0.20
${\alpha }_2$	Pedagogical efficiency-Active	0.50
${\alpha }_3$	Pedagogical efficiency-Collaborative	0.45
${\alpha }_4$	Pedagogical efficiency-Digital	0.30
${\alpha }_5$	Pedagogical efficiency-Tutoring	0.60
$\beta$	Stress attenuation coefficient	0.15
$\overline{S}$	Mean stress level	0.30
$\theta$	Critical mastery threshold	0.25
$p$	Hill coefficient	3
$\gamma$	Social facilitation intensity	0.40
$\sigma$	Social interaction range	0.15
$\lambda$	Forgetting rate	0.05

Remark 4.1. The ranking ${\alpha }_5>{\alpha }_2>{\alpha }_3>{\alpha }_4>{\alpha }_1$ reflects the meta-analytic findings of [8]: personalized tutoring has the highest effect size on student achievement, followed by active learning and collaborative work, while lecture-based teaching alone yields the lowest individual impact.

4.3. Numerical Scheme

Let $x_j=\left(j-1\right)\text{Δ}x$ for $j=1,\cdots ,N_x$ with $\text{Δ}x=1/\left(N_x-1\right)$ , and $t^n=n\text{Δ}t$ with time step $\text{Δ}t=0.05$ . Let $u_j^n\approx u\left(x_j,t^n\right)$ .

The nonlocal social facilitation term is approximated by the trapezoidal rule:

$\mathcal{W}{\left[u^n\right]}_j=\Delta x{\displaystyle \sum }_{k=1}^{N_x}w\left(x_j,x_k\right)u_k^n,w\left(x_j,x_k\right)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{{\left(x_j-x_k\right)}^2}{2{\sigma }^2}\right).$ (29)

The full IMEX scheme reads: given $u^n$ , find $u^{n+1}$ satisfying

$\begin{array}{c}\frac{u_j^{n+1}-u_j^n}{\Delta t}=\frac{1}{\Delta x^2}\left[D_{j+\frac{1}{2}}^n\left(u_{j+1}^{n+1}-u_j^{n+1}\right)-D_{j-\frac{1}{2}}^n\left(u_j^{n+1}-u_{j-1}^{n+1}\right)\right]\\ +r_j^n{u}_j^n\left(1-u_j^n\right)h\left(u_j^n\right)\left(1+\gamma \mathcal{W}{\left[u^n\right]}_j\right)-\lambda u_j^n,\end{array}$ (30)

for $j=2,\cdots ,N_x-1$ , with Neumann boundary conditions $u_1^{n+1}=u_2^{n+1}$ and $u_{N_x}^{n+1}=u_{N_x-1}^{n+1}$ , and where $D_{j\pm \frac{1}{2}}^n=\left(D_j^n+D_{j\pm 1}^n\right)/2$ . The implicit diffusion part yields a tridiagonal linear system solved at each time step by the Thomas algorithm.

Remark 4.2. The IMEX scheme is unconditionally stable with respect to the diffusion term (implicit treatment) and first-order accurate in time. The CFL condition for the explicit reaction term is automatically satisfied since the reaction term is bounded on $\left[0,1\right]$ .

4.4. Initial Condition

To represent a realistic starting scenario, we initialize with a small localized group of learners having partial mastery:

$u_0\left(x\right)=0.15+0.20\exp \left(-\frac{{\left(x-0.3\right)}^2}{2\times {0.05}^2}\right),x\in \left[0,1\right].$ (31)

This corresponds to a baseline mastery of 15% across the population, with a localized peak of 35% centered at $x=0.3$ (a small subgroup of more advanced learners). Note that $u_0\left(x\right)\in \left[0,1\right]$ for all $x$ , consistent with Theorem 3.3.

4.5. Simulation Results and Discussion

The numerical simulation of model (13) over $T=50$ weeks produces the three-panel figure (Figure 1) described below. All results are consistent with the theoretical predictions of Section 3.

4.5.1. Final Spatial Knowledge Distribution

The left panel of Figure 1 displays the final spatial distribution $u\left(x,T\right)$ for each of the five pedagogical scenarios, alongside the initial condition $u_0\left(x\right)$ (31) (black dashed curve). Several observations are in order.

First, the initial condition exhibits the expected Gaussian peak centered at $x=0.3$ with a maximum of approximately 0.35, superimposed on a baseline of 0.15, confirming that the simulation starts from a realistic heterogeneous configuration with a small subgroup of more advanced learners.

Second, at $t=T=50$ weeks, all five methods have produced a dramatic and spatially uniform elevation of the knowledge density, with all profiles converging to nearly flat distributions well above 0.80. This spatial homogenization is a direct consequence of the diffusion operator $\nabla \cdot \left[D\left(x,t\right)\nabla u\right]$ : the knowledge initially concentrated near $x=0.3$ has diffused outward across the entire learner population. The Hill threshold mechanism (Proposition 2.5) plays a key role here—once the localized peak exceeds $\theta =0.25$ , the snowball effect propagates the learning front across the domain.

Third, a clear ordering is visible between the five final profiles, consistent with the ranking of pedagogical efficiencies ${\alpha }_5>{\alpha }_2>{\alpha }_3>{\alpha }_4>{\alpha }_1$ .

4.5.2. Mean Knowledge Dynamics

The central panel of Figure 1 shows the temporal evolution of the mean knowledge level $\overline{u}\left(t\right)$ for each method. All five curves exhibit the characteristic sigmoidal shape predicted by the logistic-Hill reaction term: an initial slow phase (below the threshold $\theta =0.25$ ), followed by a rapid acceleration once the threshold is crossed, and a final plateau corresponding to the nontrivial equilibrium $u^{\text{**}}$ .

The vertical dashed line at $t=1/\lambda =20$ weeks marks the Ebbinghaus reference time: without pedagogical reinforcement, a knowledge item would decay to $e^{-1}\approx 37\text{\%}$ of its initial value by this time. The fact that all five curves have already surpassed $\overline{u}=0.80$ before $t=20$ weeks demonstrates that the pedagogical methods effectively counteract the forgetting mechanism, consistently with the condition $r^{\text{*}}>\lambda$ established in Proposition 3.6.

4.5.3. Comparative Impact of Pedagogical Methods

Table 4 summarizes the mean knowledge level $\overline{u}\left(T\right)$ at $t=T=50$ weeks for each pedagogical scenario, as well as the half-mastery time $t_{1/2}$ (first time $\overline{u}\left(t\right)\ge 0.50$ ), read directly from Figure 1.

Table 4. Simulation results: mean knowledge at $t=T=50$ weeks and half-mastery time $t_{1/2}$ for each pedagogical method.

Method	${\alpha }_i$	$\overline{u}\left(T\right)$	$t_{1/2}$ (weeks)
$M_5$ : Personalized tutoring	0.60	≈0.940	≈4
$M_2$ : Active learning	0.50	≈0.930	≈5
$M_3$ : Collaborative learning	0.45	≈0.925	≈6
$M_4$ : Digital resources	0.30	≈0.895	≈9
$M_1$ : Lecture-based teaching	0.20	≈0.850	≈18

Figure 1. Numerical simulation of model (13) over $T=50$ weeks for five pedagogical scenarios (one method active at full intensity at a time). Left: final spatial knowledge distribution $u\left(x,T\right)$ ; the black dashed curve shows the initial condition $u_0\left(x\right)$ . Center: temporal evolution of the mean knowledge level $\overline{u}\left(t\right)$ ; the vertical dashed line marks the Ebbinghaus reference time $1/\lambda =20$ weeks. Right: bar chart of $\overline{u}\left(T\right)$ ranked in decreasing order. Parameter values as in Table 3.

The results confirm the expected ranking in full agreement with the meta-analytic findings of [8]. Personalized tutoring ($M_5$ ) achieves the highest long-term mastery and the fastest convergence, with the population reaching half-mastery in approximately 4 weeks. By contrast, lecture-based teaching ($M_1$ ) yields the lowest equilibrium level ($\overline{u}\left(T\right)\approx 0.850$ ) and the slowest convergence ($t_{1/2}\approx 18$ weeks)—more than four times slower than tutoring. This gap of approximately 9% in final mastery between $M_5$ and $M_1$ , while seemingly modest in absolute terms, represents a substantial difference in long-term competence accumulation over an academic year.

4.5.4. Ranking of Pedagogical Methods

The right panel of Figure 1 presents a bar chart of $\overline{u}\left(T\right)$ ranked in decreasing order. The ranking $M_5\succ M_2\succ M_3\succ M_4\succ M_1$ is unambiguous and robust. Three clusters emerge naturally:

(i) High-impact methods ($\overline{u}\left(T\right)>0.92$ ): tutoring ($M_5$ ), active learning ($M_2$ ), and collaborative learning ($M_3$ ). These methods share a high diffusion efficiency (${\delta }_3=0.08$ , ${\delta }_5=0.06$ ) that amplifies the social contagion of knowledge and accelerates the spread of learning across the population.

(ii) Intermediate method ($0.88<\overline{u}\left(T\right)<0.92$ ): digital resources ($M_4$ ), which combine moderate pedagogical efficiency (${\alpha }_4=0.30$ ) with limited social diffusion (${\delta }_4=0.03$ ).

(iii) Low-impact method ($\overline{u}\left(T\right)\approx 0.85$ ): lecture-based teaching ($M_1$ ), characterized by the lowest values of both ${\alpha }_1=0.20$ and ${\delta }_1=0.01$ , reflecting the predominantly passive and non-interactive nature of this format.

Remark 4.3. The simulation confirms all four theoretical predictions of Section 3: (i) positivity and boundedness $u\in \left[0,1\right]$ ; (ii) convergence to nontrivial equilibria $u^{**}>0$ for all methods since $r^{*}>\lambda$ ; (iii) snowball effect triggered by the Hill threshold from the localized initial peak; (iv) Ebbinghaus forgetting visible in the slower convergence of $M_1$ .

5. Study of the Extended Two-Compartment Model

5.1. Motivation

The single-compartment model (13) captures the spatiotemporal dynamics of cognitive competences—declarative and procedural knowledge that can be directly assessed through standardized evaluations. However, a growing body of educational research emphasizes that learning outcomes are also strongly shaped by non-cognitive factors such as motivation, self-efficacy, perseverance, and metacognitive awareness [13]-[15]. These non-cognitive competences interact bidirectionally with cognitive ones: a learner who masters a concept gains confidence (motivation reinforced by cognitive success), while a motivated learner invests more effort and thus acquires knowledge more efficiently (motivation facilitates cognitive acquisition).

To capture this bidirectional coupling, we introduce a second state variable $v\left(x,t\right)\in \left[0,1\right]$ representing the local density of non-cognitive competences (motivation, self-efficacy, engagement) at position $x$ and time $t$ . The resulting two-compartment system generalizes (13) and constitutes the extended model analyzed in this section.

5.2. The Two-Compartment System

The extended model is:

$\boxed{\{\begin{array}{l}\frac{\partial u}{\partial t}=D_u\text{Δ}u+r_u{u}^a{v}^b\left(1-u\right)-{\lambda }_{u}u,\hfill \\ \frac{\partial v}{\partial t}=D_v\text{Δ}v+r_v{u}^c{v}^d\left(1-v\right)-{\lambda }_{v}v,\hfill \end{array}}\left(x,t\right)\in Q_T,$ (32)

supplemented by homogeneous Neumann boundary conditions for both $u$ and $v$ :

$\frac{\partial u}{\partial n}=\frac{\partial v}{\partial n}=0,\left(x,t\right)\in \partial \Omega \times \left(0,T\right),$ (33)

and initial conditions:

$u\left(x,0\right)=u_0\left(x\right),v\left(x,0\right)=v_0\left(x\right),x\in \Omega ,$ (34)

where $0\le u_0,v_0\le 1$ almost everywhere.

The variables and parameters are interpreted as follows:

• $u\left(x,t\right)$ : cognitive competence density (declarative and procedural knowledge);

• $v\left(x,t\right)$ : non-cognitive competence density (motivation, self-efficacy, metacognition);

• $D_u,D_v\ge 0$ : diffusion coefficients for cognitive and non-cognitive competences, respectively. In general, $D_u\ne D_v$ , as the social propagation mechanisms differ between knowledge (transmitted through instruction and peer explanation) and motivation (transmitted through social modeling and peer encouragement);

• $r_u,r_v>0$ : intrinsic growth rates for cognitive and non-cognitive competences;

• $u^a{v}^b$ : cross-coupling term in the $u$ -equation, modeling the facilitation of cognitive acquisition by non-cognitive competences. A motivated learner ($v$ high) acquires knowledge ($u$ ) more efficiently;

• $u^c{v}^d$ : cross-coupling term in the $v$ -equation, modeling the reinforcement of motivation by cognitive success. A knowledgeable learner ($u$ high) gains confidence ($v$ ) more readily;

• ${\lambda }_u,{\lambda }_v\ge 0$ : forgetting/disengagement rates for cognitive and non-cognitive competences, respectively. In general, ${\lambda }_v<{\lambda }_u$ , as motivational states tend to decay more slowly than specific knowledge items;

• $a,b,c,d\ge 0$ : coupling exponents controlling the strength and nonlinearity of the bidirectional interaction.

Remark 5.1. When $b=0$ and $d=0$ , the two equations in (32) decouple: each reduces to an independent Fisher-KPP-type reaction-diffusion equation [11] [12]. The full coupled system (32) is therefore a strict generalization of the classical Fisher-KPP model, with cross-species interaction terms analogous to those found in predator-prey or mutualism models in mathematical ecology [2].

5.3. Well-Posedness of the Extended System

Theorem 5.2 (Well-Posedness of the Extended System). Assume $D_u,D_v>0$ , $r_u,r_v>0$ , ${\lambda }_u,{\lambda }_v\ge 0$ , $a,b,c,d\ge 0$ , and let $u_0,v_0\in L^{\infty }\left(\Omega \right)$ with $0\le u_0,v_0\le 1$ a.e. Then there exists a unique weak solution $\left(u,v\right)$ of (32)-(34) satisfying:

$u,v\in L^2\left(0,T;H^1\left(\Omega \right)\right)\cap C\left(\left[0,T\right];L^2\left(\Omega \right)\right)\cap L^{\infty }\left(Q_T\right),$ (35)

and the invariance property:

$0\le u\left(x,t\right)\le 1,0\le v\left(x,t\right)\le 1,a.e.in{Q}_T.$ (36)

Proof. The proof follows the same Galerkin approximation strategy as Theorem 3.2, applied to the coupled system. The key observation is that the reaction terms $F_u\left(u,v\right)=r_u{u}^a{v}^b\left(1-u\right)-{\lambda }_{u}u$ and $F_v\left(u,v\right)=r_v{u}^c{v}^d\left(1-v\right)-{\lambda }_{v}v$ are locally Lipschitz on ${\left[0,1\right]}^2$ and satisfy the quasi-positivity conditions: $F_u\left(0,v\right)=0\ge 0$ and $F_v\left(u,0\right)=-{\lambda }_v\cdot 0=0\ge 0$ , ensuring nonnegativity of solutions by the comparison principle. The upper bound $u,v\le 1$ follows from the fact that $F_u\left(1,v\right)=-{\lambda }_u\le 0$ and $F_v\left(u,1\right)=-{\lambda }_v\le 0$ , so that $u=1$ and $v=1$ are super-solutions of the respective equations. □

5.4. Equilibrium Analysis of the Extended System

Spatially homogeneous steady states $\left(u^{\text{*}},v^{\text{*}}\right)\in {\left[0,1\right]}^2$ of (32) satisfy the algebraic system:

$\{\begin{array}{l}{r}_u{\left(u^{\text{*}}\right)}^a{\left(v^{\text{*}}\right)}^b\left(1-u^{\text{*}}\right)={\lambda }_u{u}^{\text{*}},\hfill \\ r_v{\left(u^{\text{*}}\right)}^c{\left(v^{\text{*}}\right)}^d\left(1-v^{\text{*}}\right)={\lambda }_v{v}^{\text{*}}.\hfill \end{array}$ (37)

Proposition 5.3 (Equilibria of the Extended System). System (37) always admits the trivial equilibrium $\left(u^{*},v^{*}\right)=\left(0,0\right)$ . Under the conditions

$\frac{r_u}{{\lambda }_u}>1$ and $\frac{r_v}{{\lambda }_v}>1$ ,(38)

there exists at least one nontrivial equilibrium $\left(u^{\text{**}},v^{\text{**}}\right)\in {\left(0,1\right)}^2$ .

Proof. The trivial equilibrium $\left(0,0\right)$ satisfies (37) trivially. For the nontrivial equilibrium, define:

$\Phi \left(u,v\right)=r_u{u}^{a-1}{v}^b\left(1-u\right)-{\lambda }_u,\Psi \left(u,v\right)=r_v{u}^c{v}^{d-1}\left(1-v\right)-{\lambda }_v.$

Under condition (38), one can show by a fixed-point argument (Schauder theorem applied to the map $\left(u,v\right)\mapsto \left(G_u\left(v\right),G_v\left(u\right)\right)$ where $G_u$ and $G_v$ are the respective solution maps) that there exists a fixed point $\left(u^{\text{**}},v^{\text{**}}\right)\in {\left(0,1\right)}^2$ . □

Remark 5.4. Condition (38) has a natural pedagogical interpretation: the ratio $r_u/{\lambda }_u$ (resp. $r_v/{\lambda }_v$ ) must exceed 1 for cognitive (resp. non-cognitive) competences to be sustained against forgetting (resp. disengagement). This generalizes the condition $r^{*}>\lambda$ established for the single-compartment model.

5.5. Local Stability of the Extended System

Let $\left(u^{\text{**}},v^{\text{**}}\right)$ be a nontrivial equilibrium of (37). Setting $u=u^{\text{**}}+\epsilon \tilde{u}$ and $v=v^{\text{**}}+\epsilon \tilde{v}$ and linearizing (32), the perturbation $\left(\tilde{u},\tilde{v}\right)$ satisfies:

$\frac{\partial }{\partial t}\left(\begin{array}{l}\tilde{u}\hfill \\ \tilde{v}\hfill \end{array}\right)=\left(\begin{array}{cc}{D}_u\Delta & 0\\ 0& D_v\Delta \end{array}\right)\left(\begin{array}{l}\tilde{u}\hfill \\ \tilde{v}\hfill \end{array}\right)+J^{\text{*}}\left(\begin{array}{l}\tilde{u}\hfill \\ \tilde{v}\hfill \end{array}\right),$(39)

where $J^{\text{*}}$ is the Jacobian matrix of the reaction terms evaluated at $\left(u^{\text{**}},v^{\text{**}}\right)$ :

$J^{\text{*}}=\left(\begin{array}{ll}{J}_{11}\hfill & J_{12}\hfill \\ J_{21}\hfill & J_{22}\hfill \end{array}\right),$ (40)

with:

$J_{11}=r_u{\left(v^{\text{**}}\right)}^b\left[a{\left(u^{\text{**}}\right)}^{a-1}\left(1-u^{\text{**}}\right)-{\left(u^{\text{**}}\right)}^a\right]-{\lambda }_u,$ (41)

$J_{12}=r_{u}b{\left(u^{\text{**}}\right)}^a{\left(v^{\text{**}}\right)}^{b-1}\left(1-u^{\text{**}}\right),$ (42)

$J_{21}=r_{v}c{\left(u^{\text{**}}\right)}^{c-1}{\left(v^{\text{**}}\right)}^d\left(1-v^{\text{**}}\right),$ (43)

$J_{22}=r_v{\left(u^{\text{**}}\right)}^c\left[d{\left(v^{\text{**}}\right)}^{d-1}\left(1-v^{\text{**}}\right)-{\left(v^{\text{**}}\right)}^d\right]-{\lambda }_v.$ (44)

Theorem 5.5 (Stability Criterion for the Extended System). The nontrivial equilibrium $\left(u^{**},v^{**}\right)$ of (32) is locally asymptotically stable if and only if, for each Neumann eigenvalue ${\nu }_k$ ($k\ge 0$ ), the matrix:

$J^{\text{*}}-\left(\begin{array}{cc}{D}_u{\nu }_k& 0\\ 0& D_v{\nu }_k\end{array}\right)$ (45)

has all eigenvalues with strictly negative real part, i.e.:

$\text{tr}\left(J^{\text{*}}\right)-\left(D_u+D_v\right){\nu }_k<0,\text{det}\left(J^{\text{*}}\right)-\left(D_u{J}_{22}+D_v{J}_{11}\right){\nu }_k+D_u{D}_v{\nu }_k^2>0,$ (46)

for all $k\ge 0$ .

Proof. Seeking solutions of (39) of the form $\left(\tilde{u},\tilde{v}\right)=\left({\varphi }_k,{\psi }_k\right){\text{e}}^{\mu t}$, where $\left({\varphi }_k,{\psi }_k\right)$ are Neumann eigenfunctions, the eigenvalues $\mu$ satisfy:

$\text{det}\left(J^{\text{*}}-\left(\begin{array}{cc}{D}_u{\nu }_k+\mu & 0\\ 0& D_v{\nu }_k+\mu \end{array}\right)\right)=0,$

which gives the characteristic polynomial:

${\mu }^2-\left[\text{tr}\left(J^{\text{*}}\right)-\left(D_u+D_v\right){\nu }_k\right]\mu +\left[\text{det}\left(J^{\text{*}}\right)-\left(D_u{J}_{22}+D_v{J}_{11}\right){\nu }_k+D_u{D}_v{\nu }_k^2\right]=0.$ (47)

By the Routh-Hurwitz criterion, both roots of (47) have negative real part if and only if the conditions (46) hold. □

Corollary 5.6 (Turing Instability Condition). If $\text{tr}\left(J^{*}\right)<0$ and $\text{det}\left(J^{*}\right)>0$ (stable without diffusion), but there exists ${\nu }_k>0$ such that the second condition in (46) fails, then the system (32) undergoes a Turing-type diffusion-driven instability [16]. This would correspond, in the educational context, to the spontaneous emergence of spatial heterogeneity in knowledge distribution—i.e., the formation of learning clusters—driven by the differential diffusion of cognitive and non-cognitive competences ($D_u\ne D_v$ ).

Remark 5.7. The possibility of Turing instability in the extended model (32) is pedagogically significant: it suggests that, even in a homogeneous initial population, differential diffusion of knowledge and motivation can spontaneously generate spatial patterns of competence—corresponding to the emergence of high-performing and low-performing subgroups within the learner population. This phenomenon, well-documented in educational sociology [17], here receives a mathematical underpinning.

5.6. Numerical Simulation of the Extended System

We simulate system (32) on $\Omega =\left[0,1\right]$ with the same spatial and temporal discretization as Section 4 (, $\text{Δ}t=0.05$ , $T=50$ weeks). The parameter values, updated with respect to the preliminary discussion in Section 2.7, are given in Table 5.

Table 5. Parameter values for the extended two-compartment simulation.

Parameter	Description	Value
$D_u$	Diffusion-cognitive competences	0.03
$D_v$	Diffusion-non-cognitive competences	0.08
$r_u$	Growth rate-cognitive	0.55
$r_v$	Growth rate-non-cognitive	0.40
$a$	Self-coupling exponent ($u$ )	0.5
$b$	Cross-coupling exponent ($v\to u$ )	0.5
$c$	Cross-coupling exponent ($u\to v$ )	0.5
$d$	Self-coupling exponent ($v$ )	0.5
${\lambda }_u$	Forgetting rate-cognitive	0.05
${\lambda }_v$	Disengagement rate-non-cognitive	0.02

Remark 5.8. The coupling exponents are set to $a=b=c=d=0.5$ (square-root coupling) rather than integer values, so as to ensure that the cross-coupling terms $u^{0.5}{v}^{0.5}$ remain numerically active even for moderate initial values. This choice is consistent with the Cobb-Douglas-type interaction functions used in economic growth models [2] and avoids the numerical extinction phenomenon that arises when integer exponents are combined with near-zero initial conditions. The conditions for the existence of nontrivial equilibria (38) are verified: $r_u/{\lambda }_u=11\gg 1$ and $r_v/{\lambda }_v=20\gg 1$ .

The initial conditions are set to:

$\begin{array}{l}{u}_0\left(x\right)=0.30+0.25\exp \left(-\frac{{\left(x-0.3\right)}^2}{2\times {0.08}^2}\right),\\ v_0\left(x\right)=0.35+0.20\exp \left(-\frac{{\left(x-0.6\right)}^2}{2\times {0.10}^2}\right),\end{array}$ (48)

representing a population with moderate baseline cognitive mastery (30%) and non-cognitive engagement (35%), with two spatially separated localized peaks: an advanced cognitive subgroup near $x=0.3$ and a highly motivated subgroup near $x=0.6$ .

5.7. Results of the Extended Simulation

5.7.1. Final Spatial Distributions

The left panel of Figure 2 displays the final spatial distributions $u\left(x,T\right)$ (solid blue) and $v\left(x,T\right)$ (dashed red), alongside the respective initial conditions $u_0$ (blue dotted) and $v_0$ (red dotted).

Three observations are noteworthy. First, both competence fields have achieved near-complete spatial homogenization at $t=T=50$ weeks: $u\left(x,T\right)\approx 0.911$ and $v\left(x,T\right)\approx 0.949$ uniformly across $\Omega =\left[0,1\right]$ . The two localized initial peaks —the cognitive peak at $x=0.3$ and the motivational peak at $x=0.6$ —have been fully absorbed into the spatially homogeneous equilibrium, driven by the diffusion operators $D_u\text{Δ}u$ and $D_v\text{Δ}v$ .

Second, the non-cognitive field $v$ reaches a higher equilibrium value ($v^{\text{**}}\approx 0.949$ ) than the cognitive field $u$ ($u^{\text{**}}\approx 0.911$ ). This is consistent with the parameter values: the disengagement rate ${\lambda }_v=0.02$ is four times lower than the cognitive forgetting rate ${\lambda }_u=0.05$ , and the non-cognitive diffusion coefficient $D_v=0.08$ is nearly three times larger than $D_u=0.03$ , facilitating faster and more efficient spread of motivation through the learner population.

Third, the initial spatial separation between the cognitive peak ($x=0.3$ ) and the motivational peak ($x=0.6$ ) has no visible effect on the final distribution, demonstrating the robustness of the model to the spatial arrangement of initial conditions.

5.7.2. Temporal Dynamics of Mean Competences

The central panel of Figure 2 shows the temporal evolution of $\overline{u}\left(t\right)$ (solid blue) and $\overline{v}\left(t\right)$ (dashed red). Both curves exhibit a rapid sigmoidal convergence toward their respective equilibria, with the following key features:

(i) Rapid convergence. Both competences reach their equilibrium values within approximately 10 weeks—significantly faster than the single-compartment model under the best pedagogical scenario ($M_5$ : $t_{1/2}\approx 4$ weeks for $\overline{u}$ ). This acceleration is a direct consequence of the bidirectional coupling: the initial motivational level $\overline{v}\left(0\right)\approx 0.40$ boosts the cognitive growth term $r_u{u}^{0.5}{v}^{0.5}\left(1-u\right)$ from the very start, while the cognitive progress simultaneously reinforces motivation through the term $r_v{u}^{0.5}{v}^{0.5}\left(1-v\right)$ .

(ii) Non-cognitive leadership. The non-cognitive variable $\overline{v}\left(t\right)$ consistently precedes $\overline{u}\left(t\right)$ throughout the transient phase, reaching $\overline{v}=0.90$ at approximately $t=8$ weeks while $\overline{u}$ reaches the same level only around $t=10$ weeks. This two-week lead reflects the higher diffusion coefficient $D_v=0.08>D_u=0.03$ and the lower disengagement rate ${\lambda }_v=0.02<{\lambda }_u=0.05$ : motivational states spread faster and decay slower than specific knowledge items, consistent with the educational psychology literature [13] [15].

(iii) Equilibrium values. The horizontal dotted lines confirm the asymptotic convergence to $u^{\text{**}}\approx 0.911$ and $v^{\text{**}}\approx 0.949$ , in agreement with Proposition 5.3 (existence of nontrivial equilibria under (38)) and Theorem 5.5 (local asymptotic stability).

(iv) Mutualistic amplification. Both equilibrium values ($u^{\text{**}}=0.911$ , $v^{\text{**}}=0.949$ ) substantially exceed what each compartment would achieve independently. Indeed, setting $b=d=0$ (no coupling), the single-compartment cognitive model with $r_u=0.55$ and ${\lambda }_u=0.05$ would yield a lower equilibrium, since the cross-coupling term $v^{0.5}$ provides a persistent multiplicative boost throughout the dynamics. This mutualistic amplification effect is the central pedagogical message of the extended model: interventions that simultaneously target cognitive and non-cognitive competences produce synergistic gains that exceed the sum of their individual contributions.

5.7.3. Phase Portrait Analysis

The right panel of Figure 2 displays the phase portrait of the mean trajectory $\left(\overline{u}\left(t\right),\overline{v}\left(t\right)\right)$ in the $\left(u,v\right)$ -plane. The trajectory (solid black curve) originates from the initial state $\left(\overline{u}\left(0\right),\overline{v}\left(0\right)\right)\approx \left(0.34,0.40\right)$ (green circle) and converges monotonically to the nontrivial equilibrium $\left(u^{\text{**}},v^{\text{**}}\right)\approx \left(0.911,0.949\right)$ (red star).

Several structural features of the phase portrait deserve attention:

(i) Monotone trajectory. The trajectory is strictly monotone increasing in both coordinates: $\overline{u}$ and $\overline{v}$ increase simultaneously throughout the dynamics, with no oscillation or overshoot. This is consistent with the mutualistic (cooperative) nature of the coupling: both variables reinforce each other and neither can decrease while the other increases.

(ii) Concave curvature. The trajectory is concave, bending toward the $\overline{v}$ -axis in the early phase. This reflects the faster initial growth of the non-cognitive variable ($D_v>D_u$ , ${\lambda }_v{\lambda }_u$ ): motivation rises more rapidly than cognitive mastery in the transient phase, pulling the trajectory upward before it curves toward the equilibrium.

(iii) Stability confirmation. The convergence to the single point $\left(u^{\text{**}},v^{\text{**}}\right)$ with no limit cycle or divergence confirms the local asymptotic stability of this equilibrium, as established by Theorem 5.5 under the Routh-Hurwitz conditions (46).

(iv) Absence of Turing instability. The spatial homogeneity of the final distributions (left panel) confirms that the parameter configuration used here does not trigger the Turing diffusion-driven instability described in Corollary 5.6: the diffusion ratio $D_v/D_u=8/3\approx 2.67$ is insufficient to destabilize the spatially uniform equilibrium for the chosen Jacobian values.

Figure 2. Numerical simulation of the extended two-compartment model (32) over $T=50$ weeks. Left: final spatial distributions $u\left(x,T\right)$ (solid blue, cognitive) and $v\left(x,T\right)$ (dashed red, non-cognitive); dotted curves show the respective initial conditions. Center: temporal evolution of mean competence levels $\overline{u}\left(t\right)$ and $\overline{v}\left(t\right)$ ; dotted horizontal lines mark the equilibria $u^{\text{**}}\approx 0.911$ and $v^{\text{**}}\approx 0.949$ . Right: phase portrait of the mean trajectory $\left(\overline{u}\left(t\right),\overline{v}\left(t\right)\right)$ , from initial state (green circle) to nontrivial equilibrium (red star). Parameter values as in Table 5.

6. Conclusions

6.1. Summary of Contributions

This paper has proposed and analyzed a reaction-diffusion partial differential equation framework for modeling the spatiotemporal propagation of knowledge within a structured learner population. The central object of study is the knowledge density field $u\left(x,t\right)\in \left[0,1\right]$ , governed by the nonlinear PDE (13), whose structure integrates four fundamental mechanisms of educational dynamics in a unified mathematical setting: (i) Social contagion of knowledge via the diffusion operator with pedagogically modulated coefficient; (ii) Nonlinear learning dynamics via a logistic-Hill reaction term combining saturation, threshold effects, and social facilitation; (iii) Pedagogical modulation of the learning rate encoding the comparative efficiencies of five teaching methods and the inhibitory effect of cognitive stress; (iv) Knowledge decay via a linear forgetting term consistent with the Ebbinghaus forgetting curve.

The mathematical analysis established three levels of theoretical results. At the well-posedness level, Theorem 3.2 guarantees the existence and uniqueness of a weak solution, and Theorem 3.3 confirms that the interval $\left[0,1\right]$ is a positively invariant region—a fundamental requirement for pedagogical consistency. At the equilibrium level, Propositions 3.5 and 3.6 identify the trivial equilibrium $u^{\text{*}}=0$ (total knowledge loss) and the nontrivial equilibrium $u^{\text{**}}\in \left(\theta ,1\right)$ (sustained mastery), the latter existing if and only if the net learning rate exceeds the forgetting rate. At the stability level, Theorem 3.10 provides an explicit condition for the local asymptotic stability of $u^{\text{**}}$ , with the pedagogically meaningful corollary that high-mastery populations are self-stabilizing.

The fictitious numerical simulation validated all theoretical predictions and produced quantitative insights into the comparative impact of the five pedagogical methods. The ranking—personalized tutoring > active learning > collaborative learning > digital resources > lecture-based teaching—is unambiguous and consistent with the meta-analytic literature. The extended two-compartment model, coupling cognitive competences $u$ and non-cognitive competences $v$ (motivation, self-efficacy), revealed a mutualistic amplification phenomenon: the bidirectional coupling produces equilibrium values that substantially exceed what either compartment would achieve independently, illustrating the synergistic gains of interventions that simultaneously address cognitive and motivational dimensions of learning.

6.2. Limitations and Perspectives

Despite its theoretical richness, the present model rests on several simplifying assumptions that open natural directions for future research. The most immediate perspective is the calibration of model parameters from real educational datasets. A natural extension is to formulate an optimal control problem for the spatio-temporal allocation of teaching resources. Graph-theoretic generalizations replacing the continuous Laplacian by a graph Laplacian would allow the model to capture the topology of real learner interaction networks. Stochastic PDE generalizations would allow for uncertainty quantification. Finally, the possibility of Turing-type diffusion-driven instability in the extended system provides a mathematical explanation for the spontaneous formation of learning clusters observed in educational sociology.

6.3. Concluding Remarks

The reaction-diffusion framework developed in this paper provides, to the best of our knowledge, the first unified PDE model that simultaneously accounts for social knowledge contagion, pedagogical method differentiation, cognitive stress, Hill-type threshold effects, nonlocal social facilitation, Ebbinghaus forgetting, and the coupling between cognitive and non-cognitive competences. Its mathematical tractability—well-posedness, equilibrium analysis, stability theory—combined with its numerical implementability and natural calibration pathway, makes it a versatile and operationally useful tool for educational science.