Properties of q -Differential Equations of Higher Order and Visualization of Fractal Using q -Bernoulli Polynomials

: We introduce several q -differential equations of higher order which are related to q -Bernoulli polynomials and obtain a symmetric property of q -differential equations of higher order in this paper. By giving q -varying variations, we identify the shape of the approximate roots of q -Bernoulli polynomials, a solution of q -differential equations of higher order, and ﬁnd several conjectures associated with them. Furthermore, based on q -Bernoulli polynomials, we create a Mandelbrot set and a Julia set to ﬁnd a variety of related ﬁgures.


Introduction
One of the differential equations that converts nonlinear equations into linear equations is the Bernoulli differential equation. A Bernoulli differential equation is an equation of the form dy dx + p(x)y − g(x)y m = 0, where m is any real number, p(x) and g(x) are continuous functions on the interval. If m = 0 or m = 1, the above equation is linear, and if not, the equation is nonlinear. The Bernoulli differential equation can be reduced to a linear differential equation with substitution u = y 1−m . Then, for u we obtain a linear equation du dx + (1 − m)p(x)u = (1 − m)g(x). This Bernoulli differential equation has many applications to problems modeled by nonlinear differential equations, equations about the population expressed in logistic equations or Verhulst equations, physics, etc.
If m = 0 in (1), then the Bernoulli differential equation has the solution which is the generating function of the Bernoulli polynomials. The equation is as follows.
where B n (x) is the Bernoulli polynomials, see [1,2]. For |t| < 2π, the Bernoulli numbers and polynomials can be expressed as ∞ ∑ n=0 B n t n n! = t e t − 1 , ∞ ∑ n=0 B n (x) t n n! = t e t − 1 e tx , respectively.
Based on the concept above, we can consider the q-Bernoulli differential equation of the first order D q y + p(x)y − g(x)y m = 0 in q-calculus. When m = 0 in (1), the q-Bernoulli polynomials are a solution of the following q-differential equation of the first order. D q B n,q (x) + q(B 0,q − 1) B 1,q B n,q (x) + [n] q (q −n B n,q (qx) − (x − q −2 )B n−1,q (x)) where D q is the derivative in q-calculus and B n,q (x) is the q-Bernoulli polynomials, see [3,4]. The q-Bernoulli numbers and polynomials can be expressed as ∞ ∑ n=0 B n,q t n [n] q ! = t e q (t) − 1 , ∞ ∑ n=0 B n,q (x) t n [n] q ! = t e q (t) − 1 e q (tx), respectively.
We note that (3) becomes (2) when q → 1. Many mathematicians discovered q-differential equations by using special polynomials as a solution and studying their properties and identities, see [5,6]. For example, in [5], Hermoso, Huertas, Lastra, and Soria-Lorente studied the q-differential equation based on q-Hermit polynomials. In [6], q-differential equations and properties related to Euler polynomials and Genocchi polynomials were studied.
Based on the above paper, our purpose is to find various q-differential equations of higher order that contain q-Bernoulli polynomials as a solution of the q-differential equation of higher order. In Section 3, we find q-differential equations of higher order that have q-Bernoulli polynomials as the solution by expanding the equation in (3) and check its associated symmetric properties. In Section 4, we find the values of q-Bernoulli numbers, show the approximate roots of q-Bernoulli polynomials, and organize several conjectures for q-Bernoulli numbers and polynomials. Using q-Bernoulli polynomials, we construct a Mandelbrot set and Julia set and find the results of various figures and phenomena in Section 5.

Preliminaries
To reach the goal of this paper, we will summarize the definitions and theorems and make arrangements as follows.
Let n, q ∈ R with q = 1. The number [n] q = 1 − q n 1 − q is called q-number, see [11]. We note that lim q→1 [n] q = n. In particular, for k ∈ Z, [k] q is called q-integer. The q-Gaussian binomial coefficients are defined by where m and r are nonnegative integers, see [9,11,16]. For r = 0, the value is 1 since the numerator and the denominator are both empty products. One Definition 1. Let z be any complex numbers with |z| < 1. Two forms of q-exponential functions can be expressed as We note that lim q→1 e q (z) = e z , see [10,11]. Theorem 1. From Definition 1, we note that (i) e q (x)e q (y) = e q (x + y), if yx = qxy.
From the result of using the two concepts of q-exponential functions, new types of Bernoulli, Euler, and Genocchi polynomials appear and many mathematicians have studied their properties and identities. This topic is studied in various studies via computer, see [1][2][3][4]6,12,13,17,18]. The generating functions of q-Euler polynomials and q-Genocchi polynomials used in this paper can be confirmed in Definitions 2 and 3, see [1][2][3][4]13]. Definition 2. The generating function for the q-Euler numbers and polynomials are Let q → 1 in Definition 2. Then, we can find the Euler numbers and polynomials as Definition 3. The generating function for the q-Genocchi numbers and polynomials are Setting q → 1 in Definition 3, we can find the Genocchi numbers and polynomials as ∞ ∑ n=0 G n t n n! = 2t e t + 1 , Definition 4. The q-derivative of a function f with respect to x is defined by and D q f (0) = f (0).
We can prove that f is differentiable at zero and it is clear that D q x n = [n] q x n−1 , see [7,10,14,15]. From Definition 4, we can see some formulae for q-derivative. Theorem 2. From Definition 4, we note that , (iii) for any constants a and b, D q (a f (x) + bg(x)) = aD q f (x) + bD q g(x).

