Superconvergence for optimal control problems governed by semilinear parabolic equations

1.   Introduction

Grey system theory is a method for forecasting and analyzing small sample sequences. It was first put forward by Deng ^[1,2] and has been widely used in various fields. In order to enhance the adaptability of the prediction model, Deng proposed a single-variable grey prediction model GM(2,1), which stands for one variable and second-order differential equation. However, due to defects in the structure, the prediction effect is not ideal in some cases ^[21]. In order to improve the prediction accuracy of the model, scholars put forward some improved models based on the classical GM(2,1) model. Shen et al. ^[3] proposed a method to determine the system parameters using the least square method to reduce error. Zeng et al. ^[4] proposed a new parameter estimation to GM(2,1) model based on accumulating method which reduces the matrix calculation quantity. A new healthier model was obtained and the serious morbidity problem with a multiple transformation on the initial data was resolved. In ^[5], two improved models are proposed based on modeling and predicting on non-linear rolling motion by GM(2,1), which are grey cycle extension and grey neural network respectively based on error compensation. Li et al. ^[6] introduced the GM(2,1) model and its improved algorithm and prediction steps by Matlab. By analyzing the shortcomings of GM(2,1) model, an optimized SOGM(2,1) model was constructed in ^[7] and applied to the prediction of high-speed roadbed settlement with higher prediction accuracy. In ^[8], Zhou et al. focused on the two shortcomings of GM(2,1) model and proposed appropriate measures for improvement. Moreover Zhou improved the method, which was used in solving albino equations of the classical GM(2,1) model, to solve the same equations by Laplace transform and improve the prediction accuracy. In ^[9], based on the GM(2.1) model, the boundary conditions were improved, and the effectiveness of the improved method was verified with grey data.

In ^[10], Wu et al. first discussed using perturbation theory and the least square method to address a violation of the new information priority principle in grey system theory. Ihis introduced a new grey system model with fractional order accumulation that reflects new information priority better with smaller accumulation orders. In recent years, various grey prediction models based on fractional order accumulation have been proposed in ^{[11,12,13,14,15]}. In ^[12], Wu improved the GM(2,1) using fractional order AGO to achieve an accurate prediction. Based on the GM(2,1) model with fractional order accumulation, Zeng ^[15] gave two different ways to determine the time response parameters of the fractional order GM(2,1) model, which can be flexibly selected according to the specific results in practical application.

The B-spline method has been widely applied to derive numerical solutions of partial differential equations. In ^[16], a cubic monotonicity-preserving interpolation spline was proposed to reconstruct the background value. In ^[17], Chen and Zhu developed a new GM(1,1) model based on $C^1$ monotonicity-preserving rational linear interpolation spline to provide a more reasonable method to calculate the background value. In ^[18], in order to improve the GM (1,1) model, a formula based on a $C^1$ monotonicitypreserving piecewise rational quadratic interpolation spline for calculating background value was proposed.

In the classical $GM^2(2,1)$ model, the corresponding estimated function $x^{(2)}(t)$ has the property of convexity and $C^\infty$ continuity. The classical GM$^2$(2,1) model use the piecewise linear polynomial interpolation $L_k(t) = \left({k + 1 - t} \right){x^{(2)}}(k) + \left({t - k} \right){x^{(2)}}(k + 1)$, $k = 1, 2, \ldots, n-1$ as an approximate solution of $x^{(2)}(t)$ so as to estimate background values ${z^{(2)}}(k + 1)$ and the derivative values $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k+1}}$ and $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k}}$. It is well known that the piecewise linear polynomial interpolation has convexity-preserving property, but it just has $C^0$ continuity and only has $O(h^2)$ order of convergence. In estimating $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = 1}}$, the classical GM$^2$(2,1) model simply sets ${x^{(2)}}(0) = 0$, whereas this will cause error in estimating $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = 1}}$.

