Solvability of Fuzzy Partially Differentiable Models for Caputo-Hadamard Type Goursat Problems Involving Generalized Hukuhara Difference

1. Introduction

In the 17th century, integral and differential equations emerged as powerful mathematical tools. Over nearly four centuries of development, researchers have conducted increasingly in-depth investigations into these calculus-based equations. Emerging as a product of the development of differential equations, Goursat problem is a classic hyperbolic partial differential equation model as follows:

$$z_{xy}=f(x,y,z,z_x,z_y),\forall (x,y)\in P,$$

(1)

where z is a real-valued function defined on a rectangular domain P around the origin in ${\mathbb{R}}^2$, $z_{x,y}$ denotes mixed partial derivative of z with respect to x and y, and $z_x$ and $z_y$ are partial derivatives of z with respect to x and y, respectively. These derivatives are assumed to exist and be continuous. Additionally, the function $f(x,y,z,z_x,z_y)$ is bounded and continuous on $P\times {\mathbb{R}}^3$. We note that Goursat problem (1) plays an important role in describing wave phenomena such as sound and light waves [1].

It is well known that fuzzy set theory proposed by Zadeh [2] has been attracted widespread research attention. Various extensions, including intuitionistic, hesitant, and Pythagorean fuzzy sets, have been developed, contributing to the theory’s growing maturity. And then, comparing with partial differential equations established by Bers et al. [3] on classical sets in 1996, in order to provide more comprehensive and detailed representation of uncertain information systems. Buckley et al. [4] presented fuzzy partial differential equations, which also represent a new direction in the development of traditional partial differential equations including models of the form (1).

On the other hand, to address requirements of various environments, different types of fractional derivatives have been continuously developed to characterize the dynamic processes of diverse systems. Thereupon, fractional-order differential equations (FDEs) have gained prominence due to their suitability for addressing specific problem scenarios compared to their integer-order counterparts. Fractional calculus exhibits time memory and long-range spatial correlation, enabling more accurate modeling of physical and biochemical processes with memory, inheritance, and path dependence, as noted by many researchers. Despite these advantages, such properties are frequently overlooked in classical integer-order models. Consequently, significant role of FDEs has emerged as a crucial aspect of both mathematical theory and applied research. Furthermore, FDEs provide valuable insights for researchers across various fields. Their applicability has led to implementations in multiple domains, including biotechnology, physics, material mechanics, finance and so on. For more detail, one can refer to [5,6,7,8,9] and the references therein. It is worth noting that Hadamard introduced the kernel-dependent logarithmic form of Hadamard-type derivative integral in 1892, which marks a significant milestone in FDEs. Later, to tackle related issues in nonlinear systems, Italian mathematician Caputo proposed the concept of Caputo-type derivatives in 1939. In recent years, FDEs have been undergone rapid development, giving rise to novel derivatives such as Grunwald-Letnikov, Erdelyi-Kober, and Caputo-Hadamard types [10,11,12].

However, as research on FDEs has progressed, some researchers has observed that traditional FDEs frequently fail to adequately characterize complex systems that incorporate fuzzy information [13]. To address these drawbacks, fuzzy fractional-order differential equations (FFDEs) have been proposed as a promising solution. The concept of integrating fuzzy information with FDEs has its origins in research dating back to 1989, highlighting the longstanding interest in this area. A significant advancement on FFDE is the introduction of fuzzy Laplace transform [14], which has provided a more efficient method for solving Riemann-Liouville Hukuhara (H-) differential FFDEs in the form of $L\left[{(}^{RL}{D}_{a^{+}}^{\beta }f)\left(x\right)\right]$. Furthermore, based on analytical and approximate consistent solutions, Arqub et al. [15] established the compatibility of crack conformable FDE containing the following form:

$$D^{\beta }u\left(t\right)=\underset{\epsilon \to 0^{+}}{lim}{\displaystyle \frac{u\left(t\right)\ominus u{(t+\epsilon (t-t_0))}^{1-\beta }}{-\epsilon }},$$

where $\beta \in (0,1]$ is the order of generalized fuzzy conformable fractional derivative $D^{\beta }u$, and ⊖ represent H-difference. As all we know, regarding the numerical outcomes of fuzzy linear systems, it often demonstrates strong applicability and enhances both robustness and accuracy to utilize the tau method along with Jacobi polynomials for result simplification [16]. To better adapt to fuzzy environments, researchers have developed various innovative theories and methods for FFDEs. These include Laplace transform, Fourier transform, American decomposition method, Pythagorean fuzzy partial differential equation models, Legendre spectral method, and Brownian motion analysis. More detailed information of the about work can be found in [17,18,19,20,21].

Moreover, Aliev et al. [22] developed H-difference form of Z-numbers, and Amma et al. [23] demonstrated existence and uniqueness of solutions for a class of intuitionistic fuzzy set fractional partial differential equations (FFPDEs). It can not be ignored that traditional H-difference is often criticized for its limitations. Unlike the classical H-difference, the operational criteria and properties of the interval valued function space constructed based on generalized Hukuhara (gH-) type difference ${\ominus }_{gH}$ can be applied to the following semi-linear interval differential equations [24]:

$$\begin{array}{cc}& (f+g)\left(t\right)=f\left(t\right)+g\left(t\right),\hfill \\ & \left(kf\right)\left(t\right)=kf\left(t\right),\forall k\in \mathbb{R},\hfill \\ & \left(f{\ominus }_{gH}g\right)\left(t\right)=f\left(t\right){\ominus }_{gH}g\left(t\right),\hfill \\ & \left(fg\right)\left(t\right)=f\left(t\right)g\left(t\right).\hfill \end{array}$$

The application of suitable forms of gH-fractional derivatives in FFPDEs is also prevalent (see [25]). These advancements have substantially enriched the toolbox for modeling and solving complex dynamic systems. Specially, Zhang et al. [26] studied two gH-type weak solutions for the following Caputo fractional derivative couple system in a fuzzy gH-difference environment:

