About this Journal Submit a Manuscript Table of Contents
Abstract and Applied Analysis
Volume 2012 (2012), Article ID 276080, 26 pages
http://dx.doi.org/10.1155/2012/276080
Research Article

A Note on the Inverse Problem for a Fractional Parabolic Equation

Department of Mathematics, Fatih University, 34500 Buyukcekmece, Istanbul, Turkey

Received 15 May 2012; Accepted 8 July 2012

Academic Editor: Ravshan Ashurov

Copyright © 2012 Abdullah Said Erdogan and Hulya Uygun. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

For a fractional inverse problem with an unknown time-dependent source term, stability estimates are obtained by using operator theory approach. For the approximate solutions of the problem, the stable difference schemes which have first and second orders of accuracy are presented. The algorithm is tested in a one-dimensional fractional inverse problem.

1. Introduction

Inverse problems arise in many fields of science and engineering such as ion transport problems, chromatography, and heat determination problems with an unknown internal energy source. Different typed of inverse problems have been investigated, and the main results obtained in this field of research were given by many researchers (see [110]). More than three centuries the theory of fractional derivatives developed mainly as a pure theoretical field of mathematics. Fractional integrals and derivatives appear in the theory of control of dynamical systems, when the controlled system or/and the controller is described by a fractional differential equation (see [11]). Recently, many application areas such as bioengineering applications, image and signal processing are also related to fractional calculus. Methods of solutions of problems and theory of fractional calculus have been studied by many researchers [1128]. Among them finite difference method is used for solving several fractional differential equations (see [20, 22, 23, 27] and the references therein).

1.1. Statement of the Problem

Many scientists and researchers are trying to enhance mathematical models of real-life cases for investigating and understanding the behavior of them. Therefore, some phenomena have been modeled and investigated as fractional inverse problems (see [2933] and the references therein). In this paper, we consider the fractional parabolic inverse problem with the Dirichlet condition 𝜕𝑢(𝑡,𝑥)𝜕𝜕𝑡𝑎2𝑢(𝑡,𝑥)𝜕𝑥2𝐷𝑡1/2𝑢𝑢(𝑡,𝑥)+𝜎𝑢(𝑡,𝑥)=𝑝(𝑡)𝑞(𝑥)+𝑓(𝑡,𝑥),0<𝑥<𝜋,0<𝑡𝑇,𝑢(𝑡,0)=𝑢(𝑡,𝜋)=0,0𝑡𝑇,𝑢(0,𝑥)=𝜑(𝑥),0𝑥𝜋,𝑡,𝑥=𝜌(𝑡),0<𝑥<𝜋,(1.1) where 𝑢(𝑡,𝑥) and 𝑝(𝑡) are unknown functions, 𝑎(𝑥)𝑎>0, and 𝜎>0 is a sufficiently large number. Here, 𝐷𝑡1/2=𝐷1/20+ is the standard Riemann-Liouville’s derivative of order 1/2.

Theorems on the stability of problem (1.1) are analyzed by assuming that 𝑞(𝑥) is a sufficiently smooth function, 𝑞(0)=𝑞(𝜋)=0 and 𝑞(𝑥)0.

2. Main Results

In this section, stability estimates for the solution of (1.1) are investigated. For the mathematical substantiation, we introduce the Banach space 𝐶𝛼[0,𝜋],𝛼(0,1), of all continuous functions 𝜙(𝑥) defined on [0,𝜋] with 𝜙(0)=𝜙(𝜋)=0 satisfying a Hölder condition for which the following norm is finite 𝜙𝐶𝛼[0,𝜋]=𝜙𝐶[0,𝜋]+sup0<𝑥<𝑥+<𝜋||||𝜙(𝑥+)𝜙(𝑥)𝛼,(2.1) where 𝐶[0,𝜋] is the space of all continuous function 𝜙(𝑥) defined on [0,𝜋] with the norm 𝜙𝐶[0,𝜋]=max0𝑥𝜋||||.𝜙(𝑥)(2.2) With the help of a positive operator 𝐴, we introduce the fractional spaces 𝐸𝛼,0<𝛼<1, consisting of all 𝑣 in a Banach space 𝐸 for which the following norm is finite: 𝑣𝐸𝛼=𝑣𝐸+sup𝜆>0𝜆1𝛼𝐴exp{𝜆𝐴}𝑣𝐸.(2.3) Throughout the paper, positive constants will be indicated by 𝑀𝑖(𝛼,𝛽,). Here variables are used to focus on the fact that the constant depends only on 𝛼,𝛽, and the subindex 𝑖 is used to indicate a different constant.

Theorem 2.1. Let 𝜑𝐶2𝛼+2[0,𝜋],𝐹𝐶([0,𝑇],𝐶2𝛼[0,𝜋]), and 𝜌𝐶[0,𝑇]. Then for the solution of problem (1.1), the following coercive stability estimates 𝑢𝑡𝐶([0,𝑇],𝐶2𝛼[0,𝜋])+𝑢𝐶([0,𝑇],𝐶2𝛼+2[0,𝜋])𝑥𝑀𝜌,𝑞𝐶[0,𝑇]+𝑀𝑎,𝛿,𝜎,𝛼,𝑥×,𝑞,𝑇𝜑𝐶2𝛼+2[0,𝜋]+𝐹𝐶([0,𝑇],𝐶2𝛼[0,𝜋])+𝜌𝐶[0,𝑇],𝑝𝐶[0,𝑇]𝑥𝑀𝜌,𝑞𝐶[0,𝑇]+𝑀𝑎,𝛿,𝜎,𝛼,𝑥×,𝑞,𝑇𝜑𝐶2𝛼+2[0,𝜋]+𝐹𝐶([0,𝑇],𝐶2𝛼[0,𝜋])+𝜌𝐶[0,𝑇](2.4) hold.