The main contribution of this work is to use a $C^1$ convexity-preserving rational quadratic interpolation spline given in ^[19] as approximate solution of $x^{(2)}(t)$ so as to estimate background values ${z^{(2)}}(k + 1)$ and the derivative values $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k+1}}$ and $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k}}$, $k = 1, 2, \ldots, n-1$ and thereby a new GM$^2$(2,1) model is established. The rest of this paper is structured as follows. Section 2 recall the construction of the $C^1$ convexity-preserving interpolation spline given in ^[19]. In section 3, on the basis of $C^1$ convexity-preserving rational quadratic interpolation spline, a new GM$^2$(2,1) model is constructed in detail. Several numerical examples are given in section 4. The conclusion is given in section 5.

2.   Classical $GM^2(2,1)$ model

In this subsection, we recall the classical $GM^2(2,1)$ model. Let an original non-negative and uniformly-spaced sequence be

The r-AGO sequence $X^{(r)}$ is given as follows

where

From ^[10], by using the fractional order accumulation operator, we can also rewrite ${x^{(r)}}(k)$ as ${x^{(r)}}(k) = \sum\limits_{i = 1}^k {\left({\begin{array}{*{20}{c}} {k - i + r - 1}\\ {k - i} \end{array}} \right)} {x^{(0)}}(i), k = 1, 2, \ldots, n$. Here, $\left({\begin{array}{*{20}{c}} m\\ n \end{array}} \right): = \frac{{m!(m - n)!}}{{n!}}$ and $\left({\begin{array}{*{20}{c}} 0\\ n \end{array}} \right): = 1, \left({\begin{array}{*{20}{c}} {n + 1}\\ n \end{array}} \right): = 0$. In the following discussion, we mainly focus on $r = 2$.

From Eq (2.1), we can obtain 1-AGO sequence $X^{(1)}$ and 2-AGO sequence $X^{(2)}$ as follows

where

From Eq (2.2), we can see that the 2-AGO sequence $X^{(2)}$ has convexity property, that is ${x^{(1)}}(k + 2) = {x^{(2)}}(k + 2) - {x^{(2)}}(k + 1) \ge {x^{(2)}}(k + 1) - {x^{(2)}}(k) = {x^{(1)}}(k + 1)$. Suppose that $x^{(2)}(t)$ meets the following second order grey forecasting differential equation

where $\alpha = (a_1, a_2, b)^{\mathsf{T}}$ is the parameter in the model to be estimated.

The corresponding characteristic function of Eq (2.3) is $r^2+a_1r+a_2 = 0$. Consider the case that $\Delta = a_1^2-4a_2 > 0$, then the roots of the characteristic function are $r_{1, 2} = \frac{-a_1 \pm \sqrt {\Delta}}{2}$.

Thus, the solution of Eq (2.3) with the initial condition ${\widetilde x^{(2)}}(1) = {x^{(2)}}(1), \widetilde x^{(2)}(n) = {x^{(2)}}(n)$ is as follows

where

Therefore, to obtain the prediction model of the original data sequence, we need to identify the effect of $\alpha$. For this purpose, we do the integral accumulation on both sides of Eq (2.3) in each interval $\left[{k, k + 1} \right], k = 1, 2, \ldots, n - 1$, then we can get

that is

Let background value be ${z^{(2)}}(k + 1) : = \int_k^{k + 1} {{x^{(2)}}(t)} dt$. In order to calculate the background value ${z^{(2)}}(k + 1)$, we need to integrate ${x^{(2)}}(t)$, which requires the value of $\alpha$ to be given in advance. However, the value of $\alpha$ needs to be determined from the Eq (2.4). Consequently, to estimate the value of $\alpha$, we must use some methods to estimate the background values ${z^{(2)}}(k + 1)$, $k = 1, 2, \cdots, n-1$. In the classical GM$^2$(2,1) model, for $t\in[k, k+1], k = 1, 2, \ldots, n-1$, we use the piecewise linear polynomial interpolation $L_k(t) : = \left({k + 1 - t} \right){x^{(2)}}(k) + \left({t - k} \right){x^{(2)}}(k + 1)$ as an approximate solution of $x^{(2)}(t)$, that is ${x^{(2)}}(t) \approx {L_k}(t)$ for $t \in [k, k + 1], k = 1, 2, \ldots, n - 1$. Then we get the estimated background values ${z^{(2)}}(k + 1)$, $k = 1, 2, \ldots, n-1$ as follows

And for the derivative values, we approximately have