$$\left\{\begin{array}{c}{\phantom{2}}_{gH}^C{D}_k^{\alpha }u(x,y)=f(x,y,v(x,y)),\hfill \\ {\phantom{2}}_{gH}^C{D}_k^{\beta }v(x,y)=g(x,y,u(x,y)),\hfill \end{array}\right.$$

with initial condition

$$\left\{\begin{array}{ccc}\hfill & u(x,0)={\xi }_1\left(x\right),v(x,0)={\eta }_1\left(x\right),\hfill & \hfill \forall x\in [0,a],\\ \hfill & u(0,y)={\xi }_2\left(y\right),v(0,y)={\eta }_2\left(y\right),\hfill & \hfill \forall y\in [0,b],\end{array}\right.$$

where $k=1,2$, $\alpha =({\alpha }_1,{\alpha }_2),\beta =({\beta }_1,{\beta }_2)\in (0,1]\times (0,1]$ are fractional orders of and Caputo gH-type derivative operators ${\phantom{2}}_{gH}^C{D}_k^{\alpha }$ and ${\phantom{2}}_{gH}^C{D}_k^{\beta }$, respectively.

We notice that majority of researchers focus on extension of Caputo-type and Riemann-Liouville type derivatives. Notably, Caputo-Hadamard type derivative is regarded as an extension of Riemann-Liouville type derivative, it is because their integral kernel is in the form of $log\left({\displaystyle \frac{x}{t}}\right)$ rather than the traditional $(x-t)$ form. This modification in the integral kernel reflects a broader approach to capturing the dynamics of complex systems, thereby to expanding the applicability and utility of fractional calculus in various scientific fields. And then, Gohar et al. [27] obtained existence and uniqueness of solution for the following Caputo-Hadamard fractional differential system in a non fuzzified environment:

$$\left\{\begin{array}{c}{\phantom{2}}_{CH}{D}_{a^{+}}^{\alpha }u\left(t\right)=f(t,u),t>a^{+}>0,0<\alpha <1,\hfill \\ u\left(a^{+}\right)=u_a,\hfill \end{array}\right.$$

here ${\phantom{2}}_{CH}{D}_{a^{+}}^{\alpha }u$ is $\alpha$-order Caputo-Hadamard type derivative for a given function u, f is a continuous function and $u_a$ is initial value, which equivalent Volterra integral equation is:

$$u\left(t\right)=u_a+{\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_a^t{\left(log\left({\displaystyle \frac{t}{s}}\right)\right)}^{\alpha -1}f(s,u\left(s\right)){\displaystyle \frac{ds}{s}}.$$

Taking inspiration from the work of [26,27], combining with the work of [1], we consider the following fuzzy Goursat problem of Captou-Hadamard type fractional derivative with gH-difference:

$$\left\{\begin{array}{cc}\hfill & {\phantom{2}}_{CH}^{\theta }{D}_{xy}^{i,l}u(x,y)=a_1{\phantom{2}}_{CH}^{\theta }{D}_x^{i}u(x,y)+a_2{\phantom{2}}_{CH}^{\theta }{D}_y^{l}u(x,y)+a_{3}u(x,y)+G(x,y),\hfill \\ \hfill & u(x,y_0)={\varphi }_1\left(x\right),x_0\le x\le T,\hfill \\ \hfill & u(x_0,y)={\varphi }_2\left(y\right),y_0\le y\le T,\hfill \end{array}\right.$$

(2)

where $\theta \in (0,1]$ is the order of Captou-Hadamard gH-type (mixed) partial differential derivative defined and shown later in Definition 7 and Rmark 4, constant $a_{\iota }$ pertains to the domain of $\mathbb{R}$ for any $\iota =1,2,3$, $G:\mathbf{X}\to \mathbb{R}$ is a continuously differentiable function in the closed region $\mathbf{X}=[x_0,T]\times [y_0,T]$, and for $\kappa =1,2$, initial value function ${\varphi }_{\kappa }:\mathbf{X}\to \mathbb{E}$ is also continuously bounded on fuzzy space $\mathbb{E}$. For each $i,l\in \{1,2\}$ and any $u\in C_{i,l}(\mathbf{X},\mathbb{E})$, the space of continuous functions from $\mathbf{X}$ to $\mathbb{E}$, ${\phantom{2}}_{CH}^{\theta }{D}_{xy}^{i,l}u$ is Caputo-Hadamard gH-type mixed partial differential operator, and ${\phantom{2}}_{CH}^{\theta }{D}_x^i$ and ${\phantom{2}}_{CH}^{\theta }{D}_y^l$ separately represent Caputo-Hadamard gH-type partial differential operator in regard to variables x and y, here different values of i and l correspond to distinct spaces and the type († or ‡) of the partial differentiability with respect to x and y, respectively (for instance, ${\phantom{2}}_{CH}^{\theta }{D}_x^{1}u$ and ${\phantom{2}}_{CH}^{\theta }{D}_x^{2}u$ imply that u is †-type and ‡-type partially differentiable respect to x in several).

Remark 1. (i) If the order $\theta =1$, then (1) will be degenerated to the form of integer order. Thus, we extend (1) to (2) by generalizing the classical integer-order Goursat problem to a Goursat problem with Caputo-Hadamard fractional order derivative, here the integral kernel takes the form of a logarithmic function. This extension provides a feasible direction for future research on the chaotic dynamical behavior of the equation model (2), as well as circuit design and communication applications.

(ii) To accommodate Caputo-Hadamard type derivatives in the Goursat problem, we have refined the metric fixed-point method employed in [1].

(iii) Additionally, we have introduced a fuzzy environment (i.e. gH-difference) to enhance comprehensiveness of real-world problems that the system (2) can be modeled.

The rest of this paper is scheduled as follows: In Section 2, we shall review several fundamental properties and establish some utile results in this study. By employing Banach and Schauder fixed point theorems, existence and uniqueness of solutions for (2) is provided in Section 3. By using trapezoidal and Simpson’s rule [28], we illustrate numerical approximation form of the solutions of (2) via two examples in Section 4. Finally, the work conducted in this paper is summarized, and potential directions for future research are proposed in Section 5.

2. Preliminaries

In this section, in order to better handle (2), we will recall some concepts and properties of previous work on fuzzy fractional order Caputo-Hadmard type derivatives.

Throughout of this paper, let $\Xi$ represent classical fuzzy sets in fuzzy environments. Then, $\Xi$ is fuzzy convex, upper semi-continuous and compactly supported and indicate Fuzzy number. All spaces composed of $\Xi$ are called fuzzy spaces, which is represented by $\mathbb{E}$. In the same way, for membership function L, the following conditions are met:

(i) L is upper semicontinuous.

(ii) $L(\lambda x_1+(1-\lambda )x_2)\ge min\{L\left(x_1\right)$, $L\left(x_2\right)\}$, for each $x_1,x_2\in \mathbb{R}$ and any $\lambda \in [0,1]$.

(iii) $\exists x_0\in \mathbb{R}$ such that $L\left(x_0\right)=1$.

(iv) Closure of the set $\{x\in \mathbb{R}$: $L\left(x\right)>0\}$ is compact.

Definition 1

([29]). The highest measure $d_{\mathbb{E}}$ defined on $\mathbb{E}$ is shown as

$$\begin{array}{cc}\hfill d_{\mathbb{E}}(p,q)& =\underset{r\in [0,1]}{sup}{d}_H\left({\left[p\right]}^r,{\left[q\right]}^r\right)\hfill \\ \hfill & =\underset{r\in [0,1]}{sup}max\left\{|{\underline{p}}_r-{\underline{q}}_r|,|{\overline{p}}_r-{\overline{q}}_r|\right\},\forall p,q\in \mathbb{E},\hfill \end{array}$$

(3)

where ${\left[p\right]}^r$ represents the r-level set of fuzzy number p (the r-level set of q is also represented similarly) determined by

$${\left[p\right]}^r=\left\{\begin{array}{cc}\hfill & \{z\in \mathbb{R}:p(z)\ge r\},for0<r\le 1,\hfill \\ \hfill & cl\left(\mathrm{supp}p\right),forr=0,\hfill \end{array}\right.$$

(4)

${\underline{p}}_r$ represents the left endpoint and ${\overline{p}}_r$ represents the right endpoint of p. In addition, in (4), $cl$ is closure of the set and $\mathrm{supp}p$ denotes support set of p, which is defined by $\mathrm{supp}p=cl\{x\in \mathbb{R}:p(x)>0\}$.

In the sequel, let $len{\left[p\right]}^r={\overline{p}}_r-{\underline{p}}_r$ be diameter of the r-level set of p and $h(x;r)$ denote the r-level set of the function $h\left(x\right)$ for all $x\in \mathbb{R}$.

Lemma 1

([30]). The measure space $(\mathbb{E},d_{\mathbb{E}})$ satisfying the elimination rule $a+b=a+c\Rightarrow b=c$ for each $a,b,c\in (\mathbb{E},d_{\mathbb{E}})$, is complete and semilinear, one knows that for any $\varpi \in (\mathbb{E},d_{\mathbb{E}})$,

$\left(m_1\right)$$d_{\mathbb{E}}(a+b,c+\varpi )\le d_{\mathbb{E}}(a,c)+d_{\mathbb{E}}(b,\varpi )$;

$\left(m_2\right)$ if $a\ominus b,c\ominus \varpi$ exist, then $d_{\mathbb{E}}(a\ominus b,c\ominus \varpi )\le d_{\mathbb{E}}(a,c)+d_{\mathbb{E}}(b,\varpi )$;

$\left(m_3\right)$$a\ominus (b+c)=a\ominus b\ominus c$ when $a\ominus (b+c)$ and $a\ominus b\ominus c$ exist;

$\left(m_4\right)$ if $a\ominus (b\ominus c)$ and $a\ominus b$ both exist, $a\ominus (b\ominus c)=a\ominus b+c$;

$\left(m_5\right)$$(-1)(a\ominus b)=(-1)a\ominus (-1)b$ is valid if $a\ominus b$ and $(-1)a\ominus (-1)b$ exist.

Here, ⊖ represents H-difference defined by the following

$${\left[a\ominus b\right]}^r=\left[{\underline{a}}_r-{\underline{b}}_r,{\overline{a}}_r-{\underline{b}}_r\right]$$

Definition 2

([31]). $A_1,A_2,A_3\in \mathbb{E}$, gH-difference is defined as follows:

$$A_1{\ominus }_{gH}{A}_2=A_3⟺\left\{\begin{array}{cc}\hfill & (†)A_1=A_2+A_3,\hfill \\ \hfill & or(‡)A_2=A_1+(-1)A_3.\hfill \end{array}\right.$$

(5)

Remark 2

([31]). The r-level set of $A_1{\ominus }_{gH}{A}_2$ can be represented as ${\left[A_1{\ominus }_{gH}{A}_2\right]}^r=\left[\underline{L},\overline{\mathbb{R}}\right]$, where $\underline{{\rm Y}}=min\left\{\underline{A_1}\left(r\right)-\underline{A_2}\left(r\right),\overline{A_1}\left(r\right)-\overline{A_2}\left(r\right)\right\}$, $\overline{{\rm Y}}=max\left\{\underline{A_1}\left(r\right)-\underline{A_2}\left(r\right),\overline{A_1}\left(r\right)-\overline{A_2}\left(r\right)\right\}$.

Lemma 2

([31]). $h\left(x\right)$ is a fuzzy set-valued function defined on $[x_1,\infty ]$, and further defined in terms of $\left(\underline{h}(x;r),\overline{h}(x;r)\right)$. Here, $\underline{h}(x;r)$ and $\overline{h}(x;r)$ represent the left and right endpoint functions with respect to the level set parameter r, respectively. Moreover, for any point $r\in [0,1]$, both $\underline{h}(x;r)$ and $\overline{h}(x;r)$ are Riemann integrable. Further, for any x belonging to $[x_0,x_1]$, ${\int }_{x_0}^{x_1}\left|\underline{h}(x;r)\right|dx$ and ${\int }_{x_0}^{x_1}\left|\overline{h}(x;r)\right|$ are constrained by $\underline{K}\left(r\right)$ and $\overline{K}\left(r\right)$ i.e. ${\int }_{x_0}^{x_1}\left|\underline{h}(x;r)\right|dx\le$$\underline{K}\left(r\right)$ and ${\int }_{x_0}^{x_1}\left|\overline{h}(x;r)\right|$$\le \overline{K}\left(r\right)$. We can say that $h\left(x\right)$ is an abnormal fuzzy Riemannian integrable on $[x_1,\infty ]$, whose integral result is also a fuzzy number, with further results as follows:

$${\int }_{x_0}^{\infty }h(x;r)dx=\left[{\int }_{x_0}^{\infty }\underline{h}(x;r)dx,{\int }_{x_0}^{\infty }\overline{h}(x;r)dx\right].$$

Definition 3

([31]). For point $(x_0,y_0)\in (x_1,x_2)\times (y_1,y_2)\subset P$, here P is the same as in (1), $g:(x_1,x_2)\times (y_1,y_2)\to \mathbb{E}$ is said to be first-order generalized fuzzy partial differentiable to x, if the following relationship with respect to $g(x_0,y_0)$ exists:

(i)

$$\underset{\varrho \to 0}{lim}{\displaystyle \frac{g(x_0+\varrho ,y_0){\ominus }_{gH}g(x_0,y_0)}{\varrho }}=\underset{\varrho \to 0}{lim}{\displaystyle \frac{g(x_0,y_0){\ominus }_{gH}g(x_0-\varrho ,y_0)}{\varrho }}=g_x^{\prime }(x_0,y_0)$$

for a small enough given number $\varrho >0$ and a gH-difference both $g(x_0+\varrho ,y_0){\ominus }_{gH}g(x_0,y_0)$ and $g(x_0,y_0){\ominus }_{gH}g(x_0-\varrho ,y_0)$ existing in $(\mathbb{E},d_{\mathbb{E}})$.

(ii) Given a small enough number $\hslash >0$, if a gH-difference both $g(x_0,y_0){\ominus }_{gH}g(x_0+\varrho ,y_0)$ and $g(x_0-\hslash ,y_0){\ominus }_{gH}g(x_0,y_0)$ exist in $(\mathbb{E},d_{\mathbb{E}})$, then one has

$$\underset{\hslash \to 0}{lim}{\displaystyle \frac{g(x_0,y_0){\ominus }_{gH}g(x_0+\hslash ,y_0)}{-\hslash }}=\underset{\hslash \to 0}{lim}{\displaystyle \frac{g(x_0-\hslash ,y_0){\ominus }_{gH}g(x_0,y_0)}{-\hslash }}=g_x^{\prime }(x_0,y_0).$$

We call it †-type partial differentiable for the first case (i), and ‡-type partially differentiable in the second case (ii). Similar definitions can also be extended to the †- and ‡-types partial derivatives with respect to y, as well as to higher-order generalized fuzzy partial derivatives.

Remark 3

([31]). For $g:(x_1,x_2)\times (y_1,y_2)\to \mathbb{E}$ and level set parameter r, suppose that there exists $(x_0,y_0)\in$$(x_1,x_2)\times (y_1,y_2)$ such that $\overline{g}(·,·;r)$ and $\underline{g}(·,·;r)$ are both differentiable partially at $(x_0,y_0)$. Then

(i) the following relationship holds if g is †-type differentiable partially:

$$g_x^{\prime }(x,y)(x_0,y_0;r)=\left[\left({\underline{g}}_x^{\prime }\right)(x_0,y_0;r),\left({\overline{g}}_x^{\prime }\right)(x_0,y_0;r)\right],$$

(ii) and one has a correlation as follows

$$g_x^{\prime }(x,y)(x_0,y_0;r)=\left[\left({\overline{g}}_x^{\prime }\right)(x_0,y_0;r),\left({\underline{g}}_x^{\prime }\right)(x_0,y_0;r)\right].$$

when g is ‡-type differentiable partially.

For the partial derivative of g with respect to y, similar results also apply, which are not presented here.

Lemma 3

([31]). (Schauder fixed-point theorem) Λ is a bounded closed convex set that exists in a Banach space, and assuming that there is a self-mapping S that is completely continuous, we say that the mapping S has a fixed point in Λ.

Lemma 4

([34]). (Banach fixed-point theorem) If Ω is a complete metric space (i.e., every Cauchy sequence in Ω converges to a point within Ω) and f is a contraction mapping on Ω, then f has a unique fixed point $x^{*}\in \Omega$ such that $f\left(x^{*}\right)=x^{*}$.

Here, a contraction mapping refers to a function $f:\Omega \to \Omega$ such that there exists a constant $0\le \tau <1$ for which the following inequality holds for all $u,v\in \Omega$:

$$d\left(f\left(u\right),f\left(v\right)\right)\le \tau ·d(u,v),$$

where d denotes the metric on Ω.

Definition 4

([35]). Under the condition of order $\theta >0$, the fractional integral of Hadamard gH-type for the given function f is the fuzzy value function above $f\in {\mathcal{C}}^F[a,b]\cap {\mathcal{L}}^F[a,b]$ (f is a real-valued function defined on the interval $[a,b]$ that is both fuzzy continuous and fuzzy lebesgue integrable) and is defined as follows: for each $x>a>0$,

$${\phantom{2}}^H{I}_a^{\theta }f\left(x\right)={\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_a^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}f\left(s\right){\displaystyle \frac{ds}{s}},$$

and the following results apply to the r-level sets associated with integral of Hadamard gH-type:

(i) If f is †-type differentiable, it is evident that the associated r-level set satisfies the following relationship:

$${{[}^H{I}_a^{\theta }f\left(x\right)]}^r=\left[{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_a^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}\underline{f}(s;r){\displaystyle \frac{ds}{s}},{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_a^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}\overline{f}(s;r){\displaystyle \frac{ds}{s}}\right].$$

(ii) Otherwise, it is clear to see that the related r-level set fulfills

$${{[}^H{I}_a^{\theta }f\left(x\right)]}^r=\left[{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_a^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}\overline{f}(s;r){\displaystyle \frac{ds}{s}},{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_a^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}\underline{f}(s;r){\displaystyle \frac{ds}{s}}\right],$$

when $f\left(x\right)$ is ‡-type differentiable.

Definition 5

([35]). For all $x>a>0$, the fractional derivative of Hadamard gH-type for a given function f is fuzzy value function on ${\mathcal{C}}^F[a,b]\cap {\mathcal{L}}^F[a,b]$ and is defined by

$${\phantom{2}}^H{D}_{\phantom{2}}^{\theta }f\left(x\right)={\displaystyle \frac{1}{\Gamma (n-\theta )}}{\delta }^n{\int }_a^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{n-\theta -1}f\left(s\right){\displaystyle \frac{ds}{s}}$$

and the fractional order $\theta \in (n-1,n)$, where $n\in {\mathbb{Z}}^{+}$ and ${\delta }^n={\left(x{\displaystyle \frac{d}{dx}}\right)}^n$ is n-order derivative of Hadamard-type differentiation $x{\displaystyle \frac{d}{dx}}$.

Definition 6

([36]). Under the condition of order $\theta >0$, the fractional derivative of Caputo-Hadamard gH-type for a given function f on ${\mathcal{C}}^F[a,b]\cap {\mathcal{L}}^F[a,b]$ is defined as

$${\phantom{2}}^{CH}{D}_a^{\theta }f\left(x\right)={\displaystyle \frac{1}{\Gamma (n-\theta )}}{\int }_a^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{n-\theta -1}{\delta }^{n}f\left(s\right){\displaystyle \frac{ds}{s}}$$

and is fuzzy value for $x>a>0$ and $n-1<\theta <n\in {\mathbb{Z}}^{+}$, here ${\delta }^n$ is the same as Definition 5. Furthermore, when f is differentiable of †-type and ‡-type on ${\mathcal{C}}^F[a,b]\cap {\mathcal{L}}^F[a,b]$, f is said to be differentiable of †-type and ‡-type fractional Caputo-Hadamard gH-type derivative, respectively.

Definition 7.

Let ${\delta }^n$ be the same as in Definition 5. Under the condition of the order θ satisfying $n-1<\theta <n\in {\mathbb{Z}}^{+}$, Caputo-Hadamard gH-type partial differential operator with respect to x for a given function $g:(x_0,x_1)\times (y_0,y_1)\to \mathbb{E}$ is defined as

$${\phantom{2}}_{CH}^{\theta }{D}_x^{i}g(x,y)={\displaystyle \frac{1}{\Gamma (n-\theta )}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{s}}\right)}^{n-\theta -1}{\delta }^{n}g(s,y){\displaystyle \frac{ds}{s}}$$

and is fuzzy value for $x>x_0>0$, $y>y_0>0$ and $i\in \{1,2\}$. Furthermore,

(i) When $i=1$, that is, if $g(x,y)$ satisfies $\left(\mathrm{i}\right)$ of Definition 3, then ${\phantom{2}}_{CH}^{\theta }{D}_x^1$ is referred to as the type † Caputo-Hadamard gH-type partial differential operator with respect to x.

(ii) When $i=2$, that is, if $g(x,y)$ satisfies $\left(\mathrm{ii}\right)$ in Definition 3, then ${\phantom{2}}_{CH}^{\theta }{D}_x^2$ is referred to as the type ‡ Caputo-Hadamard gH-type partial differential operator with respect to x.

Remark 4.

Similar to Definition 7, Caputo-Hadamard gH-type partial differential operator with respect to y can be defined in a similar manner. Furthermore, for all $i,l\in \{1,2\}$, Caputo-Hadamard gH-type mixed partial differential operator of $g:(x_0,x_1)\times (y_0,y_1)\to \mathbb{E}$ is defined as

$${\phantom{2}}_{CH}^{\theta }{D}_{xy}^{i,l}g(x,y)={\displaystyle \frac{1}{\Gamma (n-\theta )}}{\displaystyle \frac{1}{\Gamma (n-\theta )}}{\int }_{x_0}^x{\int }_{y_0}^y{\left(log{\displaystyle \frac{x}{s}}\right)}^{n-\theta -1}{\left(log{\displaystyle \frac{y}{t}}\right)}^{n-\theta -1}{\tilde{\delta }}^{n}g(s,t){\displaystyle \frac{ds}{s}}{\displaystyle \frac{dt}{t}},$$

where $n-1<\theta <n\in {\mathbb{Z}}^{+}$ and ${\tilde{\delta }}^n={\left(xy{\displaystyle \frac{d^2}{dxdy}}\right)}^n$ is n-order derivative of mixed Hadamard-type differentiation $xy{\displaystyle \frac{d^2}{dxdy}}$, and for a fixed $l\in \{1,2\}$, ${\phantom{2}}_{CH}^{\theta }{D}_{xy}^{1,l}$${\phantom{2}}_{CH}^{\theta }{D}_{xy}^{2,l}$ separately represent †-type and ‡-type Caputo-Hadamard gH-type partial derivative operators with respect to x, and as far as a fixed $i\in \{1,2\}$, ${\phantom{2}}_{CH}^{\theta }{D}_{xy}^{i,1}$${\phantom{2}}_{CH}^{\theta }{D}_{xy}^{i,2}$ denote †-type and ‡-type Caputo-Hadamard gH-type partial derivative operators with respect to y, respectively.

Lemma 5.

Let $f\in {\mathcal{C}}^F[a,b]\cap {\mathcal{L}}^F[a,b]$, r be the level set parameter and the order $\theta >0$.

(I) If $len\left[f(x;r)\right]\le len\left[f(a;r)\right]$, then

($a_1$) $\left({\phantom{2}}^H{I}_a^{\theta }{\phantom{2}}^{CH}{D}_a^{\theta }f\right)\left(x\right)=f\left(x\right)$${\ominus }_{gH}{\sum }_{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}f\left(a\right)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı}$ if $f\left(x\right)$ is differentiable of †-type.

($a_2$) Otherwise, $\left({\phantom{2}}^H{I}_a^{\theta }{\phantom{2}}^{CH}{D}_a^{\theta }f\right)\left(x\right)=-f\left(x\right)$${\ominus }_{gH}(-1){\sum }_{ı=0}^{n-1}$${\displaystyle \frac{{\sigma }^{ı}f\left(a\right)}{ı!}}$${\left(log{\displaystyle \frac{t}{a}}\right)}^{ı}$ when $f\left(x\right)$ be differentiable of ‡-type.

(II) When $len\left[f(x;r)\right]\ge len\left[f(a;r)\right]$, one has

($b_1$) $\left({\phantom{2}}^H{I}_a^{\theta }{\phantom{2}}^{CH}{D}_a^{\theta }f\right)\left(x\right)=-f\left(x\right){\ominus }_{gH}$$(-1){\sum }_{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}f\left(a\right)}{ı!}}$${\left(log{\displaystyle \frac{t}{a}}\right)}^{ı}$ for taking $f\left(x\right)$ be differentiable of †-type.

($b_2$) $\left({\phantom{2}}^H{I}_a^{\theta }{\phantom{2}}^{CH}{D}_a^{\theta }f\right)\left(x\right)=f\left(x\right){\ominus }_{gH}{\sum }_{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}f\left(a\right)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı}$ for setting $f\left(x\right)$ be differentiable of ‡-type.

Proof. 

We only need to demonstrate the rationality of $len\left[f(x;r)\right]\le len\left[f(a;r)\right]$, and the corresponding $len\left[f(x;r)\right]\ge len\left[f(a;r)\right]$ can be similarly obtained. Since

$$\begin{array}{cc}\hfill & \left({\phantom{2}}^H{I}_a^{\theta CH}{D}_a^{\theta }\underline{f}\right)(x;r)=\underline{f}(x;r)-\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}\underline{f}(a;r)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı},\hfill \\ \hfill & \left({\phantom{2}}^H{I}_a^{\theta CH}{D}_a^{\theta }\overline{f}\right)(x;r)=\overline{f}(x;r)-\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}\overline{f}(a;r)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı},\hfill \end{array}$$