Proof. Let us search for the solution of inverse problem (1.1) in the following form (see [8]): 𝑢(𝑡,𝑥)=𝜂(𝑡)𝑞(𝑥)+𝑤(𝑡,𝑥),(2.5) where 𝜂(𝑡)=𝑡0𝑝(𝑠)𝑑𝑠.(2.6) Using the overdetermined condition, we get 𝜌𝜂(𝑡)=(𝑡)𝑤𝑡,𝑥𝑞(𝑥),𝜌(2.7)𝑝(𝑡)=(𝑡)𝑤𝑡𝑡,𝑥𝑞(𝑥).(2.8) Using identity (2.8) and the triangle inequality, it follows that ||||=||||𝜌𝑝(𝑡)(𝑡)𝑤𝑡𝑡,𝑥𝑞(𝑥)||||𝑥𝑀||𝜌,𝑞||+||𝑤(𝑡)𝑡𝑡,𝑥||𝑥𝑀,𝑞max0𝑡𝑇||𝜌||(𝑡)+max0𝑡𝑇max0𝑥𝜋||𝑤𝑡||𝑥(𝑡,𝑥)𝑀,𝑞max0𝑡𝑇||𝜌||(𝑡)+max0𝑡𝑇𝑤𝑡(𝑡,)𝐶2𝛼[0,𝜋](2.9) for any 𝑡,𝑡[0,𝑇]. Here, 𝑤(𝑡,𝑥) is the solution of the following problem: 𝜕𝑤(𝑡,𝑥)𝜕𝜕𝑡𝑎2𝑤(𝑡,𝑥)𝜕𝑥2𝜌𝑎(𝑡)𝑤𝑡,𝑥𝑞(𝑥)𝑑2𝑞(𝑥)𝑑𝑥2𝐷𝑡1/2𝐷𝑤(𝑡,𝑥)𝑡1/2𝜌(𝑡)𝐷𝑡1/2𝑤𝑡,𝑥𝑞(𝑥)𝑞(𝑥)+𝜎𝜌(𝑡)𝑤𝑡,𝑥𝑞(𝑥)𝑞(𝑥)+𝜎𝑤(𝑡,𝑥)=𝑓(𝑡,𝑥),0<𝑥<𝜋,0<𝑡𝑇,𝑤(𝑡,0)=𝑤(𝑡,𝜋)=0,0𝑡𝑇,𝑤(0,𝑥)=𝜑(𝑥),0𝑥𝜋.(2.10) For simplicity, we assign 𝐹(𝑡,𝑥)=𝑎𝜌(𝑡)𝑞(𝑥)𝑑2𝑞(𝑥)𝑑𝑥2𝜎𝜌(𝑡)𝑞(𝑥)𝐷𝑞(𝑥)+𝑡1/2𝜌(𝑡)𝑞(𝑥)𝑞(𝑥)+𝑓(𝑡,𝑥),𝐺(𝑡,𝑥)=𝑄1𝑞,𝜌,𝑥,𝑥𝑤,𝑡𝑡,𝑥+𝑄2𝑞,𝑥,𝑥𝐷𝑡1/2𝑤𝑡,𝑥+𝐷𝑡1/2𝑤(𝑡,𝑥),(2.11) where 𝑄1𝑞,𝜌,𝑥,𝑥=1,𝑡𝑞(𝑥)𝑑𝑎2𝑞(𝑥)𝑑𝑥2,𝑄+𝜎𝜌(𝑡)2𝑞,𝑥,𝑥𝑞=(𝑥)𝑞(𝑥).(2.12) Note that functions 𝐹(𝑡,𝑥),𝑄1(𝑞,𝜌,𝑥,𝑥,𝑡) and 𝑄2(𝑞,𝑥,𝑥) only contain given functions. Then, we can rewrite problem (2.10) as 𝜕𝑤(𝑡,𝑥)𝜕𝜕𝑡𝑎2𝑤(𝑡,𝑥)𝜕𝑥2+𝜎𝑤(𝑡,𝑥)=𝐹(𝑡,𝑥)+𝐺(𝑡,𝑥),0<𝑥<𝜋,0<𝑡𝑇,𝑤(𝑡,0)=𝑤(𝑡,𝜋)=0,0𝑡𝑇,𝑤(0,𝑥)=𝜑(𝑥),0𝑥𝜋.(2.13) So, the end of proof of Theorem 2.1 is based on estimate (2.9) and the following theorem.

Theorem 2.2. For the solution of problem (2.10), the following coercive stability estimate 𝑤𝑡𝐶2𝛼[0,𝜋]𝑀𝑎,𝛿,𝜎,𝛼,𝑥×,𝑞,𝑇𝜑𝐶2𝛼+2[0,𝜋]+𝐹𝐶([0,𝑇],𝐶2𝛼[0,𝜋])+𝜌𝐶[0,𝑇](2.14) holds.