For $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = 1}}$, by setting ${x^{(2)}}(0) = 0,$ then we get $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = 1}} \approx {x^{(2)}}(1) - {x^{(2)}}(0) = {x^{(2)}}(1)$.

For each interval $[k, k+1]$, $k = 1, 2, \ldots, n-1$, by substituting the estimated background values ${z^{(2)}}(k + 1)$ given in Eq (2.5) and the approximate derivative values $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k+1}}$ and $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k}}$ into Eq (2.4), we have

Further applying the least square method, we get the estimated value of $\alpha$ by the formula as follows

where

With the estimated parameter $\alpha = (a_1, a_2, b)^{\mathsf{T}}$, we further get the solution ${\widetilde x^{(2)}}(t)$ of the second order grad forecasting differential equation given in (2.3) with the initial condition ${\widetilde x^{(2)}}(1) = {x^{(2)}}(1), \; \widetilde x^{(2)}(n) = {x^{(2)}}(n)$. Finally, we get the following classical $GM^2(2,1)$ prediction formula

3.   $C^1$ convexity-preserving interpolation spline

Given convexity data set $\{ ({x_k}, {y_k}) \in {R^2}\} _{k = 1}^n$ with ${\Delta _1} < {\Delta _2} < \ldots < {\Delta _{n - 1}}$, where $x_1 < x_2 < \ldots < x_n$ and ${\Delta _k} = [{y_{k + 1}}-{y_k}]/[{x_{k + 1}}-{x_k}], k = 1, 2, \ldots, n - 1$. For $x\in[x_k, x_{k+1}]$, $k = 1, 2, \ldots, n-1$, Let ${h_k} = {x_{k + 1}} - {x_k}, t = (x - {x_k})/{h_k}, {A_k} = {\Delta _k} - {d_k}, {B_k} = {d_{k + 1}} - {\Delta _k},$ then the $C^1$ convexity-preserving rational quadratic interpolation spline $S(x)$ developed in ^[19] is as follows

where

From Eq (3.1), it is easy to check that $S(x_k^ +) = {y_k}, S(x_{k + 1}^ -) = {y_{k + 1}}, S'(x_k^ +) = {d_k}, S'(x_{k + 1}^ -) = {d_{k + 1}}$, $k = 1, 2, \ldots, n-1$ which implies that $S(x)\in C^1[x_1, x_n]$. For ${\Delta _1} < {\Delta _2} < \ldots < {\Delta _{n - 1}}$, from ^[19], we see that $S''(x) > 0$, which means that $S(x)$ has convexity-preserving properties. Furthermore, let $f(x)\in C^3[x_1, x_n]$ be a given convex function with which $S_(x)$ is compared and $y_i = f(x_i), i = 1, 2, \ldots, n,$ from ^[19], we have $f(x) - S(x) = O({h^3})$, which indicates that interpolant has $O({h^3})$ convergence. Since the $C^1$ convexity-preserving rational quadratic interpolation spline $S(x)$ developed in ^[19] has many merits, in the following we shall use it as approximate solution of $x^{(2)}(t)$ so as to estimate background values ${z^{(2)}}(k + 1)$ and the derivative values $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k+1}}$ and $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k}}$ and a method for estimating $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = 1}}$ will be given.

4.   Establish new GM$^2$(2,1) model

For the original non-negative sequence ${X^{(0)}} = \left\{ {{x^{(0)}}(1), {x^{(0)}}(2), \ldots, {x^{(0)}}(n)} \right\}$, we first calculate its 1-AGO sequence ${X^{(1)}} = \left\{ {{x^{(1)}}(1), {x^{(1)}}(2), \ldots, {x^{(1)}}(n)} \right\}$, 2-AGO sequence ${X^{(2)}} = \left\{ {{x^{(2)}}(1), {x^{(2)}}(2), \ldots, {x^{(2)}}(n)} \right\}$. From the above discussion, the 2-AGO sequence $X^{(2)}$ is a convex data set. The corresponding estimated function $x^{(2)}(t)$ has the property of convexity and $C^\infty$ continuity. Thus, we consider to reconstruct $x^{(2)}(t)$ with the $C^1$ convexity-preserving rational quadratic interpolation spline developed in ^[19]. For the convex 2-AGO sequence data set $\left\{ {(x_k, {x^{(2)}}(k)) \in {R^2}} \right\}_{k = 1}^n$ with $x_k = k$, for each of the subintervals $\left[{k, k + 1} \right]$, by using the $C^1$ convexity-preserving rational quadratic interpolation spline given in (3.1) as approximate solution of $x^{(2)}(t)$, $t\in\left[{k, k + 1} \right]$, $k = 1, 2, \ldots, n-1$, we have