if it is †-type differentiable for f, by utilizing Lemma 2,we have

$$\begin{array}{cc}\hfill \left({\phantom{2}}^H{I}_a^{\theta CH}{D}_a^{\theta }f\right)(x;r)& =\left\{{\phantom{2}}^H{I}_a^{\theta CH}{D}_a^{\theta }\underline{f}(x;r){,}^H{I}_a^{\theta CH}{D}_a^{\theta }\overline{f}(a;r)\right\}\hfill \\ \hfill & =\{\underline{f}(x;r)-\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}\underline{f}(a;r)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı},\hfill \\ \hfill & \overline{f}(x;r)-\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}\overline{f}(a;r)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı}\}\hfill \\ \hfill & =f\left(x\right){\ominus }_{gH}\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}f\left(a\right)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı}.\hfill \end{array}$$

Further, one knows that

$$\begin{array}{cc}\hfill \left({\phantom{2}}^H{I}_a^{\theta CH}{D}_a^{\theta }f\right)(x;r)& =\left\{{\phantom{2}}^H{I}_a^{\theta CH}{D}_a^{\theta }\overline{f}(x;r){,}^H{I}_a^{\theta CH}{D}_a^{\theta }\underline{f}(a;r)\right\}\hfill \\ \hfill & =\{\overline{f}(x;r)-\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}\overline{f}(a;r)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı},\hfill \\ \hfill & \underline{f}(x;r)-\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}\underline{f}(a;r)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı}\}\hfill \\ \hfill & =-f\left(x\right){\ominus }_{gH}-\sum _{ı=0}^{n-1}{\displaystyle \frac{{\sigma }^{ı}f\left(a\right)}{ı!}}{\left(log{\displaystyle \frac{t}{a}}\right)}^{ı},\hfill \end{array}$$

