Stochastic response of the Volterra integral hysteresis model via a gaussian basis function solution

1. Introduction

Hysteretic nonlinear systems have broad applications in various engineering domains, including viscoelastic composite dampers [1] and other related fields [2-3]. For example, Humaidi et al. developed various advanced nonlinear control schemes for magnetorheological dampers (MRDs) in the seismic vibration control of buildings [4-6]. The common feature of hysteretic nonlinear systems is that their mechanical response depends not only on the current input but also on the past loading history, which results in a memory effect.

Several models have been developed for hysteretic nonlinear systems, including the Bouc-Wen model [7], the bilinear model [8], the Duhem model [9], the Preisach model [10], among others. The Volterra integral model represents one such approach and stands out because of its clear physical basis and flexible parameter tunability. It employs a convolution integral with an exponentially decaying kernel, expressing the current stress as the sum of an instantaneous nonlinear strain term and a weighted cumulative term that accounts for the history of strain rates, thereby naturally capturing the material's memory. Its integral-type constitutive relation is physically transparent, and by adjusting the memory kernel, various hysteresis loop shapes can be reproduced. Furthermore, unlike several conventional models, the model is characterized by a relatively concise set of parameters, and its formulation does not require piecewise functions.

When such systems are subjected to random environmental excitations, their responses become stochastic and must be treated within a probabilistic framework. The evolution of the response probability density function is governed by the Fokker-Planck-Kolmogorov (FPK) equation. For the Volterra integral model with an exponentially decaying kernel, the FPK equation is a standard partial differential equation, and the Volterra integral model has the advantage of being transformable into a standard three-dimensional Markov system. However, even this FPK equation remains difficult to solve analytically. Traditional approximate methods (such as equivalent linearization [11]-[13], the moment method [14], and the stochastic averaging method [15]-[17]) often introduce significant errors when dealing with such strongly memory-dependent systems because they either ignore the cumulative historical effects or oversimplify the kernel function structure. This limitation is particularly severe when the kernel decay rate is moderate or small, severely restricting the applicability of these methods.

In recent years, semi-analytical methods based on the least-squares fitting approach with radial basis functions have shown great potential in stochastic dynamics [18] and other related fields [19]-[25]. These methods represent the unknown probability density function as a weighted sum of Gaussian basis functions and determine the weight coefficients by minimizing the residual sum of squares of the FPK equation at sample points, thereby transforming the solution of the partial differential equation into an algebraic optimization problem. This approach exhibits good robustness to the type and strength of system nonlinearity. It has been successfully applied to the stochastic response analysis of various hysteretic and strongly nonlinear systems. Chen et al. used this method to analyze the response of Duhem models under harmonic and random excitations [26], as well as the stochastic response of Bouc-Wen systems when excited by Poisson white noise [27].

In this paper, we investigate the steady-state stochastic response of a hysteretic system described by the Volterra integral model under Gaussian white noise excitation. First, we establish the differential equation of motion of the system and analyze the influence of the memory kernel parameters on the area and morphology of the hysteresis loop. Then, based on the Gaussian basis function (GBF) fitting method, we derive semi-analytical expressions for the stationary joint probability density function and the marginal PDFs of the system. Finally, a comparison with MCS is performed. This study aims to 1) reveal the non-monotonic regulatory mechanisms of key parameters – such as the memory kernel decay rate, the memory damping strength, and the nonlinear stiffness – on the response probability distribution and to uncover the physical mechanism of “memory-induced stiffness” induced by the cumulative effect of historical strain rates when the kernel decay rate is extremely small; 2) verify the computational accuracy and applicability of the GBF method for the Volterra integral hysteresis model; and 3) provide methodological support for the stochastic vibration analysis and parameter optimization of hysteretic models with strong memory effects.

2. Model

2.1. Volterra integral model

This paper introduces a theoretical stress–strain model with hysteretic characteristics. The governing equation is shown in Eq. (1):

1