Definition 5.
We define the iterated maps of functions as follows: denote the iterated map of function f as f : C → C such that f n : z → f ( f (· · · ( f n (z)))).

Definition 6.
The orbit of the point z ∈ C under the action of the function of f is said to be bounded if there exists M ∈ R such that | f n (z)| < M for all n ∈ N. If not, it is said to be unbounded.
Let f : X → X be a complex function, with X as a subset of C. Point z is said to be a fixed point of f if f (z) = z. Point z is said to be a periodic point of the period n of f if n is the smallest natural number such that f n (z) = z. We also say that z = ∞ is a fixed point or a 1 is a periodic point. Definition 7 ([19]). Let c ∈ C be given. Then, the Mandelbrot sequence is We know the Mandelbrot set expected for the divergent sequence in the complex plane, i.e., M :

Definition 8 ([19]).
Let c ∈ C be the fixed point. The Julia sequence is then We also know the Julia set for the complex plane, i.e., K c :

Several Symmetric Properties of q-Differential Equations of Higher Order for q-Bernoulli Polynomials
In this section, we find that the q-Bernoulli polynomials are solutions to several qdifferential equations of higher order. We also derive symmetric properties of q-differential equations of higher order which are related to q-Bernoulli polynomials. In addition, we show q-differential equations of higher order which combine q-Euler numbers and q-Genocchi numbers.
Theorem 3. The q-Bernoulli polynomial B n,q (x) is a solution of the following q-differential equation of higher order.
q,x B n,q (x) where B n,q is the q-Bernoulli numbers.
Proof. Consider the q-derivative after substituting qx instead of x in q-Bernoulli polynomials. Then, we have Multiplying in Equation (4), we obtain By using q-Bernoulli numbers and Cauchy product in Equation (5), we find Comparing the coefficients of Equation (6), we have Here, we note a relation of B n,q (x) and D q,x B n,q (x) of (8) on the left hand side of (7), we have the following equation.
From Equation (9), we complete the required result.

Corollary 1.
Let k be a nonnegative integer. Then, the following holds: where B n is the Bernoulli numbers and B n (x) is the Bernoulli polynomials.
where B n is the Bernoulli numbers and B n (x) is the Bernoulli polynomials.
Theorem 4. The q-Bernoulli polynomial B n,q (x) satisfies the following q-differential equation of higher order.
q,x B n,q (q −1 x) where B n,q is the q-Bernoulli numbers.
Proof. In order to find another differential equation of higher order, we consider substituting qt for t this time. From q-Bernoulli polynomials, we have To make the calculations easier, we multiply t in Equation (10) as From Equation (11), we find Comparing both sides of t n [n] q ! in Equation (12), we have From Equation (13), we have We note a relation of q-Bernoulli polynomials and D Using the relation of Equation (15) in Equation (14), we can find where the required result is completed instantly.
From now on, we find some q-differential equation of higher order by combining q-Euler and q-Genocchi numbers.
Theorem 5. The q-Bernoulli polynomial B n,q (x) satisfies the following q-differential equation of higher order which combines q-Euler numbers and polynomials.
q,x B n,q (x) + (E 1,q + E 1,q (1))D (1) q,x B n,q (x) where E n,q is q-Euler numbers and E n,q (x) is the q-Euler polynomials.
Proof. We can consider the following equation from the generating function of q-Bernoulli polynomials.
From (16), we have a relation of q-Bernoulli polynomials and q-Euler numbers and polynomials.