when it is ‡-type differentiable for f. □

Lemma 6.

If $f(t,u)\in {\mathcal{C}}^F[a,b]\cap {\mathcal{L}}^F[a,b]$, then

(i) it follows from ($a_1$) and ($b_2$) of Lemma 5 that the following fuzzy initial value problem:

$$\left\{\begin{array}{cc}\hfill & {\phantom{2}}^{CH}{D}_a^{\theta }u\left(t\right)=f\left(t,u\left(t\right)\right),t>a>0,1>\theta >0,\hfill \\ \hfill & u\left(a\right)\in \mathbb{E}.\hfill \end{array}\right.$$

(6)

is equivalent to Volterra integral equation

$$u\left(t\right)=u\left(a\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_a^t{\left(log{\displaystyle \frac{t}{s}}\right)}^{\theta -1}f(s,u\left(s\right)){\displaystyle \frac{ds}{s}};$$

(ii) the equivalent Volterra integral form for (6) is derived as

$$u\left(t\right)=u\left(a\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_a^t{\left(log{\displaystyle \frac{t}{s}}\right)}^{\theta -1}f(s,u\left(s\right)){\displaystyle \frac{ds}{s}},$$

when $u\left(t\right)$ meets the conditions ($a_2$) and ($b_1$) of Lemma 5.

Proof. It follows from [[26] Lemma 3] and [[27] Lemma 2.6] that this conclusion can be obtained, therefore, it will not be discussed in detail here. □

3. Main Result

In this section, we consider existence and uniqueness of solutions to (2).

To obtain ideal results, we simplify and replace (2). Now, based on the same denotation for all symbols in (2) and using the equivalent Volterra integral equations in Lemma 6, we shall find the integral form of $u(x,y)\in C_{i,l}(\mathbf{X},\mathbb{E})\left(i,l\in \{1,2\}\right)$. Let ${\phantom{2}}_{CH}^{\theta }{D}_x^{i}u(x,y)=\nu (x,y)$, ${\phantom{2}}_{CH}^{\theta }{D}_y^{l}u(x,y)=\omega (x,y)$, and $\omega$, $\nu \in C(\mathbf{X},\mathbb{E})$ ( $C(\mathbf{X},\mathbb{E})$ represents the continuous bounded function space from $\mathbf{X}$ to $\mathbb{E}$ ), when $i,l=1$ it indicates that $u(x,y)$ satisfies case (i) of Lemma 6, and when $i,l=2$, $u(x,y)$ corresponds to case (ii) of Lemma 6. Then the fuzzy Goursat problem (2) can be transformed into the following equivalent form:

$$\left\{\begin{array}{c}{\phantom{2}}_{CH}^{\theta }{D}_x^i\omega (x,y)={\phantom{2}}_{CH}^{\theta }{D}_y^l\nu (x,y)=a_1\nu (x,y)+a_2\omega (x,y)+a_{3}u(x,y)+G(x,y)\hfill \\ \nu (x,y_0)={\phantom{2}}_{CH}^{\theta }{D}_x^i{\varphi }_1\left(x\right),x_0\le x\le T,\hfill \\ \omega (x_0,y)={\phantom{2}}_{CH}^{\theta }{D}_y^i{\varphi }_2\left(y\right),y_0\le y\le T.\hfill \end{array}\right.$$

(7)

Due to the different methods of taking i and l in (7), several different forms of solutions will be produced, and due to Lemmas 5 and 6 and (6), each type of solution also exhibits the following states.

The solutions of (7) can be earned in the form:

$$\left\{\begin{array}{c}u(x,y)={\varphi }_1\left(x\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \nu (x,y)={\phantom{2}}_{CH}^{\theta }{D}_x^1{\varphi }_1\left(x\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \omega (x,y)={\phantom{2}}_{CH}^{\theta }{D}_y^1{\varphi }_2\left(y\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x^{\prime },y){\displaystyle \frac{dx}{x^{\prime }}},\hfill \end{array}\right.$$

(8)

when $i=l=1$.

We can express the solutions of (7) as

$$\left\{\begin{array}{cc}\hfill & u(x,y)={\varphi }_1\left(x\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \hfill & \nu (x,y)={\phantom{2}}_{CH}^{\theta }{D}_x^1{\varphi }_1\left(x\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \hfill & \omega (x,y)={\phantom{2}}_{CH}^{\theta }{D}_y^2\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x^{\prime },y){\displaystyle \frac{dx}{x^{\prime }}},\hfill \end{array}\right.$$

(9)

if $i=1$ and $l=2$.

The solutions of (7) can be obtained in the form:

$$\left\{\begin{array}{cc}\hfill & u(x,y)={\varphi }_1\left(x\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \hfill & \nu (x,y)={\phantom{2}}_{CH}^{\theta }{D}_x^2{\varphi }_1\left(x\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \hfill & \omega (x,y)={\phantom{2}}_{CH}^{\theta }{D}_y^1{\varphi }_2\left(y\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x^{\prime },y){\displaystyle \frac{dx}{x^{\prime }}},\hfill \end{array}\right.$$

(10)

when $i=2$ and $l=1$.

One can obtain the solutions of (7) as follows:

$$\left\{\begin{array}{cc}\hfill & u(x,y)={\varphi }_1\left(x\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \hfill & \nu (x,y)={\phantom{2}}_{CH}^{\theta }{D}_x^2{\varphi }_1\left(x\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \hfill & \omega (x,y)={\phantom{2}}_{CH}^{\theta }{D}_y^2{\varphi }_2\left(y\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}(a_1\nu +a_2\omega +a_{3}u+G)(x^{\prime },y){\displaystyle \frac{dx}{x^{\prime }}},\hfill \end{array}\right.$$

(11)

if $i=l=2$.

Similar to [1], $u(x,y)$ satisfying integral equations (8) - (11) is the solution of (2). The current issue has been equivalently reformulated into an integral equation represented by $u(x,y)$, as presented in (8)–(11). Without loss of generality, the solution $u(x,y)$ for (7) exists in the following two equivalent integral forms

$$\begin{array}{cc}\hfill & †-\mathrm{type}\mathrm{solution}:u(x,y)={\varphi }_1\left(x\right)+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}},\hfill \\ \hfill & ‡-\mathrm{type}\mathrm{solution}:u(x,y)={\varphi }_1\left(x\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{dy}{y^{\prime }}}.\hfill \end{array}$$

(12)

In this section, we focus exclusively on the solvability of $‡-\mathrm{type}\mathrm{solution}$. By Definition 2 and (12), one can knows that $†-\mathrm{type}\mathrm{solution}$ can be viewed as a degenerate non-fuzzy form of $‡-\mathrm{type}\mathrm{solution}$, their solvability can be derived using similar methods and is therefore not elaborated in detail in this section.

Theorem 1.

If $u(x,y)\in C_{i,l}(\mathbf{X},\mathbb{E})$ in (12) satisfies (7), we say that it is the solution to the model (2).

Proof. 

Based on Lemma 1, we set up

$$\begin{array}{cc}& N=\left\{u\in C_{i,l}(\mathbf{X},\mathbb{E}):d_{\mathbb{E}}(u,{\varphi }_1\left(x\right))=\underset{r\in [0,1]}{sup}\underset{(x,y)\in \mathbf{X}}{max}\left\{|{\underline{u}}_r-\underline{{\varphi }_1}(x;r)|,\left|{\overline{u}}_r-\overline{{\varphi }_1}(x;r)\right|\right\}\le M_1\right\},\hfill \\ & \underset{u\in C_{i,l}(\mathbf{X},\mathbb{E}),\omega ,\nu \in C(\mathbf{X},\mathbb{E})}{max}\left\{(a_1\nu (x,y;r)+a_2\omega (x,y;r)+a_{3}u(x,y;r)+G(x,y;r))\in N\right\}\le M_2,\hfill \\ & \underset{(x,y)\in \mathbf{X}}{max}\left\{{\varphi }_1(x;r),{\varphi }_2(y;r)\right\}\le M_3,\hfill \end{array}$$

where $M_i\left(i\in \{1,2,3\}\right)$ is a constant. Narrow down the range of the set N as follows:

$$\begin{array}{cc}\hfill & N_h=\left\{u\in C_{i,l}(\mathbf{X},\mathbb{E}):d_{\mathbb{E}}(u_r,{\varphi }_1(x;r))=\underset{r\in [0,1]}{sup}\underset{(x,y)\in X_0}{max}\left\{|{\underline{u}}_r-\underline{{\varphi }_1}(x;r)|,\left|{\overline{u}}_r-\overline{{\varphi }_1}(x;r)\right|\right\}\le M_1\right\},\hfill \end{array}$$

(13)

where $X_0=[x_0,h]\times [y_0,h]$. If $x_0\ne y_0$, then

$$h\in \left[e^{{\displaystyle \frac{log{x}_0+log{y}_0-\sqrt{{(log{x}_0+log{y}_0)}^2-4\left(log{x}_{0}log{y}_0-A^{\frac{1}{\theta }}\right)}}{2}}},e^{{\displaystyle \frac{log{x}_0+log{y}_0+\sqrt{{(log{x}_0+log{y}_0)}^2-4\left(log{x}_{0}log{y}_0-A^{\frac{1}{\theta }}\right)}}{2}}}\right],$$

and $A=\left(A={\displaystyle \frac{\Gamma {(\theta +1)}^2}{M_2}}\left(M_1-{\displaystyle \frac{M_3}{\Gamma \left(\theta \right)}}\right)\right)$. If $x_0=y_0$, then $h=min\left\{x_0{e}^{\left[{\displaystyle \frac{M_1{\Gamma }^2(\theta +1)\Gamma \left(\theta \right)-{\Gamma }^2(\theta +1)M_3}{M_2\Gamma \left(\theta \right)}}\right]},T\right\}$. One can easily to see that $N_h$ is a smaller range of expressions within N and still maintains its original continuity, boundedness, and properties such as closed convex sets remain unchanged. In the following procedures, we exclusively employ left-endpoint functions from $N_h$, analogous conclusions hold for right-endpoint functions, which are not further elaborated here. The existence of $u(x,y)$ as ‡-type solution will be demonstrated next.

Define an operator $A:\mathbf{X}\to \mathbb{E}$ as follows:

$$\begin{array}{cc}\hfill \left(Au\right)(x,y)& ={\varphi }_1\left(x\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{\phantom{2}}_{CH}^{\theta }{D}_y^{l}u(x,y){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & ={\varphi }_1\left(x\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}.\hfill \end{array}$$

(14)

The problem (7) has been reformulated as demonstrating the existence of a fixed point for (14). We first prove that $A{N}_h\subseteq N_h$, by utilizing (5):

$$\begin{array}{cc}\hfill \left|\left(\underline{A}{u}_r\right)(x,y;r)-\underline{{\varphi }_1}(x;r)\right|& =\left|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|\hfill \\ \hfill & \le {\displaystyle \frac{M_3}{\Gamma \left(\theta \right)}}+{\displaystyle \frac{M_2}{\Gamma \left(\theta \right)\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\int }_{x_0}^x{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & ={\displaystyle \frac{M_3}{\Gamma \left(\theta \right)}}+{\displaystyle \frac{M_2}{\Gamma (\theta +1)\Gamma (\theta +1)}}{\left(log{\displaystyle \frac{y}{y_0}}\right)}^{\theta }{\left(log{\displaystyle \frac{x}{x_0}}\right)}^{\theta }.\hfill \end{array}$$

Further, one has

$$\begin{array}{c}\hfill d_{\mathbb{E}}(A\left(u\right)(x,y),{\varphi }_1\left(x\right))\le {\displaystyle \frac{M_3}{\Gamma \left(\theta \right)}}+{\displaystyle \frac{M_2}{\Gamma (\theta +1)\Gamma (\theta +1)}}{\left(log{\displaystyle \frac{y}{y_0}}\right)}^{\theta }{\left(log{\displaystyle \frac{x}{x_0}}\right)}^{\theta }\le M_1.\end{array}$$

(15)

From (15), we can immediately obtain $A{N}_h\subset N_h$.

Next, we shall prove the equicontinuity property of $A{N}_h$ and further demonstrate the existence of this solution. Take $u(x,y),u_n(x,y)\in N_h$ such that the following conditions are met:

$$\underset{n\to \infty }{lim}{d}_{\mathbb{E}}(u_n,u)=0.$$

Obviously, we have

$$\begin{array}{cc}\hfill \left|\left(\underline{A}{u}_n\right)(x,y;r)-\left(\underline{A}u\right)(x,y;r)\right|& ={\displaystyle \frac{1}{\Gamma \left(\theta \right)}}\left|{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\left({\underline{\omega }}_n(x,y;r)-\underline{\omega }(x,y;r)\right){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|\hfill \\ \hfill & \le {\displaystyle \frac{d_{\mathbb{E}}({\underline{\omega }}_n(x,y;r),\underline{\omega }(x,y;r))}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & ={\displaystyle \frac{d_{\mathbb{E}}({\underline{\omega }}_n(x,y;r),\underline{\omega }(x,y;r))}{\Gamma (\theta +1)}}{\left(log{\displaystyle \frac{y}{y_0}}\right)}^{\theta }.\hfill \end{array}$$

One can get

$$\begin{array}{c}\hfill d_{\mathbb{E}}(A{u}_n(x,y;r),Au(x,y;r))\le {\displaystyle \frac{d_{\mathbb{E}}({\underline{\omega }}_n(x,y;r),\underline{\omega }(x,y;r))}{\Gamma (\theta +1)}}{\left(log{\displaystyle \frac{y}{y_0}}\right)}^{\theta }.\end{array}$$

Due to the continuity of $\omega (x,y)$, we have

$$\underset{n\to \infty }{lim}{d}_{\mathbb{E}}({\omega }_n(x,y;r),\omega (x,y;r))=0.$$

Furthermore, we can obtain

$$\underset{n\to \infty }{lim}{d}_{\mathbb{E}}(A{u}_n(x,y;r),Au(x,y;r))=0.$$

The continuity of operator A is obtained from this, further proving its equicontinuity. Before further proof, we can easily conclude that the following conclusion holds true

$$\underset{(x,y)\to (x_0,y_0)}{lim}{\int }_{y_0}^y{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}=0.$$

(16)

From (16), one knows that there exist $h_1\in (x_0,h)$ and $h_2\in (y_0,h)$ such that the following results hold: for all $x\in (x_0,h_1)$ and any $y\in (y_0,h_2)$,

$${\displaystyle \frac{2M_2}{\Gamma \left(\theta \right)\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}<\epsilon .$$

Thus, when $y_0\le y_1<y_2\le h_2$ and $x_0\le x_1<x_2\le h_1$, we have

$$\begin{array}{cc}\hfill & \left|\underline{A}u(x_1,y_1;r)-\underline{A}u(x_1,y_2;r)\right|\hfill \\ \hfill & =\left|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_1,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}-{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_2}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_2,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|\hfill \\ \hfill & \le \left|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_1,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|+\left|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_2}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_2,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|\hfill \\ \hfill & \le {\displaystyle \frac{M_2}{\Gamma \left(\theta \right)\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}{\int }_{x_0}^{x_1}{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{x_1}{x^{\prime }}}\right)}^{\theta -1}{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & +{\displaystyle \frac{M_2}{\Gamma \left(\theta \right)\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_2}{\int }_{x_0}^{x_2}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{x_2}{x^{\prime }}}\right)}^{\theta -1}{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & \le {\displaystyle \frac{2M_2}{\Gamma \left(\theta \right)\Gamma \left(\theta \right)}}\epsilon \le \epsilon .\hfill \end{array}$$

(17)

To fully discuss the situation on the interval $[x_0,h]\times [y_0,h]$, we take $y_1,y_2\in \left[{\displaystyle \frac{y_0+h_2}{2}},h\right]$ and $x_1,x_2\in \left[{\displaystyle \frac{x_0+h_1}{2}},h\right]$, and have

$$\begin{array}{cc}\hfill & \left|\underline{A}u(x_1,y_1;r)-\underline{A}u(x_2,y_2;r)\right|\hfill \\ \hfill & =\left|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_1,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}-{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_2}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_2,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|\hfill \\ \hfill & =|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_1,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}-{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_2,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & +{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_2,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}-{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_2}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x_2,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}|\hfill \\ \hfill & \le {\displaystyle \frac{M_2}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}\left|\left[{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}-{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\right]\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}+{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_1}^{y_2}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\left|\underline{\omega }(x_2,y^{\prime };r)\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}},\hfill \end{array}$$

(18)

which first term implies that

$$\begin{array}{cc}\hfill & {\displaystyle \frac{M_2}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{y_1}\left|\left[{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}-{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\right]\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & \le {\displaystyle \frac{M_2}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{{\displaystyle \frac{y_0+h_2}{2}}}\left|\left[{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}-{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\right]\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & +{\displaystyle \frac{M_2}{\Gamma \left(\theta \right)}}{\int }_{{\displaystyle \frac{y_0+h_2}{2}}}^{y_1}\left|\left[{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}-{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\right]\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & \le {\displaystyle \frac{2M_2}{\Gamma \left(\theta \right)}}{\int }_{y_0}^{{\displaystyle \frac{y_0+h_2}{2}}}{\left(log{\displaystyle \frac{y_1}{y^{\prime }}}\right)}^{\theta -1}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & +{\displaystyle \frac{M_2}{\Gamma (\theta +1)}}\left[{\left(log{\displaystyle \frac{y_1}{\frac{a+h_2}{2}}}\right)}^{\theta }-{\left(log{\displaystyle \frac{y_2}{\frac{a+h_2}{2}}}\right)}^{\theta }+{\left(log{\displaystyle \frac{y_2}{y_1}}\right)}^{\theta }\right]\hfill \\ \hfill & \le 2\epsilon +{\displaystyle \frac{M_2}{\Gamma (\theta +1)}}\left[{\left(log{\displaystyle \frac{y_1}{\frac{a+h_2}{2}}}\right)}^{\theta }-{\left(log{\displaystyle \frac{y_2}{\frac{a+h_2}{2}}}\right)}^{\theta }+{\left(log{\displaystyle \frac{y_2}{y_1}}\right)}^{\theta }\right].\hfill \end{array}$$

(19)

Similarly, the following result holds for the second term on the right side of equation (18):

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_1}^{y_2}{\left(log{\displaystyle \frac{y_2}{y^{\prime }}}\right)}^{\theta -1}\left|\underline{\omega }(x_2,y^{\prime };r)\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}& \le {\displaystyle \frac{M_2}{\Gamma \left(\theta \right)}}{\int }_{y_1}^{y_2}{\left(log{\displaystyle \frac{y_2}{y_1}}\right)}^{\theta -1}{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & \le {\displaystyle \frac{M_2}{\Gamma (\theta +1)}}{\left(log{\displaystyle \frac{y_2}{y_1}}\right)}^{\theta }.\hfill \end{array}$$

(20)

Putting (19) and (20) into (18), then one gets

$$\begin{array}{cc}\hfill & \left|\underline{A}u(x_1,y_1;r)-\underline{A}u(x_2,y_2;r)\right|\hfill \\ \hfill & \le {\displaystyle \frac{M_2}{\Gamma (\theta +1)}}\left[{\left(log{\displaystyle \frac{y_1}{\frac{a+h_2}{2}}}\right)}^{\theta }-{\left(log{\displaystyle \frac{y_2}{\frac{a+h_2}{2}}}\right)}^{\theta }+2{\left(log{\displaystyle \frac{y_2}{y_1}}\right)}^{\theta }\right]+2\epsilon ,\hfill \end{array}$$

which means that there exist $h_1^{{\phantom{2}}^{\prime }}\in \left(0,{\displaystyle \frac{h_1-x_0}{2}}\right)$ and $h_2^{{\phantom{2}}^{\prime }}\in \left(0,{\displaystyle \frac{h_2-y_0}{2}}\right)$ such that for each $x_1,x_2\in \left[{\displaystyle \frac{x_0+h_1^{\prime }}{2}},h_1\right]$ and any $y_1,y_2\in \left[{\displaystyle \frac{y_0+h_2^{\prime }}{2}},h_2\right]$ satisfying $\left|y_1-y_2\right|<h_2$,

$$\left|\underline{A}u(x_1,y_1;r)-\underline{A}u(x_2,y_2;r)\right|<\epsilon .$$

(21)

Due to (17) and (21), it can be further concluded that $Au$ is equicontinuity. Furthermore, one can deduce that $Au$ is completely continuous. Hence, Lemma 3 ensures the existence of solutions to (2). □

In the sequence, we will use Lemma 4 to illustrate that (2) there exists only one ‡-type solution. Before proving, let us first list a hypothesis.

(H) 

Taking $\tilde{G}(x,y)$, $\tilde{u}(x,y)\in C_{i,l}(\tilde{\mathbf{X}},\mathbb{E})$, where $\tilde{\mathbf{X}}=[1,e]\times [1,e]$, we have $d_{\mathbb{E}}(G(x,y),\tilde{G}(x,y))\le d_{\mathbb{E}}(u(x,y),\tilde{u}(x,y))$.

Theorem 2.

Let $M_1$ be a constant as in the proof of Theorem 1. If the above assumption (H) is true and max{$|a_1|$, $|a_2|$, $|a_3|$, $M_1$}$\le 0.2$, then the differential equation system (2) has a unique solution $u(x,y)\in C_{i,l}(\tilde{\mathbf{X}},\mathbb{E})$.

Proof. 

Still define $N_h$ as (13) in Theorem 1, and $x_0\in [1,x]$, $y_0\in [1,y]$. The operator $\widehat{A}$ is now defined as

$$\begin{array}{c}\hfill \left(\widehat{A}u\right)(x,y)={\varphi }_1\left(x\right)\ominus (-1){\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\omega (x,y^{\prime }){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}.\end{array}$$

Now we have transformed the uniqueness problem of solution for (2) into a fixed point problem. The first step is to prove that for any $u\in N_h$,

$$\begin{array}{cc}\hfill \left|\left(\underline{\widehat{A}}u\right)(x,y;r)-\underline{{\varphi }_1}(x;r)\right|& =\left|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|\hfill \\ \hfill & \le {\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\left|\underline{\omega }(x,y^{\prime };r)\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & \le {\displaystyle \frac{M_2}{\Gamma (\theta +1)}}{\left(log{\displaystyle \frac{h}{y_0}}\right)}^{\theta }\hfill \\ \hfill & \le {\displaystyle \frac{M_2}{\Gamma (\theta +1)}}{\displaystyle \frac{M_1\Gamma (\theta +1)}{M_2}}=M_1,\hfill \end{array}$$

(22)

where $M_2$ is same as in the proof of Theorem 1.

$\widehat{A}u\in N_h$ can immediately be obtained by utilizing (22). From (3) and Lemma 1, we uncover the following findings: For $\tilde{u}\in N_h$ and $\tilde{\omega }={\phantom{2}}_{CH}^{\theta }{D}_y^l\tilde{u}$,

$$\begin{array}{cc}\hfill & d_{\mathbb{E}}(\widehat{A}u(x,y),\widehat{A}\tilde{u}(x,y))\hfill \\ \hfill & \le \left|{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\underline{\omega }(x,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}-{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\underline{\tilde{\omega }}(x,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\right|\hfill \\ \hfill & \le {\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}\left|\underline{\omega }(x,y^{\prime };r)-\underline{\tilde{\omega }}(x,y^{\prime };r)\right|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}},\hfill \end{array}$$

and through the utilization of (11), one gets

$$\begin{array}{cc}\hfill |\underline{\omega }& (x,y^{\prime };r)-\underline{\tilde{\omega }}(x,y^{\prime };r)|\hfill \\ \hfill =& |{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}(a_1\underline{\nu }+a_2\underline{\omega }+a_3\underline{u}+\underline{G})(x^{\prime },y^{\prime };r){\displaystyle \frac{dx}{x^{\prime }}}\hfill \\ \hfill & -{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}(a_1\tilde{\underline{\nu }}+a_2\tilde{\underline{\omega }}+a_3\tilde{\underline{u}}+\tilde{\underline{G}})(x^{\prime },y^{\prime };r){\displaystyle \frac{dx}{x^{\prime }}}|\hfill \\ \hfill \le & {\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}\left(a_1\left|\underline{\nu }-\tilde{\underline{\nu }}\right|+a_2\left|\underline{\omega }-\tilde{\underline{\omega }}\right|+a_3\left|\underline{u}-\tilde{\underline{u}}\right|+\left|\underline{G}-\tilde{\underline{G}}\right|\right)(x^{\prime },y^{\prime };r){\displaystyle \frac{dx}{x^{\prime }}}\hfill \\ \hfill =& {\int }_{x_0}^x{a}_1{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}\left|\underline{\nu }(x^{\prime },y^{\prime };r)-\underline{\tilde{\nu }}(x^{\prime },y^{\prime };r)\right|{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}+{\int }_{x_0}^x{a}_2{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}\left|\underline{\omega }(x^{\prime },y^{\prime };r)-\underline{\tilde{\omega }}(x^{\prime },y^{\prime };r)\right|{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}\hfill \\ \hfill & +{\int }_{x_0}^x{a}_3{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}\left|\underline{u}(x^{\prime },y^{\prime };r)-\underline{\tilde{u}}(x^{\prime },y^{\prime };r)\right|{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}}+{\int }_{x_0}^x{\left(log{\displaystyle \frac{x}{x^{\prime }}}\right)}^{\theta -1}\left|\underline{G}(x^{\prime },y^{\prime };r)-\underline{\tilde{G}}(x^{\prime },y^{\prime };r)\right|{\displaystyle \frac{d{x}^{\prime }}{x^{\prime }}},\hfill \end{array}$$

here $\tilde{\nu }={\phantom{2}}_{CH}^{\theta }{D}_x^i\tilde{u}$ and

$$\begin{array}{cc}\hfill |\underline{\nu }& (x,y;r)-\underline{\tilde{\nu }}(x,y;r)|\hfill \\ \hfill \le & {\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_{y_0}^y{\left(log{\displaystyle \frac{y^{\prime }}{y}}\right)}^{\theta -1}\left(a_1|\underline{\nu }-\underline{\tilde{\nu }}|+a_2|\underline{\omega }-\underline{\tilde{\omega }}|+a_3|\underline{\mu }-\underline{\tilde{\mu }}|+|\underline{G}-\underline{\tilde{G}}|\right)(x,y^{\prime };r){\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill \le & {\displaystyle \frac{1}{\Gamma \left(\theta \right)}}\left[{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{a}_1\right|\underline{\nu }(x,y;r)-\underline{\tilde{\nu }}(x,y;r)|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}+{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{a}_2|\underline{\omega }(x,y;r)-\underline{\tilde{\omega }}(x,y;r)|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}\hfill \\ \hfill & +{\int }_{y_0}^y{\left(log{\displaystyle \frac{y}{y^{\prime }}}\right)}^{\theta -1}{a}_3|\underline{u}(x,y;r)-\underline{\tilde{u}}(x,y;r)|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}+{\int }_a^y|\underline{G}(x,y;r)-\underline{\tilde{G}}(x,y;r)|{\displaystyle \frac{d{y}^{\prime }}{y^{\prime }}}].\hfill \end{array}$$

Thus, it follows that

$$\begin{array}{cc}\hfill & d_{\mathbb{E}}(\nu (x,y;r),\tilde{\nu }(x,y;r)){\left(\Gamma \left(\theta +1\right)-\left|a_1\right|log{\displaystyle \frac{h}{y_0}}\right)}^{\theta }\hfill \\ \hfill & \le d_{\mathbb{E}}(\underline{\omega }(x,y;r),\underline{\tilde{\omega }}(x,y;r))|a_2|{\left(log{\displaystyle \frac{h}{y_0}}\right)}^{\theta }+d_{\mathbb{E}}(\underline{u}(x,y;r),\underline{\tilde{u}}(x,y;r))\left|a_3\right|{\left(log{\displaystyle \frac{h}{y_0}}\right)}^{\theta }+M_1{\left(log{\displaystyle \frac{h}{y_0}}\right)}^{\theta }.\hfill \end{array}$$

(23)

Next, we get

$$\begin{array}{cc}\hfill d_{\mathbb{E}}& (\omega (x,y;r),\tilde{\omega }(x,y;r))\hfill \\ \hfill \le & {\displaystyle \frac{d_{\mathbb{E}}(\underline{\omega }(x,y;r),\underline{\tilde{\omega }}(x,y;r))\left|a_2\right|{(log\frac{h}{x_0})}^{\theta }\Gamma (\theta +1)}{\Gamma (\theta +1)-|a_1|{(log\frac{h}{x_0})}^{\theta }}}+{\displaystyle \frac{d_{\mathbb{E}}(\underline{u}(x,y;r),\underline{\tilde{u}}(x,y;r))\left(\right|a_3\left|\right|a_1|{(log\frac{h}{x_0})}^{2\theta }}{\Gamma (\theta +1)-|a_1|{(log\frac{h}{x_0})}^{\theta }}}\hfill \\ \hfill & +{\displaystyle \frac{|a_3|{(log\frac{h}{x_0})}^{\theta }\Gamma (\theta +1)-|a_1{a}_3\left|{(log\frac{h}{x_0})}^{2\theta }\right)}{\Gamma (\theta +1)-|a_1|{(log\frac{h}{x_0})}^{\theta }}}\hfill \\ \hfill \le & \left[{\displaystyle \frac{|a_1|d_{\mathbb{E}}(\underline{u}(x,y;r)-\tilde{\underline{u}}(x,y;r))+\left|a_1\right|M_1}{{\Gamma }^2(\theta +1)-|a_1|-|a_2|\Gamma (\theta +1)}}+{\displaystyle \frac{|a_3|\Gamma (\theta +1)d_{\mathbb{E}}(\underline{u}(x,y;r)-\tilde{\underline{u}}(x,y;r))+M_1\Gamma (\theta +1)}{{\Gamma }^2(\theta +1)-|a_1|-|a_2|\Gamma (\theta +1)}}\right].\hfill \end{array}$$

(24)

Ultimately, combining (23) and (24), we achieve the outcome

$$\begin{array}{cc}\hfill & d_{\mathbb{E}}(\widehat{A}u(x,y),\widehat{A}\tilde{u}(x,y))\hfill \\ \hfill & \le \left[{\displaystyle \frac{|a_3|{\left(log\frac{h}{y_0}\right)}^{\theta }{d}_{\mathbb{E}}(\underline{u}(x,y;r),\underline{\tilde{u}}(x,y;r))+M_1{\left(log\frac{h}{x_0}\right)}^{\theta }{\left(log\frac{h}{y_0}\right)}^{\theta }}{{\Gamma }^2(\theta +1)-\left(|a_1|{\left(log\frac{h}{y_0}\right)}^{\theta }\right)\Gamma (\theta +1)-\left(|a_2|{\left(log\frac{h}{y_0}\right)}^{\theta }\right)\Gamma (\theta +1)}}\right]\hfill \\ \hfill & \le K{d}_{\mathbb{E}}(u(x,y),\tilde{u}(x,y)),\hfill \end{array}$$

where

$$K=\left[{\displaystyle \frac{|a_3|{\left(log\frac{h}{y_0}\right)}^{\theta }+M_1{\left(log\frac{h}{x_0}\right)}^{\theta }{\left(log\frac{h}{y_0}\right)}^{\theta }}{{\Gamma }^2(\theta +1)-\left(|a_1|{\left(log\frac{h}{y_0}\right)}^{\theta }\right)\Gamma (\theta +1)-\left(|a_2|{\left(log\frac{h}{y_0}\right)}^{\theta }\right)\Gamma (\theta +1)}}\right].$$

Hence, drawing from (H) and max{$|a_1|$, $|a_2|$, $|a_3|$, $M_1$}$\le 0.2$, one can conclude that $K<1$ and knows that there exists a unique solution $u(x,y)\in C_{i,l}(\tilde{\mathbf{X}},\mathbb{E})$ for $i,l\in \{1,2\}$ by Lemma 4. □

4. Numerical Examples

In this section, we present two specific examples and will obtain numerical approximations of their solutions by using methods due to Dhali et at. [28]. Additionally, we provide 3D visualizations of the approximate solutions.

Definition 8. ([28]) Assuming $x\in [x_0,x_m]$, and $x_0<x_1<x_2<\dots <x_{m-1}<x_m$, then trapezoidal formula for numerical integration in uneven spaces is given by:

$$\begin{array}{cc}\hfill {\int }_{x_0}^{x_m}f\left(x\right)dx\approx & {\displaystyle \frac{1}{2}}[(x_1-x_0)(f\left(x_0\right)+f\left(x_1\right))+(x_2-x_1)(f\left(x_1\right)+f\left(x_2\right))\hfill \\ \hfill & +\dots +(x_{m-1}-x_{m-2})(f\left(x_{m-1}\right)+f\left(x_{m-2}\right))\hfill \\ \hfill & +(x_m-x_{m-1})(f\left(x_{m-1}\right)+f\left(x_m\right))].\hfill \end{array}$$

Definition 9. ([28]) Assuming $x\in [x_1,x_2]$ and $y\in [y_1,y_2]$, the following numerical integration formula (also known as Simpson’s rule) is applied:

$$\begin{array}{cc}\hfill & {\int }_{y_1}^{y_2}{\int }_{x_1}^{x_2}f(x,y)dxdy\hfill \\ \hfill & \approx {\displaystyle \frac{(x_2-x_1)}{6}}{\displaystyle \frac{(y_2-y_1)}{6}}(f(x_1,y_1)+f(x_2,y_1)+f(x_1,y_2)+f(x_2,y_2)\hfill \\ \hfill & +4[f\left({\displaystyle \frac{x_1+x_2}{2}},y_1\right)+f\left(x_1,{\displaystyle \frac{y_1+y_2}{2}}\right)+f\left(x_2,{\displaystyle \frac{y_1+y_2}{2}}\right)\hfill \\ \hfill & +f\left({\displaystyle \frac{x_1+x_2}{2}},y_2\right)]+16f\left({\displaystyle \frac{x_1+x_2}{2}},{\displaystyle \frac{y_1+y_2}{2}}\right)).\hfill \end{array}$$

(25)

Example 1.

Consider the following non abstract Caputo-Hadamard type fractional derivative equation for $0<\theta \le 1$:

$$\left\{\begin{array}{cc}\hfill & {\phantom{2}}_{CH}^{\theta }{D}_{xy}^{1,1}u(x,y)={\displaystyle \frac{\sqrt{2}}{7}}log\left(x\right)log\left(y\right),\hfill \\ \hfill & u(x,y_0)={\varphi }_1\left(x\right)=\sqrt{2},1⩽x⩽e,\hfill \\ \hfill & u(x_0,y)={\varphi }_2\left(y\right)=\sqrt{3},1⩽y⩽e.\hfill \end{array}\right.$$

(26)

For initial conditions $u(x,y_0)$ and $u(x_0,y)$, and the right-hand side of the equation of (26), it is clear to conclude that (26) exist solutions $u(x,y)\in C_{1,1}(\tilde{\mathbf{X}},\mathbb{E})$, here $\tilde{\mathbf{X}}=[1,e]\times [1,e]$. The following numerical method ensures the uniqueness of the solution to (26). We shall obtain solutions for Example 1 at different fractional orders based on the trapezoidal formula in Definition 8. We divide the interval $[1,e]$ for x and the interval $[1,e]$ for y into equally spaced subintervals, and collect the division points into sets $X=\{x_1,x_2,\dots ,x_m\}$ and $Y=\{y_1,y_2,\dots ,y_m\}$, respectively, where $1<x_1<x_2<\dots <x_m<e$, $1<y_1<y_2<\dots <y_m<e$ and $x_{\ell +1}-x_{\ell }=y_{j+1}-y_j={\displaystyle \frac{e-1}{m}}\left(\ell ,j=1,2,\dots ,m\right)$. Utilizing Definition 8, for any $x,y\in [1,e]$, one obtains that

$$\begin{array}{cc}\hfill & {\int }_1^x{\int }_1^y{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)log\left(t\right)}{st}}dsdt\hfill \\ \hfill & \approx \sum _{\ell =2}^{m}log{\left({\displaystyle \frac{x_{\ell -1}}{x_{\ell }}}\right)}^{\theta -1}{\displaystyle \frac{log\left(x_{\ell }\right)}{x_{\ell }}}\sum _{j=2}^{m}log{\left({\displaystyle \frac{y_{j-1}}{y_j}}\right)}^{\theta -1}{\displaystyle \frac{log\left(y_j\right)}{y_j}}.\hfill \end{array}$$

Then one can theoretically get solutions of (26) as

$$\begin{array}{cc}\hfill u(x,y)& =\sqrt{2}+\sqrt{3}+{\displaystyle \frac{\sqrt{2}}{7}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_1^x{\int }_1^y{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)log\left(t\right)}{st}}dsdt.\hfill \end{array}$$

(27)

For $\ell ,j=1,2,\dots ,m$, we have

$$\begin{array}{cc}\hfill u(x_{\ell },y_j)& \approx \sqrt{2}+\sqrt{3}+{\displaystyle \frac{\sqrt{2}}{7}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}\sum _{\ell =2}^{m}log{\left({\displaystyle \frac{x_{\ell -1}}{x_{\ell }}}\right)}^{\theta -1}{\displaystyle \frac{log\left(x_{\ell }\right)}{x_{\ell }}}\sum _{j=2}^{m}log{\left({\displaystyle \frac{y_{j-1}}{y_j}}\right)}^{\theta -1}{\displaystyle \frac{log\left(y_j\right)}{y_j}}.\hfill \end{array}$$

Take $m=100$ and the order $\theta =0.3,0.4,0.5,1$ ($\theta =1$ is the case of integer-order for (2)), respectively. For the †-type solution of different orders in Example 1, numerical integration is utilized to obtain approximate solutions, and 3D visualizations of these approximations are presented (see Figure 1, Figure 2, Figure 3 and Figure 4). Moreover, it can be observed that when $\theta$ is a fractional order ($\theta =0.3,0.4,0.5$), $u(x,y)$ exhibits a broader coverage compared to the integer-order case (i.e., $\theta =1$).

Example 2.

For $i=l=2$, the following Caputo-Hadamard fractional differential equation with triangular fuzzy number $M=(1,2,3)$ as the initial value is viewed:

$$\left\{\begin{array}{cc}\hfill & {\phantom{2}}_{CH}^{\theta }{D}_{xy}^{2,2}u(x,y)={\displaystyle \frac{\sqrt{7}}{15}}[log\left(x\right)+log\left(y\right)],\hfill \\ \hfill & u(x,y_0)={\varphi }_1\left(x\right)=-M,1⩽x⩽e,\hfill \\ \hfill & u(x_0,y)={\varphi }_2\left(y\right)=4M,1⩽y⩽e.\hfill \end{array}\right.$$

(28)

Similarly, we can easily obtain that the system of equations satisfies relevant theorems, so that (28) has solutions $u(x,y)\in C_{2,2}(\tilde{\mathbf{X}},\mathbb{E})$. Similar to Example 1, we utilize a numerical method to demonstrate that (28) has a unique solution and leting $M=(1,2,3)$, then one can easily obtain

$${\left[M\right]}^r=\left[{\underline{M}}_r,{\overline{M}}_r\right]=[r+1,3-r],{\left[4M\right]}^r=\left[4{\underline{M}}_r,4{\overline{M}}_r\right]=[4r+4,12-4r].$$

Based on Lemma 6 provided earlier, similar to Example 1, we can obtain fuzzy solution:

$$\begin{array}{cc}\hfill u(x,y)& =u(x,0)+u(0,y)\ominus (-1){\displaystyle \frac{15}{\sqrt{7}}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_1^x{\int }_1^y{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)+log\left(t\right)}{st}}dsdt.\hfill \end{array}$$

Next, we will characterize existence and uniqueness of solution to (28) based on r level set of fuzzy solutions:

$$\begin{array}{cc}\hfill & {\left[u(x,y)\right]}^r=\hfill \\ \hfill & \left[3r+3+{\displaystyle \frac{15}{\sqrt{7}}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_1^x{\int }_1^y{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)+log\left(t\right)}{st}}dsdt,\right.\hfill \\ \hfill & \left.9-3r+{\displaystyle \frac{15}{\sqrt{7}}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\displaystyle \frac{1}{\Gamma \left(\theta \right)}}{\int }_1^x{\int }_1^y{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)+log\left(t\right)}{st}}dsdt\right].\hfill \end{array}$$

By utilizing the trapezoidal method for numerical integration defined in Definition 9, We achieve

$$\begin{array}{cc}\hfill & {\int }_1^x{\int }_1^y{\left(log{\displaystyle \frac{x}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)+log\left(t\right)}{st}}dsdt\hfill \\ \hfill & \approx \sum _{\ell =1}^m\sum _{j=1}^m{\int }_{x_{\ell }}^{x_{\ell +1}}{\int }_{y_j}^{y_{j+1}}{\left(log{\displaystyle \frac{x_{\ell +1}}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y_{j+1}}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)+log\left(t\right)}{st}}dsdt,\hfill \end{array}$$

and for $\ell ,j=1,2,\dots ,m$,

$$f(s,t)={\left(log{\displaystyle \frac{x_{\ell +1}}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y_{j+1}}{l}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)+log\left(t\right)}{st}}.$$

According to (25), one can get for $\ell ,j=1,2,\dots ,m$,

$$\begin{array}{cc}\hfill \sum _{\ell =1}^m\sum _{j=1}^m& {\int }_{x_{\ell }}^{x_{\ell +1}}{\int }_{y_j}^{y_{j+1}}{\left(log{\displaystyle \frac{x_{\ell +1}}{s}}\right)}^{\theta -1}{\left(log{\displaystyle \frac{y_{j+1}}{t}}\right)}^{\theta -1}{\displaystyle \frac{log\left(s\right)+log\left(t\right)}{st}}dsdt\hfill \\ \hfill \approx & \sum _{\ell =1}^m\sum _{j=1}^m{\displaystyle \frac{(x_{\ell +1}-x_{\ell })(y_{j+1}-y_j)}{6}}[f(x_{\ell },y_j)+f(x_{\ell },y_{j+1})+f(x_{\ell +1},y_j)\hfill \\ \hfill & +f(x_{\ell +1},y_{j+1})]+4[f\left({\displaystyle \frac{x_{\ell }+x_{\ell +1}}{2}},y_j\right)+f\left({\displaystyle \frac{x_{\ell }+x_{\ell +1}}{2}},y_{j+1}\right)\hfill \\ \hfill & +f\left(x_{\ell },{\displaystyle \frac{y_j+y_{j+1}}{2}}\right)+f\left(x_{\ell +1},{\displaystyle \frac{y_j+y_{j+1}}{2}}\right)]+16\left[f\left({\displaystyle \frac{x_{\ell }+x_{\ell +1}}{2}},{\displaystyle \frac{y_j+y_{j+1}}{2}}\right)\right].\hfill \end{array}$$

Similar to Example 1, we divide the interval $[1,e]$ in the same manner, with the set of partition points denoted as $X=\{x_1,x_2,\dots ,x_m\}$ and $Y=\{y_1,y_2,\dots ,y_m\}$, and the partition distance as $x_{\ell +1}-x_{\ell }=y_{j+1}-y_j={\displaystyle \frac{e-1}{m}}\left(\ell ,j=1,2,3,\dots ,m\right)$. Based on Definition 9 and taking $m=100$, numerical integration is applied at each mesh point to obtain approximations. Figure 5, Figure 6, Figure 7 and Figure 8 with $\theta =0.4,0.5$ as $r=0.5$ and $r=0.1,0.5$ as $\theta =0.3$, respectively. Furthermore, 3D visualizations were generated for different orders on the same level sets (see Figure 5 and Figure 6), as well as for the same order on different level sets (see Figure 7 and Figure 8).

5. Concluding Remarks

The primary purpose of this paper is to focus on existence and uniqueness of solutions for the following FFPDE of Caputo-Hadamard type involving gH-difference:

$$\left\{\begin{array}{cc}\hfill & {\phantom{2}}_{CH}^{\theta }{D}_{xy}^{i,l}u(x,y)=a_1{\phantom{2}}_{CH}^{\theta }{D}_x^{i}u(x,y)+a_2{\phantom{2}}_{CH}^{\theta }{D}_y^{l}u(x,y)+a_{3}u(x,y)+G(x,y),\hfill \\ \hfill & u(x,y_0)={\varphi }_1\left(x\right),x_0\le x\le T,\hfill \\ \hfill & u(x_0,y)={\varphi }_2\left(y\right),y_0\le y\le T,\hfill \end{array}\right.$$

(29)

under some suitable conditions. The classic Caputo FDE had been innovatively extended to Caputo-Hadamard FDE, which held significant meaning due to its integral kernel being in the form of a logarithm rather than the traditional sense of subtraction.

Building upon the work of previous researchers, the main innovations and objectives of this paper are as follows:

(i) In the context of Goursat problems, we extended traditional Caputo FDEs to Caputo-Hadamard FDEs, thereby altering the form of the integral kernel.

(ii) We established the uniqueness and existence of solution to problem (29) using Banach fixed point theorem and Schauder fixed point theorem, rather than relying on metric fixed point theory for verification.

(iii) We presented relevant numerical examples applicable to existence and uniqueness theorem, and further obtained an approximate numerical solution for (29) by discretizing $u(x,y)$.

In disciplines such as finance, operations research, optics, medicine, biology, engineering, and other practical applications, time delays and pulse effects often arise due to instrument measurements or human interventions. A valuable potential direction for future work is to investigate whether the solution of (29) maintains existence and uniqueness under time delay $(x-t,y-\sigma )$ and during the occurrence of the pulse phenomenon.