$\sigma =E\epsilon +A{\epsilon }^3+{\int }_0^t{e}^{-\lambda (t-\tau )}\dot{\epsilon }\left(\tau \right)d\tau ,$

where $\sigma$ represents stress, $\epsilon$ represents strain, $E$ is the linear elastic modulus, $A$ represents the nonlinear elastic modulus, ${\int }_0^t{e}^{-\lambda (t-\tau )}\dot{\epsilon }\left(\tau \right)d\tau$ represents the Volterra integral. $e^{-\lambda (t-\tau )}$ is the Volterra integral memory kernel function, $\lambda$ represents the decay rate of the Volterra kernel function, $t$ is the current time, and $\tau$ represents a past moment. The integral indicates the weighted influence of the strain rate at past moments on the current stress. The images of the stress-strain model are shown in Fig. 1.

In Fig. 1(b), the stress-strain relationship diagram exhibits a hysteresis loop, which can be interpreted as the current stress always being related to the strain over a past period.

Fig. 1Volterra integral of the stress-strain model

a) Strain-time and Stress-time curve

b) Stress-strain relationship curve (schematic)

2.2. Parameter influence on the Volterra integral

With the Volterra kernel function decay rate $\lambda =$ 0.15 and the linear elastic modulus $E$ taken as 100 and 75, the Volterra integral hysteresis loops are plotted in Fig. 2. When the linear elastic modulus $E$ is taken as 100 and the Volterra kernel function decay rate $\lambda$ is taken as 0.15 and 0.075, the Volterra integral hysteresis loops are presented in Fig. 3(a); when $\lambda$ is taken as 1.25 and 2, the hysteresis loops are presented in Fig. 3(b); when $A$ is taken as 0 and 7e3, the hysteresis loops are presented in Fig. 4.

Fig. 2Volterra integral hysteresis loops under different linear elastic moduli

Fig. 3Volterra integral hysteresis loops under different Volterra kernel function decay rates

a)

b)

Fig. 4Volterra integral hysteresis loops under different nonlinear elastic moduli

Figs. 2-4 show that when the linear elastic modulus $E$ decreases, the hysteresis loop undergoes a marked change in shape, the material stiffness decreases, and the overall vertical position of the loop decreases (because the elastic part diminishes), but the hysteresis loop exhibits no change in width. That is, when $E$ decreases, the hysteresis loop exhibits a constant area, and the loop shifts downward overall in the stress direction, resulting in a flatter profile, but the dissipated energy remains unchanged. When $\lambda$ is very small, as the Volterra kernel function decay rate $\lambda$ increases, the hysteresis loop substantially increases in area. However, as $\lambda$ continues to increase, the hysteresis loop area decreases instead. As the nonlinear elastic modulus $A$ increases, the area of the hysteresis loop remains unchanged, but its shape transforms into an “S” shape. It is evident that the linear elastic modulus $E$ and the nonlinear elastic modulus $A$ are key factors determining the shape of the Volterra integral hysteresis loop. The area of the Volterra integral hysteresis loop is determined by the Volterra kernel function decay rate $\lambda$, and it also depends on loading characteristics.

3. Methodology

We employ the RBF method [28], where the output is formed by summing $N$ Gaussian basis functions each multiplied by its corresponding weight coefficient. This sum is used to represent the steady-state solution of Eq. (9), denoted as $\tilde{P}$:

where the subscript $j$ denotes the $j$-th Gaussian basis function. $w$ represents the weight coefficient, $\sigma$ denotes the standard deviation, $\varphi \left(\mathbf{q},{\mathit{\mu }}_j,{\sigma }_j\right)$ denotes a Gaussian basis function:

where ${q=\left[x,y,z\right]}^T$, ${\mathbf{\mu }}_j={\left[{\mu x}_j\text{,}{\mu y}_j\text{,}{\mu z}_j\right]}^T$, $‖\cdot ‖$ is the Euclidean norm, representing the distance between the input coordinate $\mathbf{q}$ and the center point coordinate ${\mathbf{\mu }}_j$. We select a domain ${\mathrm{\Omega }}_{\varphi }=\left[-d_c,d_c\right]\times \left[-d_c,d_c\right]\times \left[-d_c,d_c\right]={\left[-\mathrm{1,1}\right]}^3$ in the $x-y-z\in R^3$ space, then divide it into 18×18×18 grid points, with each grid point serving as a GBF center, i.e., $N$ = ${18}^3$. ${\sigma }_j$ represents the variance, and the smaller it is, the sharper the function graph becomes. Setting the basis function variance to ${\sigma }_j=\frac{d_{max}}{\sqrt{2N}}$ yields a moderate shape, where $d_{max}$ denotes the maximum distance between any two center points. This choice ensures sufficient overlap between adjacent basis functions, thereby avoiding numerical ill-conditioning.

Fig. 5Flowchart of the proposed solution procedure

Define a function as the residual between $\frac{\partial P}{\partial t}$ and $\frac{\partial \tilde{P}}{\partial t}$:

Plugging Eq. (10) into Eq. (12), the residual $\xi \left(\mathbf{q},\mathbf{w}\right)$ can be rewritten as:

where $s_j$ denotes the error function.

The following normalization condition holds for all $j$ between $\tilde{P}$ and the GBF:

The conditions on the weight coefficients are:

We define a sampling domain ${\mathrm{\Omega }}_S=\left[-d_s,d_s\right]\times \left[-d_s,d_s\right]\times \left[-d_s,d_s\right]$. Here $N_s$ denotes the total number of sampling points. Typically, $d_s$ is between $d_c$ and 2.5$d_c$ to ensure that the sampling points completely cover the effective support regions of all basis functions. This is a common practice to guarantee a well-conditioned interpolation matrix, so ${\mathrm{\Omega }}_S={\left[-\mathrm{2,2}\right]}^3$. We set $N_s$ = ${36}^3$. The sampling density is far higher than the theoretical threshold required to ensure convergence of the discrete approximation to the continuous solution, which guarantees the accuracy of the approximation. The Lagrangian function $L_s$ is then constructed based on the FPK residual evaluated at randomly sampled collocation points:

The necessary condition for minimizing ${\mathrm{L}}_s\left(\mathbf{w}{,k}^{\text{'}}\right)$ is:

They are expanded into matrix form as:

where the error matrix $\mathbf{S}{\mathbf{S}}^T$ is expressed as:

The optimal weight coefficient is represented in vector form:

Substituting the $N_s$ into the error function $s_j\left({\mathbf{q}}_i\right)$ in Eq. (21) and calculating the vector in Eq. (22), we obtain the weight coefficients.

We acquire the approximate PDF $\tilde{P}$. Then, integrating $\tilde{P}$ yields the marginal PDF:

Here, we adopt the root mean square error (RMS) of the PDF solution as the metric to quantify the error of the Gaussian basis function fitting:

where $p_i$ is the PDF value from MCS (a 4th-order stochastic Runge-Kutta algorithm [29]) at point $i$; $\stackrel{-}{p_i}$ is the corresponding value from the RBF method. For Eq. (2), 3000 independent noise realizations are generated. Each noise sequence contains 10,000 points, with the last 6,000 points of the response considered the steady-state response. Then, the simulation time step is $\mathrm{\Delta }t=$0.04. The final 3,000×6,000 $\left(x,y\right)$ simulated points serve to statistically determine the probability of each grid falling within the $x$-$y$ plane. The following sections present comparisons between the GBF and MCS results. The computed results are listed in Table 1.

4. Numerical simulation and error analysis

4.1. Effect of the Volterra kernel decay rate

With parameters fixed at $\alpha$ = 1, $\beta$ = 0.15, $\lambda$ is taken as 1.25 and 1.35 to examine how $\lambda$ influences the system when it takes large values.

From Fig. 7, when the decay rate $\lambda$ of the Volterra memory kernel is relatively large, a further increase in $\lambda$ leads to consistent variations between the PDF profiles and the hysteresis loops. The underlying physical mechanism is that excessively rapid decay of the kernel restricts significant contributions to the current stress solely to strain rates at very recent past instants. Consequently, the equivalent damping supplied by the memory damping term diminishes, reducing the system’s overall energy dissipation capacity, which in turn amplifies the amplitudes of both the displacement and velocity responses. On the steady-state probability density function, this behavior manifests as a lowered peak and broadened tails. Thus, it can be concluded that in this case, the Volterra kernel function decay rate $\lambda$ is negatively correlated with the hysteresis loop area.

Fig. 6Joint probability density plots for different values of the Volterra kernel function decay rate

a)$\lambda$ = 1.25 RBFNN outcomes

b)$\lambda$ = 1.25 Monte Carlo simulation

c)$\lambda$ = 1.35 RBFNN outcomes