Corollary 3.
Setting q → 1 in Theorem 5, the following holds: where E n is the Euler numbers and E n (x) is the Euler polynomials.
Theorem 6. The q-Bernoulli polynomial B n,q (x) satisfies the following q-differential equation of higher order which combines q-Genocchi numbers and polynomials.
q,x B n,q (x) + (G 1,q + G 1,q (1))D (1) q,x B n,q (x) where G n,q is the q-Genocchi numbers and G n,q (x) is the q-Genocchi polynomials.
Proof. From the generating function of q-Bernoulli polynomials, we have a relation between q-Bernoulli polynomials and q-Genocchi numbers and polynomials as By a similar method from Theorem 5, we can find the following equation.
Hence, we complete the proof of Theorem 6.

Corollary 4.
Putting q → 1 in Theorem 6, it holds that where G n is the Genocchi numbers and G n (x) is the Genocchi polynomials.
Proof. To find a symmetric property of q-differential equations of higher order for q-Bernoulli polynomials, suppose a form of A as .

From form A, we can have
and (20) Comparing Equations (19) and (20), we find Applying Equation (15) to Equation (21), we obtain and complete the proof of Theorem 7.

Corollary 6.
Consider q → 1 in Theorem 7. Then, the following holds: where B n is the Bernoulli numbers and B n (x) is the Bernoulli polynomials.

The Observation of q-Bernoulli Numbers and Scattering Zeros of the q-Bernoulli Polynomials
In this section, we try to find approximate values of q-Bernoulli numbers and approximate roots of q-Bernoulli polynomials which appear with changes in q. We use MATHEMATICA to pile up the structure and make some conjectures based on this.
Based on the generating function of q-Bernoulli numbers, several q-Bernoulli numbers B n,q are found as follows: , , , From q-Bernoulli numbers, Table 1 shows the approximate values of B n,q which appear with changes in q. In Table 1, we can see that if q decreases, there are approximate values of q-Bernoulli numbers near the absolute value of zero. From Table 1, we can see the location of B n,q shown by varying q and n, as shown in Figure 1. In Figure 1, nonnegative integers of the x-axis represent the value of n. For example, it means that 0 is B 0,q , 1 is B 1,q , · · · , and 15 is B 15,q . The blue dots, the yellow squares, the green rhombuses, and the red triangles in Figure 1 are the approximate values of q-Bernoulli numbers when q = 0.99, 0.6, 0.4, 0.1, respectively. The lines represent variations of the approximate values for q-Bernoulli numbers. Here, from Table 1 and Figure 1, we can observe the following.

Conjecture 1.
If the value of n increases, then the value of q-Bernoulli numbers approaches zero when q = 0.1.
Next, several q-Bernoulli polynomials B n,q (x) are shown in the following: Let us fix 0 ≤ n ≤ 50. Then, we can find the structure of the approximate roots of q-Bernoulli polynomials shown by varying q as in Figure 2 From Figure 2, it can be inferred that the structures of approximate roots for q-Bernoulli polynomials change from elliptical to circular form as q becomes smaller and n increases. In addition, it can be seen that the numbers of real roots among the approximate roots decrease as q becomes smaller in Figure 2. Table 2 shows the approximate real zeros of B n,0.2 and B n,0.5 with changes in n. In Table 2, we can see that if n ≥ 8, one of the approximate real root values is 1 in B n,0.2 . Additionally, if n ≥ 10, one of the approximate real root values is 1 in B n,0.5 . Therefore, we can organize the following. To find the exact results of the experiment, Figure 3 shows the results of the changed value of q. Let us fix n = 50. The left figure (a) in Figure 3 shows the location for approximate roots of q-Bernoulli polynomials under the condition of q = 0.1, the condition of the middle figure (b) is q = 0.5, and the condition of the right figure (c) is q = 0.9. Panels (a), (b), and (c) in Figure 3 contain the approximate real roots and the approximate imaginary roots for q-Bernoulli polynomials. Here, we expect the approximate root location for q-Bernoulli polynomials to approach any approximate diagram. Hence, we remove the approximate real roots of B 50,q of a changed value of q because these polynomials have some approximate real roots.
Let us consider n = 50 for q-Bernoulli polynomials. In (d), (e), and (f) of Figure 3, the red dots are the approximate imaginary roots for B 50,q , the blue dots are the circle's centers, and the blue lines represent the circles closest to the approximate imaginary roots. The lower left figure (d) in Figure 3 is when q = 0.1, the lower middle figure (e) is when q = 0.5, and the lower right figure (f) is when q = 0.9. Figure 3 shows the location of the approximate roots of q-Bernoulli polynomials approximated in circular form. Table 3 Table 3, we suggest a conjecture as follows.

Conjecture 3.
Let us fix q = 0.1. If n ≥ 50, then the location of approximate roots for q-Bernoulli polynomials can be found on the circle.