Proof. In a Banach space 𝐸=𝐶[0,𝜋], with the help of the positive operator 𝐴 defined by 𝜕𝐴𝑢=𝑎(𝑥)2𝑢(𝑡,𝑥)𝜕𝑥2+𝜎𝑢,(2.15) with 𝐷𝑢(𝐴)=(𝑥)𝑢,𝑢,𝑢[]𝐶0,𝜋,𝑢(0)=𝑢(𝜋)=0,(2.16) where 𝜎 is a positive constant, problem (2.10) can be written in the abstract form as an initial-value problem 𝑤𝑡𝑤+𝐴𝑤=𝐹(𝑡)+𝐺(𝑡),0<𝑡𝑇,(0)=𝜑.(2.17) By the Cauchy formula, the solution can be written as 𝑤(𝑡)=𝑒𝑡𝐴𝜑𝑡0𝑒(𝑡𝑠)𝐴(𝐹(𝑠)+𝐺(𝑠))𝑑𝑠.(2.18) Applying the formula 𝐷𝑡1/2𝑢(𝑡)=𝑡0𝑢(𝜁)𝑑𝜉𝜋(𝑡𝜉)1/2,(2.19) we get the following presentation of the solution of abstract problem (2.17): 𝐷1/2𝑤(𝑡)=𝑡0𝐴𝑒𝜉𝐴𝜑𝜋(𝑡𝜉)1/2𝑑𝜉𝑡0𝐹(𝜉)𝜋(𝑡𝜉)1/2𝑑𝜉𝑡0𝐺(𝜉)𝜋(𝑡𝜉)1/2𝑑𝜉+𝑡0𝜉0𝐴𝑒(𝜉𝑠)𝐴𝐹(𝑠)𝜋(𝑡𝜉)1/2+𝑑𝑠𝑑𝜉𝑡0𝜉0𝐴𝑒(𝜉𝑠)𝐴𝐺(𝑠)𝜋(𝑡𝜉)1/2𝑑𝑠𝑑𝜉.(2.20) Changing the order of integration, we obtain that 𝐷1/2𝑤(𝑡)=𝑡0𝐴𝑒𝜉𝐴𝜑𝜋(𝑡𝜉)1/2𝑑𝜉𝑡0𝐹(𝜉)𝜋(𝑡𝜉)1/2+𝑑𝜉𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐹(𝑠)𝑑𝑠𝑡0𝐺(𝜉)𝜋(𝑡𝜉)1/2+𝑑𝜉𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐺(𝑠)𝑑𝑠=5𝑘=1𝐽𝑘,(2.21) where 𝐽1(𝑡)=𝑡0𝐴𝑒𝜉𝐴𝜑𝜋(𝑡𝑝)1/2𝐽𝑑𝜉,2(𝑡)=𝑡0𝐹(𝜉)𝜋(𝑡𝜉)1/2𝐽𝑑𝜉,3(𝑡)=𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝐽𝑑𝜉𝐹(𝑠)𝑑𝑠,4(𝑡)=𝑡0𝐺(𝜉)𝜋(𝑡𝜉)1/2𝐽𝑑𝜉,5(𝑡)=𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐺(𝑠)𝑑𝑠.(2.22) Now, we estimate 𝐽𝑘(𝑡),𝑘=1,2,3,4,5 separately. It is known that [13] 𝐴𝛼𝑒𝑡𝐴𝐸𝐸𝑀,0𝛼1.(2.23) Since operators 𝐴 and exp(𝑡𝐴) commute, 𝐴𝑒𝑡𝐴𝜑𝐸𝛼𝑒𝑡𝐴𝐸𝛼𝐸𝛼𝐴𝜑𝐸𝛼𝑒𝑡𝐴𝐸𝐸𝐴𝜑𝐸𝛼.(2.24) Applying the definition of norm of the spaces 𝐸𝛼 and (2.23) and (2.24), we get 𝐽1(𝑡)𝐸𝛼=𝑡0𝐴𝑒𝜉𝐴𝜑𝜋(𝑡𝑝)1/2𝑑𝜉𝐸𝛼𝑀1𝐴𝜑𝐸𝛼(2.25) for any 𝑡,𝑡[0,𝑇]. Estimation of 𝐽2(𝑡) is as follows: 𝐽2(𝑡)𝐸𝛼=𝑡0𝐹(𝜉)𝜋(𝑡𝜉)1/2𝑑𝜉𝐸𝛼𝐹(𝑡)𝐶(𝐸𝛼)𝑡01𝜋(𝑡𝜉)1/2𝑑𝜉𝑀2𝐹𝐶(𝐸𝛼).(2.26) Let us estimate 𝐽3(𝑡): 𝐽3(𝑡)𝐸𝛼=𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐹(𝑠)𝑑𝑠𝐸𝛼𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝑑𝑠𝐸𝛼𝐸𝛼𝐹𝐶(𝐸𝛼).(2.27) It is proven that (see [28]) 𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐸𝐸𝑀.𝑡𝑠(2.28) Using the definition of norm of the spaces 𝐸𝛼, we can obtain that 𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝑑𝑠𝐸𝛼𝐸𝛼=𝑡0𝑡𝑠𝐴𝑒(𝑡𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐸𝐸𝑑𝑠+sup𝜆>0𝑡0𝑡𝑠𝜆1𝛼𝐴𝑒𝜆𝐴𝐴𝑒(𝑡𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐸𝐸𝑑𝑠.(2.29) Using estimates (2.23) and (2.28), we get 𝐽3(𝑡)𝐸𝛼𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝑑𝑠𝐸𝛼𝐸𝛼𝐹𝐶(𝐸𝛼)𝑀3𝐹𝐶(𝐸𝛼).(2.30) Expanding 𝐺(𝑠), estimation of 𝐽4(𝑡) is as follows: 𝐽4(𝑡)𝐸𝛼𝑡0𝑄1𝑞,𝜌,𝑥,𝑥𝑤,𝑡𝜉,𝑥𝜋(𝑡𝜉)1/2𝐸𝛼+𝑑𝜉𝑡0𝑄2𝑞,𝑥,𝑥𝐷𝑡1/2𝑤𝜉,𝑥𝜋(𝑡𝜉)1/2𝐸𝛼𝑑𝜉+𝑡0𝐷𝑡1/2𝑤(𝜉,𝑥)𝜋(𝑡𝜉)1/2𝐸𝛼𝑑𝜉.(2.31) It is known that (see [34]) 𝑤𝐸𝛼𝐷𝑀𝑡1/2𝑤𝐸𝛼.(2.32) Since 𝑄1(𝑞,𝜌,𝑥,𝑥,𝑡) and 𝑄2(𝑞,𝑥,𝑥) are known functions, it is easy to see that 𝐽4(𝑡)𝐸𝛼𝑀4𝑞,𝜌,𝑥,𝑥,𝑇𝑡01𝜋(𝑡𝜉)1/2𝐷𝑡1/2𝑤(𝜉)𝐸𝛼𝑑𝜉.(2.33) Estimation of 𝐽5(𝑡) can be given similar to the estimation of 𝐽4(𝑡). By (2.23) and (2.32), 𝐽5(𝑡)𝐸𝛼𝑡0𝑡𝑠𝐴𝑒(𝜉𝑠)𝐴𝜋(𝑡𝜉)1/2𝑑𝜉𝐺(𝑠)𝑑𝑠𝑀5𝑞,𝜌,𝑥,𝑥,𝑇𝑡0𝐷𝑡1/2𝑤(𝑠)𝐸𝛼𝑑𝑠.(2.34) Finally combining estimates (2.25), (2.26), (2.30), (2.33), and (2.34), we get 𝐷𝑡1/2𝑤𝐸𝛼𝑀1𝐴𝜑𝐸𝛼+𝑀2+𝑀3𝐹𝐶(𝐸𝛼)+𝑡0𝑀4𝜋(𝑡𝜉)1/2+𝑀5𝐷𝑡1/2𝑤(𝑠)𝐸𝛼𝑑𝑠.(2.35) Using the Gronwall’s inequality, we can write 𝐷𝑡1/2𝑤𝐸𝛼𝑒𝑀6𝑀1𝐴𝜑𝐸𝛼+𝑀7𝐹𝐶(𝐸𝛼).(2.36) From the last estimate, we can obtain the estimate for 𝑤𝑡(𝑡) by using problem (2.17) and well-posedness of the Cauchy problem in 𝐶(𝐸𝛼) (see [35]). So the following theorem finishes the proof of Theorem 2.2.

Theorem 2.3 (see, [36]). For 0<𝛼<1/2, the norms of the spaces 𝐸𝛼(𝐶[0,𝜋],𝐴) and 𝐶2𝛼[0,𝜋] are equivalent.

3. Numerical Results

We have not been able to obtain a sharp estimate for the constants figuring in the stability inequalities. So we will provide the following results of numerical experiments of the following problem: 𝜕𝑢(𝑡,𝑥)=𝜕𝜕𝑡2𝑢(𝑡,𝑥)𝜕𝑥2𝑢(𝑡,𝑥)+𝐷𝑡1/21𝑢(𝑡,𝑥)+𝑝(𝑡)sin𝑥+𝑓(𝑡,𝑥),𝑓(𝑡,𝑥)=3𝑡𝜋𝑡1/2+2𝜋𝑡1/2],[],[],𝑢𝜋sin𝑥,𝑥(0,𝜋),𝑡(0,1𝑢(0,𝑥)=sin𝑥,𝑥0,𝜋𝑢(𝑡,0)=𝑢(𝑡,𝜋)=0,𝑡0,1𝑡,2=1𝑡.(3.1) The exact solution of the given problem is 𝑢(𝑡,𝑥)=(1𝑡)sin𝑥 and for the control parameter 𝑝(𝑡) is 1+𝑡.

3.1. The First Order of Accuracy Difference Scheme

For the approximate solution of the problem (3.1), the Rothe difference scheme 𝑢𝑘𝑛𝑢𝑛𝑘1𝜏=𝑢𝑘𝑛+12𝑢𝑘𝑛+𝑢𝑘𝑛12𝑢𝑘𝑛+𝐷𝜏1/2𝑢𝑘𝑛+𝑝𝑘𝑞𝑛𝑡+𝑓𝑘,𝑥𝑛,𝑓𝑡𝑘,𝑥𝑛=3𝑡𝑘1𝜋𝑡𝑘1/2+2𝜋𝑡𝑘1/2sin𝑥𝑛,𝑝𝑘𝑡=𝑝𝑘,𝑞𝑛=sin𝑥𝑛,𝑥𝑛=𝑛,𝑡𝑘𝑢=𝑘𝜏,1𝑘𝑁,1𝑛𝑀1,𝑀=𝜋,𝑁𝜏=1,0𝑛=sin𝑥𝑛𝑢,0𝑛𝑀,𝑘0=𝑢𝑘𝑀𝑢=0,0𝑘𝑁,𝑘𝑠𝑡=𝜌𝑘𝑡,𝜌𝑘=1𝑡𝑘𝜋,0𝑘𝑁,𝑠=,2(3.2) where 𝑥 denotes greatest integer less than 𝑥 is constructed. Throughout the paper, let us denote 𝜌𝑡𝑘=1𝑡𝑘,𝑞𝑛=sin𝑥𝑛,𝑡𝑘=𝑡𝑘,𝑥=𝑘𝜏,0𝑘𝑁,𝑁𝜏=1𝑛=𝑥𝑛,𝑓𝑡=𝑛,0𝑛𝑀1,𝑀=𝜋𝑘,𝑥𝑛=3𝑡𝑘1𝜋𝑡𝑘1/2+2𝜋𝑡𝑘1/2sin𝑥𝑛,𝐹𝑡𝑘,𝑥𝑛=𝜌𝑡𝑘𝑥sin𝑠𝑥sin𝑛+1𝑥2sin𝑛𝑥+sin𝑛12𝑥sin𝑛1𝜋𝑘𝑚=1Γ(𝑘𝑚+(1/2))𝜏(𝑘𝑚)!1/2𝑥sin𝑠𝑥sin𝑛𝑡+𝑓𝑘,𝑥𝑛.(3.3) We search the solution of (3.2) in the following form: 𝑢𝑘𝑛=𝜂𝑘𝑞𝑛+𝑤𝑘𝑛,(3.4) where 𝜂𝑘=𝑘𝑖=1𝑝𝑖𝜏,1𝑘𝑁,𝜂0=0.(3.5) Moreover for the interior grid point 𝑢𝑘𝑠, we have that 𝑢𝑘𝑠=𝜂𝑘𝑞𝑠+𝑤𝑘𝑠𝑡=𝜌𝑘.(3.6) From (3.4), (3.5), and the condition 𝑢𝑘𝑠=𝜌(𝑡𝑘), it follows that 𝜂𝑘=𝜌𝑡𝑘𝑤𝑘𝑠𝑞𝑠,𝑝(3.7)𝑘=𝜂𝑘𝜂𝑘1𝜏𝑢,1𝑘𝑁,(3.8)𝑘𝑛=𝜌𝑡𝑘𝑤𝑘𝑠𝑞𝑠𝑞𝑛+𝑤𝑘𝑛,0𝑘𝑁,0𝑛𝑀,(3.9) where 𝑤𝑘𝑛,0𝑘𝑁,0𝑛𝑀 is the solution of the difference scheme 𝑤𝑘𝑛𝑤𝑛𝑘1𝜏=𝑤𝑘𝑛+12𝑤𝑘𝑛+𝑤𝑘𝑛12𝑤𝑘𝑛𝑤𝑘𝑠𝑥sin𝑠𝑥sin𝑛+1𝑥2sin𝑛𝑥+sin𝑛12𝑥sin𝑛1𝜋𝑘𝑚=1Γ(𝑘𝑚+(1/2))𝑤(𝑘𝑚)!𝑚𝑠𝑤𝑠𝑚1𝑥sin𝑠𝜏1/2𝑥sin𝑛𝑤𝑚𝑛𝑤𝑛𝑚1𝜏1/2𝑡+𝐹𝑘,𝑥𝑛𝑤,1𝑘𝑁,1𝑛𝑀1,𝑘0=𝑤𝑘𝑀𝑤=0,0𝑘𝑁,0𝑛𝑥=sin𝑛,0𝑛𝑀.(3.10) First, applying the first order of accuracy difference scheme (3.10), we obtain (𝑀+1)×(𝑀+1) system of linear equations and we write them in the matrix form 𝐴𝑤𝑘+𝑘1𝑗=0𝐵𝑗𝑤𝑗=𝐷𝜑𝑘,1𝑘𝑁,𝑤0=𝑥sin𝑛𝑀𝑛=0,(3.11) where 𝐴=1000.0.000𝑥𝑦𝑥0.𝑧1.0000𝑥𝑦𝑥.𝑧2.000..........0000.𝑧𝑀1.𝑥𝑦𝑥0000.0.001(𝑀+1)×(𝑀+1),𝐵0=000.0.000𝑎0.𝑓1.0000𝑎.𝑓2.00........000.𝑓𝑠+𝑎.00........000.𝑓𝑀1.0𝑎000.0.00(𝑀+1)×(𝑀+1),𝐵𝑗=000.0.000𝑐0.𝑑1.0000𝑐.𝑑2.00........000.𝑑𝑠+𝑐.00........000.𝑑𝑀1.0𝑐000.0.00(𝑀+1)×(𝑀+1),(3.12) for any 𝑗=1,2,,𝑘2, and 𝐵𝑘1=000.0.000𝑣0.𝑐1.0000𝑣.𝑐2.00........000.𝑐𝑠+𝑣.00........000.𝑐𝑀1.0𝑣000.0.00(𝑀+1)×(𝑀+1).(3.13) Here, for any 𝑛=1,2,,𝑀1, 1𝑥=21,𝑦=𝜏+221+1+𝜏,𝑧𝑛=𝑥sin𝑛+1𝑥2sin𝑛𝑥+sin𝑛1𝑥sin𝑠2𝑥sin𝑛𝑥sin𝑠𝑥sin𝑛𝑥sin𝑠𝜏Γ,in(𝑠+1)thcolumn,𝑎=(𝑘1/2)𝜏𝜋(𝑘1)!,𝑓𝑛𝑥=sin𝑛Γ(𝑘1/2)𝑥𝜏𝜋sin𝑠,1(𝑘1)!𝑐=𝜏𝜋Γ(𝑘𝑗1/2)(𝑘𝑗1)!Γ(𝑘𝑗+1/2),𝑑(𝑘𝑗)!𝑛=𝑥sin𝑛𝑥𝜏𝜋sin𝑠Γ(𝑘𝑗1/2)(𝑘𝑗1)!Γ(𝑘𝑗+1/2),1(𝑘𝑗)!𝑣=𝜏1𝜏,𝑐𝑛=𝑥sin𝑛𝑥𝜏sin𝑠,𝑤𝑟=𝑤𝑟0𝑤𝑟𝑀(𝑀+1)×1𝜑forany𝑟=0,1,,𝑘,𝑘=0𝜙𝑘1𝜙𝑘𝑀10(𝑀+1)×1,𝜙𝑘𝑛=𝜌𝑡𝑘𝑥sin𝑠𝑥sin𝑛+1𝑥2sin𝑛𝑥+sin𝑛12𝜌𝑡𝑘𝑥sin𝑠𝑥sin𝑛1𝜋𝑘𝑚=1Γ(𝑘𝑚+1/2)𝜏(𝑘𝑚)!1/2𝑥sin𝑠𝑥sin𝑛𝑡+𝑓𝑘,𝑥𝑛,(3.14) and 𝐷 is (𝑀+1)×(𝑀+1) identity matrix. Using (3.11), we can obtain that 𝑤𝑘=𝐴1𝐷𝜑𝑘𝑘1𝑗=0𝐵𝑗𝑤𝑗,𝑘=1,2,,𝑁,𝑤0=sin𝑥𝑛𝑀𝑛=0.(3.15) To solve the resulting difference equations, we apply the method given in (3.15) step by step for 𝑘=1,2,,𝑁. For the evaluation of 𝑤𝑟,𝑟=2,3,,𝑁,𝑤𝑟1 is needed. It is obtained in the previous step. Then, the solution pairs (𝑢,𝑝) are obtained by using the last formulas (3.9) and (3.8).

3.2. The Second Order of Accuracy Difference Scheme

For the approximate solution of the problem (3.1), the Crank-Nicholson difference scheme 𝑢𝑘𝑛𝑢𝑛𝑘1𝜏=𝑢𝑘𝑛+12𝑢𝑘𝑛+𝑢𝑘𝑛122+𝑢𝑘1𝑛+12𝑢𝑛𝑘1+𝑢𝑘1𝑛122𝑢𝑘𝑛+𝑢𝑛𝑘12+𝑝𝑘+𝑝𝑘12𝑞𝑛+𝐷𝜏1/2𝑢𝑡𝑘𝜏2,𝑥𝑛𝑡+𝑓𝑘𝜏2,𝑥𝑛,𝑝𝑘𝑡=𝑝𝑘𝑢,1𝑘𝑁,1𝑛𝑀10𝑛𝑥=sin𝑛𝑢,0𝑛𝑀,𝑘0=𝑢𝑘𝑀𝑢=0,0𝑘𝑁,𝑘𝑠+𝑢𝑘𝑠+1𝑢𝑘𝑠𝑥𝑡𝑠=𝜌𝑘𝜋,0𝑘𝑁,𝑠=2(3.16) is constructed.

Here, Γ1𝑘𝑟+2=0𝑡𝑘𝑟+1/2𝑒𝑡𝑑𝑡.(3.17) Moreover, applying the second order of approximation formula for 𝐷𝑡1/2𝑢𝑡𝑘𝜏2=1Γ(1/2)𝑡𝑘0𝜏/2𝑡𝑘𝜏2𝑠1/2𝑢(𝑠)𝑑𝑠,(3.18) it is obtained (see [27]) 𝐷𝜏1/2𝑢=𝑑23𝑢0+𝑑23𝑢1+𝑑𝜏32𝑥sin𝑛,𝑘=1,2𝑑65𝑢0+𝑑65𝑢1+𝑑65𝑢2𝑑𝜏6𝑥10sin𝑛𝑑,𝑘=2,𝑘1𝑚=2(𝑘𝑚)𝑏1+𝑏2𝑢𝑚2+(2𝑚2𝑘1)𝑏12𝑏2𝑢𝑚1+(𝑘𝑚+1)𝑏1+𝑏2𝑢𝑚+𝑑62𝑢𝑘24𝑢𝑘1+5𝑢𝑘,3𝑘𝑁.(3.19) Here and throughout the paper, 𝐷𝑡1/2𝑢=𝐷𝜏1/2𝑢𝑡𝑘𝜏2,𝑥𝑛,2𝑑=𝜋𝜏,𝑏1=1𝑘𝑚+21𝑘𝑚2,𝑏21=31𝑘𝑚+23/21𝑘𝑚23/2.(3.20) We search the solution of (3.16) in the following form: 𝑢𝑘𝑛=𝜂𝑘𝑞𝑛+𝑤𝑘𝑛,(3.21) where 𝜂𝑘=𝑘𝑖=1𝑝𝑖+𝑝𝑖12𝜏,1𝑘𝑁,𝜂0=0.(3.22) We have that 𝑢𝑘𝑠+𝑢𝑘𝑠+1𝑢𝑘𝑠𝑥𝑠=𝜂𝑘𝑥1𝑠𝑞𝑠+𝑥𝑠𝑞𝑠+1+𝑥1𝑠𝑤𝑘𝑠+𝑥𝑠𝑤𝑘𝑠+1𝑡=𝜌𝑘.(3.23) Let us denote 𝑥𝑦=𝑠=𝑥𝑥,(3.24) where 0𝑦<1. Then, one can write 𝜂𝑘=𝜌𝑡𝑘(1𝑦)𝑤𝑘𝑠𝑦𝑤𝑘𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1.(3.25) So the values of (𝑝(𝑡𝑘)+𝑝(𝑡𝑘1))/2,1𝑘𝑁 can be obtained by the following formula: 𝑝𝑘+𝑝𝑘12=𝜌𝑡𝑘𝑡𝜌𝑘1𝑤/𝜏(1𝑦)𝑘𝑠𝑤𝑠𝑘1𝑤/𝜏𝑦𝑘𝑠+1𝑤𝑘1𝑠+1/𝜏(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1.(3.26) Let 𝑤𝑟 denote 𝑤𝑟=𝑤𝑟0𝑤𝑟𝑀(𝑀+1)×1for𝑟=0,1,,𝑁.(3.27) For 𝑘=1, one can show that 𝑤1 is the solution of the difference scheme 𝑤1𝑛𝑤0𝑛𝜏=𝑤1𝑛+12𝑤1𝑛+𝑤1𝑛122+𝑤0𝑛+12𝑤0𝑛+𝑤0𝑛122𝑤1𝑛+𝑤0𝑛2+𝑞𝑛+12𝑞𝑛+𝑞𝑛122𝑞𝑛2×𝜌𝑡1(1𝑦)𝑤1𝑠𝑦𝑤1𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝜌𝑡0(1𝑦)𝑤0𝑠𝑦𝑤0𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑23𝜌𝑡1(1𝑦)𝑤1𝑠𝑦𝑤1𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤1𝑛𝑑23𝜌𝑡0(1𝑦)𝑤0𝑠𝑦𝑤0𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤0𝑛+𝑑𝜏32𝑡𝑞(𝑛)+𝑓1𝜏2,𝑥𝑛𝑤,1𝑛𝑀1,10=𝑤1𝑀𝑤=0,0𝑛𝑥=sin𝑛,0𝑛𝑀.(3.28) We have the system of linear equations and we write them in the matrix form 𝐴1𝑤1+𝐵1𝑤0=𝐷𝜑1,(3.29) where 𝐴1=10000000𝑎𝑦1𝑎0𝑙1𝑐10000𝑎𝑦1𝑎𝑙2𝑐2000000𝑎𝑙𝑠+𝑎𝑐𝑠0000000𝑙𝑠+1+𝑦1𝑐𝑠+1+𝑎0000000𝑙𝑀1𝑐𝑀1𝑎𝑦1𝑎000000001(𝑀+1)×(𝑀+1),𝐵1=00000000𝑎𝑣1𝑎0𝑑1𝑒10000𝑎𝑣1𝑎𝑑2𝑒2000000𝑎𝑑𝑠+𝑎𝑒𝑠0000000𝑑𝑠+1+𝑣1𝑒𝑠+1+𝑎0000000𝑑𝑀1𝑒𝑀1𝑎𝑣1𝑎000000000(𝑀+1)×(𝑀+1).(3.30) Here, for any 𝑛=1,2,,𝑀1, 1𝑎=22,𝑦1=1𝜏+12+12𝑑23,𝑣1=1𝜏+12+12+𝑑23,𝑙𝑛=𝑞𝑛+12𝑞𝑛+𝑞𝑛1(1𝑦)22(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞𝑛(1𝑦)2(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑23𝑞𝑛,𝑐𝑛=𝑞𝑛+12𝑞𝑛+𝑞𝑛1𝑦22(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞𝑛𝑦2(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑23𝑞𝑛,𝑑𝑛=𝑞𝑛+12𝑞𝑛+𝑞𝑛1(1𝑦)22(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞𝑛(1𝑦)2(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑑23𝑞𝑛,𝑒𝑛=𝑞𝑛+12𝑞𝑛+𝑞𝑛1𝑦22(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞𝑛𝑦2(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑑23𝑞𝑛,𝜑1=0𝜙11𝜙1𝑀10(𝑀+1)×1,𝜙1𝑛=𝑥sin𝑛+1𝑥2sin𝑛𝑥+sin𝑛122𝑥sin𝑛2𝜌𝑡1𝑡+𝜌0(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑2𝑞𝑛3𝜌𝑡1𝑡𝜌0(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑𝜏32𝑞𝑛𝑡+𝑓1𝜏2,𝑥𝑛(3.31) and 𝐷 is (𝑀+1)×(𝑀+1) identity matrix. Using (3.29), we can obtain that 𝑤1=𝐴11𝐷𝜑1𝐵1𝑤0,𝑤0=sin𝑥𝑛𝑀𝑛=0.(3.32) For 𝑘=2,  𝑤2 is the solution of the difference scheme 𝑤2𝑛𝑤1𝑛𝜏=𝑤2𝑛+12𝑤2𝑛+𝑤2𝑛122𝑤2𝑛+𝑤1𝑛2+𝑤1𝑛+12𝑤1𝑛+𝑤1𝑛122+𝑞𝑛+12𝑞𝑛+𝑞𝑛122𝑞𝑛2𝜌𝑡2(1𝑦)𝑤2𝑠𝑦𝑤2𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝜌𝑡1(1𝑦)𝑤1𝑠𝑦𝑤1𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑65𝜌𝑡2(1𝑦)𝑤2𝑠𝑦𝑤2𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤2𝑛+𝑑65𝜌𝑡1(1𝑦)𝑤1𝑠𝑦𝑤1𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤1𝑛2𝑑65𝜌𝑡0(1𝑦)𝑤0𝑠𝑦𝑤0𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤0𝑛𝑑𝜏6𝑡10𝑞(𝑛)+𝑓2𝜏2,𝑥𝑛𝑤,1𝑛𝑀1,20=𝑤2𝑀𝑤=0,0𝑛𝑥=sin𝑛,0𝑛𝑀.(3.33) The system of linear equations given above can be written in the matrix form 𝐴2𝑤2+𝐵2𝑤1+𝐶2𝑤0=𝐷𝜑2,(3.34) where 𝐴2=10000000𝑎𝑦2𝑎0𝑔110000𝑎𝑦2𝑎𝑔220000000𝑔𝑠+𝑎𝑠0000000𝑔𝑠+1+𝑦2𝑠+1+𝑎0000000𝑔𝑀1𝑀1𝑎𝑦2𝑎000000001(𝑀+1)×(𝑀+1),𝐵2=00000000𝑎𝑣2𝑎0𝑔110000𝑎𝑣2𝑎𝑔220000000𝑔𝑠+𝑎𝑠0000000𝑔𝑠+1+𝑣2𝑠+1+𝑎0000000𝑔𝑀1𝑀1𝑎𝑣2𝑎000000000(𝑀+1)×(𝑀+1),𝐶2=000000000𝑧00𝑖1𝑗100000𝑧0𝑖2𝑗20000000𝑖𝑠+𝑧𝑗𝑠0000000𝑖𝑠+1𝑗𝑠+1+𝑧0000000𝑖𝑀1𝑗𝑀10𝑧0000000000(𝑀+1)×(𝑀+1).(3.35) Here, for any 𝑛=1,2,,𝑀1, 1𝑎=22,𝑦2=1𝜏+12+12𝑑65,𝑣21=𝜏+12+12𝑑65,𝑔𝑛=𝑞𝑛+12𝑞𝑛+𝑞𝑛1(1𝑦)22(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞𝑛(1𝑦)2(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑6𝑞𝑛(1𝑦)5(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1,in(𝑠+1)thcolumn,𝑛=𝑞𝑛+12𝑞𝑛+𝑞𝑛1𝑦22(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞𝑛𝑦2(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑6𝑞𝑛𝑦5(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑖,in(𝑠+2)thcolumn,𝑛=2𝑑6𝑞𝑛(1𝑦)5(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1,𝑗𝑛2=6𝑞𝑛𝑦5(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1,𝑧=2𝑑65,𝜑2=0𝜙21𝜙2𝑀10(𝑀+1)×1,𝜙2𝑛=𝑥sin𝑛+1𝑥2sin𝑛𝑥+sin𝑛122𝑥sin𝑛2+𝑑6𝑞𝑛5×𝜌𝑡2𝑡+𝜌1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+12𝑑6𝑞𝑛5𝜌𝑡0(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑑𝜏6𝑞10𝑛𝑡+𝑓2𝜏2,𝑥𝑛.(3.36) Using (3.34), we can obtain that 𝑤2=𝐴21𝐷𝜑2𝐵2𝑤1𝐶2𝑤0,𝑤0=sin𝑥𝑛𝑀𝑛=0.(3.37) For 3𝑘𝑁, we can obtain the following difference scheme: 𝑤𝑘𝑛𝑤𝑛𝑘1𝜏=𝑤𝑘𝑛+12𝑤𝑘𝑛+𝑤𝑘𝑛122+𝑤𝑘1𝑛+12𝑤𝑛𝑘1+𝑤𝑘1𝑛122𝑤𝑘𝑛+𝑤𝑛𝑘12+𝑞𝑛+12𝑞𝑛+𝑞𝑛122𝑞𝑛2×𝜌𝑡𝑘(1𝑦)𝑤𝑘𝑠𝑦𝑤𝑘𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞𝑛+12𝑞𝑛+𝑞𝑛122𝑞𝑛2×𝜌𝑡𝑘1(1𝑦)𝑤𝑠𝑘1𝑦𝑤𝑘1𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑑𝑘1𝑚=2(𝑘𝑚)𝑏1+𝑏2×𝜌𝑡𝑚2(1𝑦)𝑤𝑠𝑚2𝑦𝑤𝑚2𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞(𝑛)(2𝑚2𝑘1)𝑏12𝑏2×𝜌𝑡𝑚1(1𝑦)𝑤𝑠𝑚1𝑦𝑤𝑚1𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞(𝑛)(2𝑚2𝑘1)𝑏12𝑏2𝑤𝑛𝑚1×𝜌𝑡𝑚(1𝑦)𝑤𝑚𝑠𝑦𝑤𝑚𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1+𝑞(𝑛)(𝑘𝑚)𝑏1+𝑏2𝑤𝑛𝑚2+(𝑘𝑚1)𝑏1+𝑏2𝑤𝑚𝑛+(2𝑚2𝑘1)𝑏12𝑏2𝑤𝑛𝑚1𝑑62𝜌𝑡𝑘2(1𝑦)𝑤𝑠𝑘2𝑦𝑤𝑘2𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤𝑛𝑘24𝑑62𝜌𝑡𝑘1(1𝑦)𝑤𝑠𝑘1𝑦𝑤𝑘1𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤𝑛𝑘1+5𝑑62𝜌𝑡𝑘(1𝑦)𝑤𝑘𝑠𝑦𝑤𝑘𝑠+1(1𝑦)𝑞𝑠+𝑦𝑞𝑠+1𝑞(𝑛)+𝑤𝑘𝑛𝑡+𝑓𝑘𝜏2,𝑥𝑛,𝑤1𝑛𝑀1,𝑘0=𝑤𝑘𝑀𝑤=0,3𝑘𝑁,0𝑛𝑥=sin𝑛,0𝑛𝑀.(3.38) This system can be written in matrix form as 𝐴3𝑤𝑘+𝐵3𝑤𝑘1+𝐶3𝑤𝑘2+𝑘3𝑗=0𝐸𝑗𝑤𝑗=𝐷𝜑𝑘,(3.39) where 𝐴3=10000000𝑎𝑦3𝑎0𝑟1𝑚10000𝑎𝑦3𝑎𝑟2𝑚20000000𝑟𝑠+𝑎𝑚𝑠0000000𝑟𝑠+1+𝑦3𝑚𝑠+1+𝑎0000000𝑟𝑀1𝑠𝑀1𝑎𝑦3𝑎000000001(𝑀+1)×(𝑀+1),𝐵3=00000000𝑎𝑋𝑎0𝑌1𝑍10000𝑎𝑋𝑎𝑌2𝑍20000000𝑌𝑠+𝑎𝑍𝑠0000000𝑌𝑠+1+𝑋𝑍𝑠+1+𝑎0000000𝑌𝑀1𝑍𝑀1𝑎𝑋𝑎00000000(𝑀+1)×(𝑀+1),𝐶3=000000000𝑆00𝑅1𝑇100000𝑆0𝑅2𝑇20000000𝑅𝑠+𝑆𝑇𝑠0000000𝑅𝑠+1𝑇𝑠+1+𝑆0000000𝑅𝑀1𝑇𝑀10𝑆000000000(𝑀+1)×(𝑀+1),𝐸0=000000000𝐺00𝐻1𝐼100000𝐺0𝐻2𝐼20000000𝐻s+𝐺𝐼s0000000𝐻s+1𝐼s+1+𝐺0000000𝐻𝑀1𝐼𝑀10𝐺000000000(𝑀+1)×(𝑀+1),𝐸1=000000000𝐺100𝐽1𝐾100000𝐺10𝐽2𝐾20000000𝐽𝑠+𝐺1𝐾𝑠0000000𝐽𝑠+1𝐾𝑠+1+𝐺10000000𝐽𝑀1𝐾𝑀10𝐺1000000000(𝑀+1)×(𝑀+1),𝐸𝑗=000000000𝐺200𝑃1𝐿100000𝐺20𝑃2𝐿20000000𝑃𝑠+𝐺2𝐿𝑠0000000𝑃𝑠+1𝐿𝑠+1+𝐺20000000𝑃𝑀1𝐿𝑀10𝐺2000000000(𝑀+1)×(𝑀+1),(3.40) for 𝑗=2,3,𝑘3.

Here, for any 𝑛=1,2,,𝑀1, 1𝑎=222,𝑑=𝜋𝜏,𝑦3=1𝜏+1