d)$\lambda$ = 1.35 Monte Carlo simulation

Fig. 7Influence of the Volterra kernel function decay rate on the displacement and velocity probability densities: a) shows px, b) shows py. Solid lines represent RBF, circles represent MCS

a)

b)

Parameters are fixed at $\alpha$ = 1, $\beta$ = 0.15, with $\lambda$ taken as 0.15 and 0.03 to examine how the decay rate of the Volterra kernel function influences the system when $\lambda$ is small.

From Fig. 9, it can be seen that when $\lambda$ is very small, the system has a long memory time, and the stress is influenced by the deformation history from much earlier times, potentially retaining “memory” of strain rates from multiple past cycles. This leads to stress responses that are not only related to current strain but also to multiple historical cycles. The hysteresis loop becomes wider or exhibits complex time-delayed responses, and the energy dissipation characteristics vary more gradually with frequency. That is, a smaller $\lambda$ leads to lower peaks and a wider distribution range in the displacement and velocity probability density images. For an extremely small $\lambda$ of the Volterra memory kernel, the memory term ceases to function solely as a damping element and approximates a history-dependent elastic constraint. Thus, it can be concluded that in this case, the Volterra kernel function decay rate $\lambda$ is positively correlated with the hysteresis loop area.

The decay rate of the Volterra memory kernel governs the fading rate of the material’s memory of past strain rates, and it exerts an inherently non-monotonic influence on both the hysteresis loop area and the system’s steady-state response. For a single exponential kernel, this non-monotonic behavior is related to its frequency-dependent equivalent stiffness and damping. This non-monotonic effect, by modifying the probability density function shape, can provide a basis for parameter design in structural vibration control and also serve as an indicator for assessing material aging conditions.

Fig. 8Joint probability density plots for different values of the Volterra kernel function decay rate

a)$\lambda$ = 0.15 RBFNN outcomes

b)$\lambda$ = 0.15 Monte Carlo simulation

c)$\lambda$ = 0.03 RBFNN outcomes

d)$\lambda$ = 0.03 Monte Carlo simulation

Fig. 9Influence of the Volterra kernel function decay rate on the displacement and velocity probability densities: a) shows px, b) shows py. Solid lines represent RBF, circles represent MCS

a)

b)

4.2. Effect of memory damping strength coefficient

Parameters are fixed at $\lambda$ = 0.15, $\alpha$ = 1, with $\beta$ taken as 0.15 and 0.25 to examine how the memory damping strength coefficient influences the system.

From Figs. 10-11, with increasing memory damping strength coefficient $\beta$, the peaks of the displacement and velocity PDFs increase, whereas the bell-shaped distribution ranges become narrower. Thus, it can be concluded that a larger memory damping strength coefficient $\beta$ results in higher peaks.

Fig. 10Steady-state probability density plots for two different memory damping strength coefficient values

a)$\beta$= 0.15 RBF outcomes

