In late December 2019, the novel coronavirus COVID-19 emerged in Wuhan, the captial of Hubei province and one of largest cities in the central part of China. One distinct feature of the COVID-19 virus is that it can spread from person to person as confirmed in [1] . There were more than 74,000 confirmed infections and more than two thousands people died until February 20, 2020 [2] . Moreover, the situations in Japan and South Korea are becoming serious as well. Hence, predicting the further spread of the COVID-19 epidemic attracts a major public attention.
Considering the epidemic's feature of spreading during the latent period, Chen [3, 4] et al. applied the time delay process to describe the typical featrue and proposed a novel dynamical system to predicte the outbreak of COVID-19. Based on the daynamical system, Chen et al. [5] also proposed a time delay dynamic system with external source to describe the trend of local outbreak for the COVID-19 and provided numerical simulations. Mathematically speaking, in their models, the time delay term involving an integral with smooth kernel could be viewed as a special kind of fractional integral in some sense.
The fractional derivatives [6] provide useful tools for the description of memory and hereditary properties of different materials. And the differential equations with fractional derivatives stand out in studying fractal geometry and fractal dynamics [7, 8] . As the fractional derivative can describe the anomalous diffusion, in this paper, we propose a novel time delay dynamic system with fractional order, which is based on the model in [3] . In the newly proposed system, the Riemann-Liouville derivative is added which can describe the confirmed and cured people's growth process.
The rest of paper is organized as follows. Section 2 presents notaions, assumptions and the corresponding fractional time delay dynamic system. The approach for estimating the future number of diagnosed people with this system is provided in section 3. Based on the offical data given by CCDC everyday, serveral numerical examples are exhibted in section 4 to verify the rationality of the fractional model and effectiveness of the estimation scheme, and to give the resonable range of the fractional order. Finally, we present some concluding remarks in section 5.
Before proceeding to a fractional time delay dynamic system, we first introduce some definitions, notaions and assumptions.
Definition 2.1 The left side Riemann-Liouville derivative is defined as [6] 
where α is the fractional order and Γ(·) denotes the standard Gamma function.
Note: In this paper, if we take fractional order α = 1, the fractional derivative reduces to the classical first order derivative.
The following semigroup property could be utilized to design a numerical algorithm later. Typically speaking for general cases, semigroup property is unvalid for Riemann-Liouville fractional operators. However, if the fractional orders are restricted in some special range, the property is valid.
Property 2.1 [9] If f (t) ∈ C 1 [0, T ], α i ∈ (0, 1) (i = 1, 2), and α 1 + α 2 ∈ (0, 1], then 0 D α1
Notations:
• I(t): cumulative infected people at time t;
• J(t): cumulative confirmed people at time t;
• G(t): currently isolated people who are infected but still in latent period at time t;
• R(t): cumulative cured people at time t;
• β: spread rate;
• γ: morbidity;
• : the rate of isolation;
• κ: cure rate.
1. The infected people averagely experience a latent period of τ 1 days before they have obvious symptoms. It is assumed that once symptoms appearing, those infected people will seek treatment and therefore become diagnosed people.
2. Because of the local government intervention, some of the infected people while still in the latent period would be isolated before they have significant symptoms. Assuming that they would be diagnosed in the next τ 1 days, which means that the average exposed period of these people are τ 1 − τ 1 ;
3. We assume that a person once isolated or in treatment, the individual would no longer transmit the coronavirus to others;
4. The confirmed people would take τ 2 days in average to be cured with rate κ;
5. We assume that growing process of J(t) and R(t) is a sub-diffusion process which means they depend on the corresponding history data from the start time to the moment.
Based on the above notations and assumptions, the dynamics of these popultaions I(t), J(t), G(t), R(t) are described by the following fractional time delay dynamic (FTDD) system:
For easiness of understanding the system (2.2)-(2.5), we present some explanations.
1. The cumulative confirmed people J(t) all come from the infected people in the latent period τ 1 ;
2. The instant change in G(t) has two factors. One is that some of infected people are isolated by the local government. Another factor is some isolated people are diagnosed after τ 1 days;
Here we take the normal distribution
with c i and c i+3 be constants.
In this section, a reconstruction algorithm of two parameters (spread rate and rate of isolation) of the dynamic system (2.2)-(2.5) is introduced based on the optimization method, utilizing the official data.
According to Property 2.1, we take 1-α derivative on both side of the equation (2.3) and 1-λ derivative on both side of the equation (2.5). Then we can transform the system (2.2)-(2.5) into
For given parameters {β, γ, , κ, τ 1 , τ 1 , τ 2 , α, λ} and the initial conditons
, R(t 0 )} of the above dynamic system (3.6), we intend to apply the Matlab innner-embedded program dde23 to solve the dynamic system numerically. Furthermore, the cumulative diagnosed people J(T ) and cumulative cured people R(T ) at any given time T are easily to be obtained.
We can also adpot implicit Euler scheme to solov the dynamic system (2.2)-(2.5) directly. Let t 0 < t 1 < · · · < t N = T , where N is the total number of intervals. This is a uniform time step τ = T /N . The shifted Grünwald formula [10] for approximationg the Riemann-Liouville derivative at time t = t n , (1 ≤ k ≤ n) is given by
Let µ = τ α+1 γβ. The following linear equations hold:
(1 − τ β)I n + τ βJ n + τ βG n = I n−1 , µh n 1 I n + (w 0 + µh n 1 )J n + µh n
where I n , J n , G n , R n are the numerical solutions of I(t n ), J(t n ), G(t n ), R(t n ), and h n i = h i (t n ). In addition, on the initial day t 0 , it is assumed that 5 people were infected the COVID-19 from unknown sources, i.e., I(t 0 ) = 5. The confirmed, isolted and recoverer people is 0, i.e., J(t 0 ) = G(t 0 ) = R(t 0 ) = 0. Moreover, it is assumed that no isolation exists before T = t + 15. According to the public data, we suppose a relatively high morbidity γ = 0.99 and cure rate as κ = 0.97. The average latent period τ 1 and treatment τ 2 are also regarded as known, i.e., τ 1 = 8 and τ 2 = 12. The average period between getting isolated and diagnosed τ 1 = 5. Therefore, the parameters that need to be estimated are The identification of parameter θ * comes to the following optimization problem:
where J Obs is the daily official data reported by the Chinese Center for Disease Control and Prevention.
The pseudo code for solving this optimization problem is given in Algorithm 1.
Input: Initial conditons {I(t 0 ), G(t 0 ), J(t 0 ), R(t 0 )}, parameters {γ, κ, τ 1 , τ 1 , τ 2 , α, λ}, initial guess θ 0 = [β 0 , 0 ], observation data J Obs , and maximum number iteration M ; Output: Predicted value I(t), J(t), G(t), R(t) and optimization parameter θ * ; 1: for 1 ≤ i ≤ M do 2:
Compute J(θ i−1 ; t) by solving (3.6) with Matlab inner-embedded program dde23 or (3.7) with Euler scheme;
Caculate θ i by solving equation (3.8) with Levenberg-Marquad (LM) method; 4: end for 5: Compute the predicted value I(t), J(t), G(t), R(t) by repeating step 2 with θ * obtained by step 3.
In this section, we shall present some numerical experiments to verify the rationality and effectiveness of our proposed model. The data including the cumulative diagnosed people J(t) and the cumulative cured people R(t) is available from January 31 untill February 14, 2020. Moreover, the data is acquired from Chinese Center for Disease Control and Prevention.
Example 1 It is significantly important to study the impact of different parameter ∆ on the dynamic model (3.6) . Therefore, we shall provide different ∆ to solve the optimization problem and analyze sentativity of the estimated parameters. The following three cases are considered. Case 1. When α > λ, Figure 1 shows the prediction results with different λ and fixed α. The reconstructed parameters are exhibited in Table 1 .
Case 2. If α = λ, the trends of cumulative confirmed and cured people are depicted in Figure 2 , and reconstructed parameters in Table 2 .
Case 3. Setting α < λ, we get the tendency of the J(t) and R(t) in Figure 3 and recovered parameters in Table 3 . RReconstr.
In treatment RReconstr.
In treatment (c) λ = 0.8 Figure 1 : Prediction with different λ when α = 1. RReconstr.
In treatment RReconstr.
In treatment (c) α = 0.96, λ = 0.96 RReconstr.
In treatment (c) α = 0.98 Figure 1-3 and Table 1-3 show α is quite sensitive as the cumulative confirmed people are sharply increasing with slight decrease of α. The estimated number of confirmed people still increase and cannot be in a stable state on early March when α < 0.965. It is worthy to mention that there is no change on the value of β and with the reduction of λ shown in Table 1 . However, the estimated number of cured people is reducing to the value that even less the people who are in treatment when λ < 0.85. Figure 2 (a) describes the predicted trend with integer order dynamic model. Comparing this results with other fractional order dynamic model, e.g., Figure 1 (a), Figure 2 (b) and Figure 3 (b), we are easy to find that the estimated cumulative confirmed people J Reconstr. with fractional order are more than that with interger order. On the other hand, the estimated cuumulative cured people R Reconstr. with fractional order are less than that with interger order.
In fact, many infected people are not easy to be diagnosed by Nucleic acid detection, which leads to the number of confirmed people less than real number of people who should be diagnosed. People who are cured may still relapse. The predicted results of fractional order dynamic model can exactly explain this phenomenon. In addition, we infer that 0.965 < α ≤ 1 and 0.84 < λ ≤ 1 are more appropriate.
Example 2 Let α = 0.97, λ = 0.87. We provide some numerical results and reconstructed parameters with official data of two areas. One area is the Hubei Province. Another area is the outside of Hubei Provice while still in China. Besides, in order to verify the accuracy of the chosen fractional orders, we use the latest official data (from February 15 to February 19) to compare the results of prediction. The predicted values of confirmed and cured people are shown in Figure 4 . In this figure, J Reconstr. and R Reconstr. are the estimated values of cumulative diagnosed and cured people for prediction. J obs and R obs correspond to the data of cumulative confirmed and cured people for reconstructions. J test and R test are displayed for comparison with the prediction. 
In this paper, we propose a novel time delay dynamic system with fractional derivative, where the growth of cumulative confirmed and cured people at time t are both considered as a sub-diffusion process. Numerical simulations are carried out to verify the rationality and effectivness of the dynamic system.In addition, two proper interval are given for the fractional orders α and λ, respectively. Numerical examples show that the proposed fractional time delay dynamic system can generate prediction which agrees the public data well.