Numerical Solutions for Fuzzy Stochastic Ordinary Differential Equations Using Heun’s Method with a Dual-Wiener Process Framework

In 1905, Einstein first introduced the concept of stochastic differential equations (SDEs), providing a mathematical description of Brownian motion by demonstrating the connection between the random behavior of small particles and the diffusion equation on a larger scale. Since then, these equations have become an indispensable tool in various fields such as physics, chemistry, biology and microelectronics, as well as in economics and finance. Traditionally, solutions of these equations have relied heavily on the Ito integration technique to achieve accuracy when dealing with stochastic phenomena. However, the exact methods face some challenges in solving non-trivial problems, which led to the resort to approximation methods as an alternative to overcome these difficulties, as discussed in study [1].

In many practical applications, data on parameters may be incomplete or distorted by measurement errors or mathematical approximations. These uncertainties lead to ambiguities in the differential equations used, which require more sophisticated methods for estimating the solutions. These uncertainties can be expressed using probabilistic or fuzzy for the simulation of complex natural phenomena [2-4].

Heun’s method, or Heun’s consistency method, is a numerical procedure for solving ordinary differential equations (ODEs) and is ideal for the order initial value problem. Compared to the simple Euler method, it is superior because of basic prediction and correct step strategies, where an interval with basic gradient is used and another interval is used to adjust constant gradients of that respective interval. Here, the constructive solution and a method that fits well are Heun’s method because of the lack of constructive solutions to fuzzy stochastic ordinary differential equations (FSODEs). In order to manage the existing uncertainty due to fuzzy parameters, the Heun’s method is integrated into the dual Wiener process within the course of the study. This makes it possible as an efficient numerical solution can be achieved, which could be utilized in several practical area of random occurrences and vagueness. Instead of that, this approach enhances the effectiveness of the solutions as well as the applicability of the numerical, as well as ansatz, methods in complex modeling situations [5-14].

In the past years, the combining of fuzziness with stochastic problems modeled by differential equations has earned attention in several fields of science, engineering, and finance. The fuzzy differential equations (FDEs) process the uncertainty that is found in real-world problems and can't be sufficiently represented by classical probabilistic models. Such models merge the flexibility of fuzzy theory and the randomness in the stochastic models.

The solutions key of FSDEs pass essentially through numerical methods that are used to derive approximate solutions. Among these methods, we focus our attention on the Heun's method, one of the second order Runge-Kutta schemes, which shows effectiveness in tackling generalized FSDEs. With low cost and high accuracy, this method is considered as an improvement of the Euler method. As well, the Heun's method assists in approximating the evolution of fuzzy quantities over time.

A critical component of stochastic models is the dual Wiener process, which expands the classical Brownian motion. In the study of complex systems, the dual Wiener process is an appropriate in modeling systems when multiple uncertainties need to be hold. In this analysis, we aim to examine the applicability of the Heun's method to FDEs obtained by dual Wiener processes.

3.1 Fuzzy sets and numbers

The fuzzy set $\check{A}$ contains pairs of elements and their associated membership functions. The membership function $\mu_{\tilde{A}}$ is defined as follows:

$\mu_{\tilde{A}}: X \rightarrow[0,1]$              (1)

where, X is the universal set. A fuzzy number is defined as a convex natural set and is given as follows: The fuzzy number z in parametric form is a pair of values $[\underline{k}, \bar{k}]$ of functions $\underline{k}(\alpha), \bar{k}(\alpha)$, where $\alpha \in[0,1]$, which meets the following conditions [1]:

1. $\underline{k}(\alpha)$ it is a left-continuous, non-decreasing function bounded in the interval [0, 1], and right continuous at 0.

2. $\bar{k}(\alpha)$ it is a left-continuous, non-increasing function bounded in the interval (0,1], and right continuous at 0.

3. $\underline{k}(\alpha) \leq \bar{k}(\alpha), \quad 0 \leq \alpha \leq 1$.

A fuzzy number $\check{A}=[a, b, c]$ is said to be triangular when its membership function is given by:

$\mu_{\check{A}}(x)=\left\{\begin{array}{cl}0 & \text { if } x \leq a \\ \frac{x-a}{b-a} & \text { if } a<x \leq b \\ \frac{c-x}{c-b} & \text { if } b<x<c \\ 1 & \text { if } \quad x \geq c\end{array}\right.$                 (2)

Consider the two fuzzy numbers $\tilde{L}=[\underline{L}(\alpha), \bar{L}(\alpha)]$ and $\widetilde{D}=[\underline{D}(\alpha), \bar{D}(\alpha)]$ and a scalar d then:

1. $\widetilde{L}=\widetilde{D}$ if and only if $\underline{L}(\alpha)=\underline{D}(\alpha)$ and $\bar{L}(\alpha)=\bar{D}(\alpha)$.

2. $\widetilde{L}+\widetilde{D}=[\underline{L}(\alpha)+\underline{D}(\alpha), \bar{L}(\alpha)+\bar{D}(\alpha)]$.

3. $\widetilde{L}-\widetilde{D}=[\underline{L}(\alpha)-\underline{D}(\alpha), \bar{L}(\alpha)-\bar{D}(\alpha)]$.

4. $d \tilde{L}=\left\{\begin{array}{l}{[d \underline{L}(\alpha), d \bar{L}(\alpha)], \text { if } \mathrm{d} \geq 0} \\ {[d \underline{D}(\alpha), d \bar{D}(\alpha)], \text { if } d<0}\end{array}\right.$

3.2 SDEs

Let's review a classical SDE [1]:

$d x_t=a d t+b d w_1+k d w_2$                 (3)

where, Eq. (3) is expressed in differential form:

$x\left(t_2\right)-x\left(t_1\right)=\int_{t_1}^{t_2} a(t, x) d s+\int_{t_1}^{t_2} b(t, x) d w_1(s)+\int_{t_1}^{t_2} k(t, x) d w_2(s)$                (4)

The last term on the right-hand side of Eq. (2) is known as the Ito integral. We will now take $0=t_0<t_1<t_2<\cdots<t_n=T$, be a grid of points on an interval [0, T], the Ito integrals are determined for each component separately, depending on the different components on the right-hand side of Eq. (4) is being:

$\int_{t_1}^{t_2} b(t, x) d w(s)=\int_{t_0}^t b^{i, j}(s, x s) \Delta w_i$                   (5)

where, $\Delta w_i=w_{t_i}-w_{t_{i-1}}$, transition in Brownian motion over a time interval.

Let us consider Eq. (5), which is first solved analytically using Ito's formula. According to this formula, if $X_t$ represents the Ito process, then:

$d \xi(t)=a d t+b d w(t)$                 (6)

and let f(x, t) the function is continuous in $(x, t) \in \times[\infty, 0)$, along with its partial derivatives $f x, f x x, f_t$. Then the process $f(\xi(t), t)$. It is characterized by a random differential defined as follows:

$d f(\xi(t), t)=\left[f_t(\xi(t), t)+f x(\xi(t), t) a(t, x)+\frac{1}{2} f x x(\xi(t), t) b_2(t, x)\right] d t+f x(\xi(t), t) b(t, x) d w(t)$                  (7)

3.3 Computational techniques for SDEs

To calculate the numerical solution of a SDEs, we define a grid of points, $0=t_0<t_1<t_2<\cdots<t_n=T$, and approximate x values $w_0<w_1<w_2<\cdots<w_n$, it is determined at specific values of t. Consider the initial value problem of a SDE [1]:

$\left\{\begin{array}{c}d X_t=a\left(t, x_t\right) d t+b\left(t, x_t\right) d w_1+k\left(t, x_t\right) d w_2 \\ X(0)=X_0\end{array}\right.$               (8)

where, $a\left(t, x_t\right), b\left(t, x_t\right)$ and $k\left(t, x_t\right)$ are all continuous functions which are defined on interval $\left[t_0, T\right]$. Where $d w_1, d w_2$ are wiener process with components $w_t^1, w_t^2, \ldots, w_t^m$. Then Eq. (8) is solved numerically as follows.

3.4 Adaptation of Heun’s method for fuzzy parameters and dual-Wiener processes

The uncertainties associated with real-word applications of Heun’s method are handled by modifications to the method for fuzzy parameters and dual-Wiener processes. Heun’s method is a powerful predictor-corrector technique for solving numerical ODEs; the approach developed in this paper makes use of fuzzy numbers to describe uncertain parameters. This shift moves the original deterministic model to stochastic fuzzy environment and initial conditions and parameters are presented as triangular or trapezoidal fuzzy numbers [16-20].

Another improvement to the method is the addition of a combined dual-Wiener process approach that models stochastic system action using two separates but interacting Wiener processes. The predictor step of Heun’s method is modified to find solutions according to the fuzzy parameters and the corrector step amortize the solution by considering the stochasticity out of the dual-Wiener processes.

This combined adaptation enhances numerical precision and guarantees that the solutions captured consider uncertainty as well as stochastic nature making it advantageous in specific applications such as financial models, engineering failure analysis or environmental impact assessments.

Let us consider using a discrete-time approximation to estimate the Ito process in solving SDEs. Let us apply a discrete-time approximation to Eq. (8) to estimate the Ito process in solving SDEs.

$x\left(x_1\right)=x\left(x_0\right)+\int_{x_0}^{x_1} a\left(s, x_s\right) d s+\int_{t_0}^t b\left(s, x_{s_1}\right) d w_{s_1}+\int_{t_0}^t k\left(s, x_{s_2}\right) d w_{s_2}$                  (9)

Also let

$I_1=\int_{t_0=x_0}^{t=x_1} a\left(s, x_s\right) d s=\int_{t_0=x_0}^{t=x_1} f(x) d x$                  (10)

where, trapezoidal rule, $x_0=t_0, x_1=t, h=x_0-x_1=t_0-t$, and

$f(x)=P_n(x)+\frac{f^{n+1}(\delta)}{(n+1)!} \prod_{i=0}^n\left(x-x_i\right)$                (11)

Provided that

$P_n(x)=\sum_{\substack{i=0 \\ i \neq k}}^n f\left(x_i\right) * L_i(x)$               (12)

where, the Lagrange equation is $L_i(x)=\frac{X-X_i}{X_k-X_i}$. And, also if n=1. Then:

$P_1(x)=\sum_{\substack{i=0 \\ i \neq k}}^1 f\left(x_i\right) * \frac{X-X_i}{X_k-X_i}=\frac{X-X_1}{X_0-X_1} * f\left(x_0\right)+\frac{X-X_0}{X_1-X_0} * f\left(x_1\right)$               (13)

Hence:

$f(x)=\frac{x-x_1}{x_0-x_1} * f\left(x_0\right)+\frac{x-x_0}{x_1-x_0} * f\left(x_1\right)+\frac{f^{(n+1)}(\varepsilon)}{(n+1)!} \prod_{i=0}^n\left(x-x_i\right)$                    (14)

By substituting Eq. (14) into Eq. (10), we obtain:

$I_1=\int_{t_0=x_0}^{t=x_1}\left(\frac{x-X_1}{X_0-X_1} * f\left(x_0\right)+\frac{x-X_0}{X_1-X_0} * f\left(x_1\right)\right) d x+\int_{t_0=x_0}^{t=x_1} \frac{f^{(n+1)}(\varepsilon)}{(n+1)!} \prod_{i=0}^n\left(x-x_i\right) d x$                 (15)

Simplifying Eq. (14), we will get:

$I_1=\int_{t_0=x_0}^{t=x_1} f(x) d x=\int_{t_0=x_0}^{t=x_1}\left\{\frac{\left(x-x_1\right)^2}{2\left(x_0-x_1\right)} * f\left(x_0\right)+\frac{\left(x-x_0\right)^2}{2\left(x_1-x_0\right)} * f\left(x_1\right)\right\} d x=\frac{x_1-x_0}{2}\left[f\left(x_0\right)+f\left(x_1\right)\right]=\frac{h}{2}\left[f\left(x_0\right)+f\left(x_1\right)\right]$                   (16)

and

$I_2=\int_{t_0}^t b\left(x_{s_1}\right) d w_{s_1}=\int_{t_0=x_0}^{t=x_1} b\left(s, x_{s_1}\right) d w_{s_1}=\int_{t_0=x_0}^{t=x_1} b(x) d w_{s_1}=\frac{x_1-x_0}{2}\left[b\left(x_0\right)+b\left(x_1\right)\right]$                 (17)

Additionally,

$I_3=\int_{t_0}^t k\left(x_{s_2}\right) d w_{s_2}=\int_{t_0=x_0}^{t=x_1} k\left(s, x_{s_1}\right) d w_{s_2}=\int_{t_0=x_0}^{t=x_1} k(x) d w_{s_2}=\frac{X_1-X_0}{2}\left[k\left(x_0\right)+k\left(x_1\right)\right]$               (18)

The Heun’s approx. is defined as cont. Time stochastic process. $y=\left\{y(T) ; t_0 \leq t<T\right\}$ satisfying the iterative scheme: $X_0=x_0$:

$x_{n+1}=x_n+\frac{h}{2}\left[f\left(x_n\right)+f\left(x_n+a \Delta n+b \Delta w_n\right)\right] \Delta n+\frac{h}{2}\left[b\left(x_n\right)+b\left(x_n+a \Delta n+b \Delta w_n\right) b\left(x_n\right)+b\left(x_n+a \Delta n+b \Delta w_n\right)\right] \Delta w_{n_1}$

$+\frac{h}{2}\left[k\left(x_n\right)+k\left(x_n+a \Delta n+b \Delta w_n\right)\right] \Delta w_{n_2}$                  (19)

where,

$\begin{gathered}\Delta t_{i+1}=t_{i+1}-t_i, \\ \Delta W_{i+1}=\mathrm{W}\left(t_{i+1}\right)-W\left(t_i\right) \sim \sqrt{t_{i+1}-t_i} N(0,1)\end{gathered}$

where, N(0, 1) represents a standard normally distributed random variable with mean zero and variance one. The function randn (1, N), N random variables will be generated according to the standard normal distribution. To obtain a random variable with a given variance $\Delta t_{i+1}$, random variables with standard normal distribution are generated by a custom MATLAB function. randn (1, N). These variables are then multiplied by the result to obtain random increments in $\Delta W_{i+1}$.

3.5 Computational solutions for FSDEs

Suppose there are fuzzy parameters in a SDE; then Eq. (3) can be rewritten as:

$d[\underline{X}(\alpha), \bar{X}(\alpha)]=[\underline{a}(\alpha), \bar{a}(\alpha)] d t+[\underline{b}(\alpha), \bar{b}(\alpha)] d w_1+[\underline{k}(\alpha), \bar{k}(\alpha)] d w_2$              (20)

Eq. (10) is now solved using analytical and numerical methods respectively. By applying limit method, FSDE Eq. (10) can be reformulated explicitly and modified as follows [20]:

$d\left[\lim _{s \rightarrow \infty} X(\alpha), \lim _{s \rightarrow 1} X(\alpha)=\left[\lim _{s \rightarrow \infty} a(\alpha), \lim _{s \rightarrow 1} a(\alpha)\right] d t+\left[\lim _{s \rightarrow \infty} b(\alpha), \lim _{s \rightarrow 1} b(\alpha)\right] d w_1+\left[\lim _{s \rightarrow \infty} k(\alpha), \lim _{s \rightarrow 1} k(\alpha)\right] d w_2\right.$                (21)

where,

$\begin{gathered}X(\alpha)=\underline{X}(\alpha)+\frac{\bar{X}(\alpha)-\underline{X}(\alpha)}{s}, \\ a(\alpha)=\underline{a}(\alpha)+\frac{\bar{a}(\alpha)-\underline{a}(\alpha)}{s} \text { and } b(\alpha)=\underline{b}(\alpha)+\frac{\bar{b}(\alpha)-\underline{b}(\alpha)}{s}\end{gathered}$

First, for the precise case, we use explicit representations of $X(\alpha), a(\alpha)$  and $b(\alpha)$, use Ito's integral to solve the problem. This version is more concise and clearly conveys the use of crisp representations and the Ito integral for solving the problem. If we apply the fuzzy concept we discussed earlier to the Heun’s method, Eq. (9) can be reformulated as follows.

$\tilde{x}\left(x_1\right)=\tilde{x}\left(x_0\right)+\int_{x_0}^{x_1} a\left(s, \tilde{x}_s\right) d s+\int_{t_0}^t b\left(s, \tilde{x}_{s_1}\right) d w_{s_1}+\int_{t_0}^t k\left(s, \tilde{x}_{s_2}\right) d w_{s_2}$               (22)

The Heun’s approximation is defined as a continuous-time stochastic process, $y=\left\{y(T) ; t_0 \leq t<T\right\}$, that satisfies an iterative scheme. When applied to SDEs, the Heun's method is adapted to handle fuzzy parameters by incorporating triangular and trapezoidal fuzzy numbers. This approach enables the approximation of solutions to fuzzy SDEs, ensuring that the process accommodates the inherent uncertainty in the system. The resulting fuzzy Heun's method is essential for obtaining accurate solutions in systems where both randomness and fuzziness are presented:

$X_0(\alpha)=w_0(\alpha)$

$w_{i+1}(\alpha)=w_i(\alpha)+\frac{h}{2}\left[f\left(w_i\right)+f\left(w_i+a\left(t_i, w_i, \alpha\right) \Delta t_{i+1}\right)+b\left(t_i, w_i, \alpha\right) \Delta w_i\right] \Delta t_{i+1}$

$\begin{aligned} & +\frac{h}{2}\left[b\left(w_i\right)+b\left(w_i+a\left(t_i, w_i, \alpha\right) \Delta t_{i+1}\right)+b\left(t_i, w_i, \alpha\right) \Delta w_n\right] \Delta w_{i_1} \\ & +\frac{h}{2}\left[k\left(w_i\right)+k\left(w_i+a\left(t_i, w_i, \alpha\right) \Delta t_{i+1}\right)+b\left(t_i, w_i, \alpha\right) \Delta w_i\right] \Delta w_{i_2}\end{aligned}$                    (23)

where,

$\begin{gathered}X_0(\alpha)=\underline{X_0}(\alpha)+\frac{\overline{X_0}(\alpha)-\underline{X_0}(\alpha)}{S}, \\ w_0(\alpha)=\underline{w_0}(\alpha)+\frac{\overline{w_0}(\alpha)-\underline{w_0}(\alpha)}{S}, \\ w_{i+1}(\alpha)=\underline{w_{i+1}}(\alpha)+\frac{\overline{w_{i+1}}(\alpha)-\underline{w_{i+1}}(\alpha)}{S}, \\ \mathrm{a}\left(t_i, w_i, \alpha\right)=\underline{\mathrm{a}}\left(t_i, w_i, \alpha\right)+\frac{\overline{\mathrm{a}}\left(t_i, w_i, \alpha\right)-\underline{\mathrm{a}}\left(t_i, w_i, \alpha\right)}{S}, \\ \underline{\mathrm{~b}}\left(t_i, w_i, \alpha\right)+\frac{\overline{\mathrm{b}}\left(t_i, w_i, \alpha\right)-\underline{\mathrm{b}}\left(t_i, w_i, \alpha\right)}{S} \\ \underline{\mathrm{k}}\left(t_i, w_i, \alpha\right)+\frac{\overline{\mathrm{k}}\left(t_i, w_i, \alpha\right)-\underline{\mathrm{k}}\left(t_i, w_i, \alpha\right)}{S} .\end{gathered}$

Applying $\lim _{s \rightarrow \infty}$ and $\lim _{s \rightarrow 1}$ when we solve an equation, we get a left and a right term. We can also extract a number of different solutions by using different values of the membership function when $\alpha \in[0,1]$. Sometimes these solutions can be weak, leading to an overlap or intersection between the left and right terms due to the random nature of the system. This overlap can be clearly seen in the problems in the following examples.

3.6 Numerical stability analysis for Heun’s methods in fuzzy SODEs

We will discuss the two most common measures of stability: mean square and asymptotic [18-20]. Assuming that $X_0 \neq 0$ with probability 1, solutions of FSDE is:

$d \tilde{X}_t=a\left(t, \tilde{x}_t\right) d t+b\left(t, \tilde{x}_t\right) d w_1+k\left(t, \tilde{x}_t\right) d w_2$               (24)

Meeting the requirements of:

$\lim _{t \rightarrow \infty} E X^2(t)=0 \leftrightarrow R\{\lambda\}+\frac{1}{2}|\mu|^2<0$                (25)

$\lim _{t \rightarrow \infty}|X(t)|=0 \leftrightarrow R\left\{\lambda-\frac{1}{2} \mu^2\right\}<0$                 (26)

The left side of Eq. (25) defines the concept of mean-square stability, while the right side of Eq. (25) describes this property in detail using the FSDE function. Similarly, Eq. (26) defines and describes asymptotic stability.

In this study, a new numerical technique is developed for approximation of FSODEs with an integration of Heun’s method and dual Wiener processes. The major contributions include: Using triangular and trapezoidal fuzzy numbers to inject fuzziness, solving SODEs that cannot be solved to any productivity level using other approach. The method enhances both accuracy of solutions as well as the computational results and the potential has been established by convergence analysis and empirical evidence and it has been shown to be useful when solving more realistic problems such as fuzzy Langevin equations. This work fills voids that exist in earlier literature by presenting a feasible approach to model uncertainty in complicated stochastic structures, and has implications for practice in areas such as finance, engineering, and physics.

For possible future works, the method should be applied to higher-dimensional systems. Several nonlinearities including rational, periodic, and logarithmic nonlinearities should be included. In addition, such works may be considered to tackle complex nonlinear systems and integrating fuzzy logic with machine learning techniques, potentially artificial neural networks, to enhance solution accuracy.