b)$\beta$= 0.15 Monte Carlo simulation

c)$\beta$ = 0.25 RBF outcomes

d)$\beta$= 0.25 Monte Carlo simulation

Fig. 11Probability density of displacement and velocity for two different memory damping strength coefficient values: a) shows px, b) shows py. Solid lines represent RBF, circles represent MCS

a)

b)

4.3. Effect of the nonlinear stiffness coefficient

Parameters are fixed at $\lambda$ = 0.15 and $\beta$ = 0.15, with $\alpha$ taken as 1 and 1.8 to examine the effect of the nonlinear stiffness coefficient on the steady-state response.

The nonlinear stiffness coefficient $\alpha$ is the dimensionless form of the cubic nonlinear coefficient $A$ defined in Section 2.1. It changes the stiffness distribution of the system, without altering the dissipative characteristics of the Volterra hysteresis term, but it significantly affects the probabilistic distribution characteristics of the response. Fig. 13 shows that as the nonlinear stiffness coefficient $\alpha$ increases, the peaks of the curves increase significantly, and the coverage range of the curves becomes narrower. Thus, it can be concluded that a larger nonlinear stiffness coefficient $\alpha$ results in higher peaks and has a relatively large influence on the steady-state PDF of displacement and velocity.

Fig. 12Joint probability density plots for different values of the nonlinear stiffness coefficient

a)$\alpha$= 1 RBF outcome

b)$\alpha$= 1 Monte Carlo simulation

c)$\alpha$= 1.8 RBF outcomes

d)$\alpha$= 1.8 Monte Carlo simulation

Fig. 13Probability densities p(x) and py for different values of the nonlinear stiffness: a) shows px, b) shows py. Solid lines represent RBF, circles represent MCS

a)

b)

5. Conclusions

Addressing the response problem of hysteretic nonlinear systems under random excitation, this paper describes their constitutive relationship based on the Volterra integral model. The GBF fitting method is extended to this case. The summary is as follows:

1) In Table 1, the proposed GBF approach demonstrates reliability in analyzing Volterra integral hysteretic nonlinear oscillators.

2) The physical meaning of the hysteresis loop parameters in the Volterra integral model is clear. The decay rate $\lambda$ plays a key role in determining the hysteresis loop area. Changes in the linear elastic modulus $E$ do not alter the dissipated energy, leaving the hysteresis loop area constant, but the loop shifts overall in the stress direction; the nonlinear elastic modulus $A$ does not change the area of the hysteresis loop, but it does change its shape.

3) The decay rate $\lambda$ of the Volterra kernel exerts a non-monotonic regulatory effect on the system response. When $\lambda$ is relatively large, the material exhibits short memory retention, causing the hysteresis loop to contract, so the peak value of the PDF decreases. Conversely, when $\lambda$ is small, long-term memory effects become significant: the accumulation of historical strain rates effectively enhances the system’s equivalent stiffness, which suppresses the response amplitude and leads to an increase in the PDF peak. This phenomenon is referred to as the “memory-induced stiffness” effect.

4) $A$ represents the cubic nonlinear term in the constitutive equation, and $\alpha$ is its dimensionless nonlinear stiffness coefficient. They do not alter the dissipative characteristics of the system. However, it does change the stiffness distribution of the system, thereby significantly affecting the probabilistic characteristics of the response. As the nonlinear stiffness coefficient $\alpha$ increases, the peak becomes higher, and the distribution range becomes narrower. As $\alpha$ increases from 1 to 1.8, $p\left(x\right)$ increases from 0.87 to 0.96, and $p\left(y\right)$ increases from 0.74 to 0.78. When $\beta$ increases from 0.15 to 0.25, $p\left(x\right)$ increases from 0.87 to 0.94, and $p\left(y\right)$ increases from 0.74 to 0.8.

Future work will extend the proposed method to the dynamic analysis of nonlinear mechatronic systems with hysteresis and uncertainties [30].