Visualization of the Mandelbrot Set and the Julia Set for q-Bernoulli Polynomials
In this section, we introduce the Mandelbrot set and the Julia set for q-Bernoulli polynomials and show their special properties according to the change in q. Figures of this section are investigated using the C program.
Let us consider B 4,0.5 (x) = 0.00371224 + 0.0416667x + 0.208333x 2 − 1.25x 3 + x 4 . Then, the Mandelbrot set is B 4,0.5 (x) + c. By using the escape time algorithm, the Mandelbrot set iterated 64 times can be seen in the left picture in Figure 4. As shown in (a) in Figure 4, the range of convergence is the absolute value of 2. Panel (a) in Figure 4 shows that the top and bottom are symmetrical. In Figure 4, the image center of the figure in (b) is (0.05687500, 0.56000000). Panel (b) is an enlarged part of (a) because the boundary of the Mandelbrot set shows more intricate detail as one looks closer or magnifies the image. The line in (c) in Figure 4 represents periodic critical orbit with period 4. The point (0.97500000, −0.32000000) in the Mandelbrot set is a periodic point with period 4 if its critical orbit is periodic with period 4.    Figure 5 shows the information about the accumulated points. As the point colors go from red to blue, we can see that there are more accumulated points. We show some fixed points which are periodic points with periods equal to one in (b) in Figure 5. In this figure, the green dots represent fixed points, and we can see that it consists of four dots. Additionally, it can be seen that there is a portion in which the fixed points are clustered, since the sizes of the points are different. In (c) in Figure 5, a point (0.97500000, −0.32000000) is connected to one of the attracting fixed points of the function. Figures 6 and 7 show the visualization of the Julia set's various shapes for B 4,0.5 (x) by setting the number of iterations to 128 or 64. Additionally, we consider the convergence range as 2. Panel (a) of Figure 6 shows the range of the real axis from −1.2050 to 1.7950 and the range of the imaginary axis from −1.43 to 1.57. The polynomial here is B 4,0.5 (x), the offset is shown if −0.535 − 0.27i, and the center of the image is 0.295 + 0.07i. To show this figure, we use HSV image, and light blue means low iterations and yellow means the 128th iteration. If we look closely in this figure, we can see that the 128th appearing area is indicated by a dot. Panel (b) of Figure 6 shows the range of the real axis from −1.2050 to 1.7950 and the range of the imaginary axis from −1.430 to 1.570, the offset is shown in 1.02 − 0.18i, and the convergence range of figure is 2. This figure shows a continent which is one of the representative shapes of the Julia set. The blue color shows low iterations and the green color shows the 128th iteration. In (c) of Figure 6, the range of the real axis and the range of the imaginary axis are the same condition as (b), the offset is shown in 0.6 − 0.71i, and the image center is 0.295 + 0.07i. Here, we use the gray image color and as we get closer to the black color, this means 128 iteration was complete.
Panel (a) of Figure 7 shows the range of the real axis from −0.455 to 1.045 and the range of the imaginary axis from −0.230 to 0.82 when the offset is 0.965 − 0.665i. Here, the image center of (a) in Figure 7 is 0.295 + 0.070i. We use the gray image color and as we get closer to the black color, this means the 128th iteration was complete. To find the property of self-similarity in this figure, we change the scale from 0.0005 to 0.025.
Panel (b) in Figure 7 shows the range of the real axis from −1.2050 to 1.7950 and the range of the imaginary axis from −1.430 to 1.570, the offset is shown in 0.65 + 0.835i, the maximum iteration is 64, and the convergence range of figure is 2. The red color shows low iterations and the blue color shows the 64th iteration.
In panel (b) of Figure 7, it can be seen that various periodic points appear. Panel (c) of Figure 7 shows the range of the real axis from −0.4550 to 1.0450 and the range of the imaginary axis from −0.680 to 0.820, the offset is shown if −0.30 + 0.395i, and the center of the image is 0.295 + 0.070i. To show this figure, we used HSV image. In addition, the red color means few iterations were carried out and the blue color means 128 iterations were carried out. In (c), one can find that the 128th appearing area is indicated by blue dots.

Conclusions
We find several q-differential equations of higher order which are related to q-Bernoulli polynomials and obtain a symmetric property of q-differential equations of higher order in Section 3. Since we show q-differential equations of higher order based on q-Bernoulli polynomials, we find some values of q-Bernoulli numbers by varying q and show approximate roots for q-Bernoulli polynomials. We also guess properties of q-Bernoulli polynomials in Section 4. By using the property of polynomials in which the q-Bernoulli polynomials contain various orders, the relevant Mandelbrot set and Julia set are shown in Section 5. We think it is necessary to establish a theoretical content for the conjectures and results obtained using computers.