where $s = k+1-t\in[0, 1]$. From (3.2), we have

and

where $k = 2, 3, \ldots, n-1$.

Then we estimate the background values ${z^{(2)}}\left({k + 1} \right) = \int_k^{k + 1} {{x^{(2)}}\left(t \right)} dt$, $k = 1, 2, \ldots, n-1$, by the following method

Moreover, we further estimate the derivative values $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k+1}}$ and $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k}}$, $k = 1, 2, \ldots, n-1$ by the following method

where the derivative values $d_k$, $k = 1, 2, \ldots, n$ are given by Eq (4.2).

Then, by substituting the estimated background values ${z^{(2)}}\left({k + 1} \right) = \int_k^{k + 1} {{x^{(2)}}\left(t \right)} dt$ and the estimated the derivative values $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k+1}}$ and $\frac{{d{x^{(2)}}(t)}}{{dt}}{|_{t = k}}$, $k = 1, 2, \ldots, n-1$ into the grey differential Eq (2.4) and applying the following least square method to solve Eq (2.4), we get

where

with $d_k$, $k = 1, 2, \ldots, n$ are given by Eq (4.2) and ${z^{(2)}}\left({k + 1} \right) = \int_k^{k + 1} {{x^{(2)}}\left(t \right)} dt$ $k = 1, 2, \ldots, n-1$ are given by Eq (4.3).

With the estimated parameter $\alpha = (a_1, a_2, b)^{\mathsf{T}}$ given by (4.5), we further solve the second order grad forecasting differential equation given in (2.3) with the initial condition ${\widetilde x^{(2)}}(1) = {x^{(2)}}(1), \; \widetilde x^{(2)}(n) = {x^{(2)}}(n)$, whose solution is denoted as ${\widetilde x^{(2)}}(t)$. Finally, we get the following new $GM^2(2,1)$ grey prediction formula

Summarize the above discussion, we give the following algorithm for establish our new grey prediction equation.

5.   Numerical examples

We shall give several examples to show that the new GM$^2$(2,1) model based on $C^1$ convexity-preserving rational quadratic interpolation spline has better predict accuracy than the classical GM(2,1) model and GM$^2$(2,1) model, Zeng's GM$^2$(2,1) model given in ^[15] and Xu and Dang's GM(2,1) model given in ^[7]. In the following examples, the relative error is computed by

Example 1. We take the total water consumption of Beijing from 2003 to 2019 in ^[20] as an example. We compare our new GM$^2$(2,1) model with the classical GM(2,1) model, the classical GM$^2$(2,1) model and Zeng's GM$^2$(2,1) model given in ^[15]. We use the data from 2003 to 2015 as the training data and the data from 2016 to 2019 as the testing data. Table 1 gives the numerical results. The results show that the new GM$^2$(2,1) has improved prediction accuracy compared to the classical GM(2,1) model, the classical GM$^2$(2,1) model and Zeng's GM$^2$(2,1) model.

Example 2. We take the precipitation data of Xichang City from 1986 to 1990 in ^[20] as the second example. We use the data from 1986 to 1989 as the training data and the data in 1990 as the testing data. From the numerical results presented in Table 2, it can be seen that compared with the classical GM(2,1) model, the classical GM2 (2,1) model and Zeng's GM2 (2,1) model given in ^[15], our new GM2 (2.1) has the best performance.

Example 3. In this example, we consider the historical syphilis incidence data of China from 2000 to 2008 given in ^[22].We use the data from 2000 to 2005 as the training data and the data from 2006 to 2008 as the testing data. We compare our new GM$^2$(2,1) model with the classical GM(2,1) model, the classical GM$^2$(2,1) model and Zeng's GM$^2$(2,1) model given in ^[15]. Table 3 gives the numerical results. From the results, compared with the classical GM(2,1) model, the classical GM$^2$(2,1) model and Zeng's GM$^2$(2,1) model, the new GM$^2$(2,1) has better performance.

Example 4. We take the data of inbound tourists (in 10,000 people) in Beijing from 2012 to 2019 as an example. The data comes from the Beijing Statistical Yearbook in Beijing Municipal Bureau of Statistics and Survey Office of the National Bureau of Statistics in Beijing. The unit is 10,000 people. We use the data from 2012 to 2015 as the training data and the data from 2016 to 2019 as the testing data. We compare our new GM$^2$(2,1) model with the classical GM(2,1) model, the classical GM$^2$(2,1) model and Zeng's GM$^2$(2,1) model given in ^[15]. Table 4 gives the numerical results. From the results, compared with the classical GM(2,1) model, the classical GM$^2$(2,1) model and Zeng's GM$^2$(2,1) model, the new GM$^2$(2,1) has the best prediction accuracy.

Example 5. In this example, we consider the data of tourists staying overnight in Guangzhou (in 10,000 people) from 2010 to 2019, which is from the Guangzhou Statistical Yearbook of Guangzhou Municipal Bureau of Statistics.The data from 2010 to 2015 are taken as the training data and the data from 2016 to 2019 are used as the testing data. Best prediction accuracy using our new GM2(2,1) model can be founded from the numerical results presented in Table 5, compared with the classical GM(2,1) model, the classical GM2 (2,1) model and Xu and Dang's GM(2,1) model given in ^[7].

Example 6. We take the data on China's foreign exchange reserves (in trillion dollars) from 2010 to 2021. The data is from the China Statistical Yearbook of National Bureau of Statistics of China. We use the data from 2010 to 2018 as the training data and the data from 2019 to 2021 as the testing data. We compare our new GM$^2$(2,1) model with the classical GM(2,1) model, the classical GM$^2$(2,1) model and Xu and Dang's GM(2,1) model given in ^[7]. Table 6 gives the numerical results. From the results, compared with the classical GM(2,1) model, the classical GM$^2$(2,1) model and Xu and Dang's GM(2,1) model, the new GM$^2$(2,1) has the best prediction accuracy.

Example 7. We take the data of inbound tourists (in 10,000 people) in Zaozhuang from 2005 to 2019 as an example. The data comes from the Ministry of Culture And Tourism of the People's Republic of China. The unit is 10,000 people. We use the data from 2005 to 2017 as the training data and the data from 2018 to 2019 as the testing data. We compare the new GM$^2$(2,1) model with the classical GM(2,1) model, the classical GM$^2$(2,1) model and the OFOPGM model given in ^[23]. Table 7 gives the numerical results. From the results, the new GM2(2,1) model has better prediction accuracy.

According to the results of the above numerical examples1–6, the prediction accuracy of our new GM$^2$(2,1) model is improved for all the numerical examples compared to the classical GM(2,1) model, the classical GM$^2$(2,1) model, Zeng's GM$^2$(2,1) model given in ^[15], Xu and Dang's GM(2,1) model given in ^[7] and the OFOPGM model given in ^[23]. There are different degrees of improvement for different data features. Based on the above data characteristics, we make the conclusion that the more convexity the original datas presents the higher the prediction accuracy of the new GM$^2$(2,1) model, the better the actual prediction effect.

6.   Conclusions

Using the $C^1$ convexity-preserving rational quadratic interpolation spline developed in ^[19] to estimate the background values and derivative values, we established a new GM$^2$(2,1) model. Numerical examples show that the new GM$^2$(2,1) model can improve the forecasting quality, especially in prediction reliability and this model performs better when the original data are presented with convexity in time series. There is a limitation in the new method: The computations concerning the background values and derivative values take more time than the methods included in classical GM$^2$(2,1) model. Future work will concentrate on applying $C^1$ convexity-preserving rational quadratic interpolation spline to establish a new GM$^r$(2,1) model with fractional $r$.

Use of AI tools declaration

The author declares he has not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research is supported by the National Natural Science Foundation of China (No. 61802129).

Conflict of interest

The author declares no conflict of interest.