This paper assesses the supply chain of Yartsagunbu (Caterpillar Fungus) in Darchula district of Nepal to identify who holds the power and how they gain power for management and marketing. We recorded two types of supply chain: (ⅰ) open supply chain, driven by open market, where the product is transported to Kathmandu before export to international market, and (ⅱ) close chain practiced by indigenous Shauka community following customary trade route to Tibet. The open chain is longer with higher number of actors compared to the close chain. This study observed that actors have intensive horizontal competition in the open chain to collect and purchase maximum quantity. Therefore, profit is disproportionately distributed to the actors in higher level of the supply chain. The profit is based on the price the actor receives, which is determined by their bargaining power. An actor's bargaining power is determined by the capital holding capacity, market information, risk appetite, networking and social ties. The study suggests that Government's interventions such as providing security to traders, access to finance, organizing auction and providing market information can help to increase the bargaining power of lower level actors. The study also suggests to minimize the disturbance in the collection site through limiting the collection permit and revising the revenue based on the market price.
1.
Introduction
It is known that fractional calculus is widely used, especially in engineering and science, and it is mainly described and realized through fractional differential equations (FDEs). For example: finance, biology, automatic control theory, mechanics, information theory and other fields have real applications. The Caputo-type fractional derivatives are divided into right Caputo derivative and left Caputo derivative. According to the research [1], by means of piecewise quadratic interpolation, a FODE high-order numerical scheme with consistent precision is given, and the unconditionally stable and consistent sharp numerical order is strict. A general high order method was used to solve FODE based on block-by-block method in [2]. In [3], by means of the L2 formula, a high-order compact finite difference method with respect to time FPDE is constructed. The high-order finite-difference methods with respect to time FPDE is introduced, moreover the first proof the stability of 3−α was presented in [4]. The finite-difference method introduced in [5] is a meshless generalized type and is applied to deal with 3D variable-order time FPDE. In [6], the Caputo fractional derivative is approximated by a fractional numerical differentiation method and uses a high-order numerical approach to deal with time-dependent FPDEs. In [7], in order to obtain the fully discrete form of the spatial fractional diffusion wave equation, the time-dependent second-order difference method and the spatial finite element method are proposed. According to the study [8], the time fractional diffusion equation in regard to the Mittag-Leffler matrix function is solved by the line method. In [9], a high-order numerical scheme is introduced to infinitely approximate the Riemann-Liouville fractional derivative. In [10], some numerical methods are proposed to solve fractional calculus, such as piecewise interpolation and Simpson's method. The paper [11] presented a fractional Taylor-type formula for the right Caputo fractional derivative, meanwhile the derivative is characterized by a fractional differential formula and a fractional integral remainder. In [12], the right Caputo fractional derivative is obtained by making use of the singularity-preserving spectrum configuration approach of fractional differential equations, the convergence analysis to obtain and the best error estimate is obtained. In [13], based on the derived operation matrix, using the Gauss-Lobatto quadrature formula is not only effective for studying fractional optimal control problems, but also for fractional order variational problems. At the same time, it is also verified that the Lagrange multiplier method is still valid for them. Discuss two second-order numerical differential equations of Caputo derivative operator, moreover, proving its unconditional convergence and stability with the help of discrete energy approach in [14]. In [15], some numerical approaches to fractional derivatives are constructed, which are Caputo-type derivatives with finite real-valued lower bounds and Riemann-Liouville-type derivatives with infinite lower bounds. In [16], based on cubic Lagrangian interpolation polynomials, proposed a high-order numerical method to approximate Caputo fractional derivatives. In [17], in order to determine the spatially related source terms in the opposite problem of the time-fractional diffusion equation, the time finite difference method and the spatially local discontinuous Galerkin method are used to construct a numerical scheme. A high-order form of the Caputo-type fractional convection-diffusion equation is constructed under the Dirichlet boundary condition in [18]. In [19], the numerical solution for the coupled nonlinear time-dependent FPDE with Caputo derivatives was given by the implicit trapezoidal method. In [20], the high-order scheme of Caputo fractional derivative approximation was developed and applied to solve the time-dependent FPDEs. With smooth conditions on the solution, piecewise polynomials were popularly used to solve fractional differential equations; refer to the literature [21,22,23,24,25] for details. For no-smooth solutions, one should introduce the non-uniform meshes for the fractional differential equations, and readers can refer to [26,27,28] for details.
In the current literature, we learn that there are limited high-order numerical schemes for studying right Caputo FDEs with uniform accuracy. Therefore, based on the idea of [1], this paper constructs a high-order numerical scheme with uniform precision convergence for right Caputo FDEs.
The main content of this paper is organized as follows: We adopt a high-order numerical scheme to solve the right Caputo FODE, and Section 2 analyzes the local truncation error and stability of this scheme in detail. In Section 3, we study time-dependent FPDE with right Caputo derivative, implementing discrete-time FPDE by finite-difference methods in time, and second-order central-difference methods in space. In Section 4, some numerical examples are given to verify efficient high-order numerical methods and support our theoretical findings. Finally, some conclusions from this work are drawn in Section 5.
2.
The high-order numerical scheme of the right Caputo FODE
2.1. Construction of the high-order numerical scheme
Consider the following FODE
where the initial condition is m(b)=m0, and m0 is a known constant, tDλbm(t) in (2.1) defined as the right Caputo fractional derivative of order λ which is given by
here Γ(⋅) means Gamma function in [29].
Now, we will construct a high-order scheme of (2.1) and divide the interval [a,b] into W uniform sub-intervals. Suppose a=t0<…<tq<…<tW=b, tq=qη, q=0,1,2,…,W, and η=b−aW is the step size, mq is the numerical solution of (2.1) at tq, and gq=g(tq,mq).
In order to discretize the right Caputo fractional derivative of (2.2), one can firstly determine the values of m(t) at tW−2,tW−1,tW, and employ quadratic Lagrange interpolation to m(t) on [tW−2, tW] as follows
where the quadratic Lagrangian interpolation basis function at points tW−2,tW−1 and tW are κi,W−2,i=0,1,2, which are defined as follows
For q=W−1,W−2, substituting (2.3) into (2.2), we can obtain
where the expression of Ei,0W−1, Ei,0W−2 are as follows
When q≤W−3, firstly in the interval of [tl−1,tl], the approximation solution of m(t) at points tl−1, tl, tl+1 is
where
Substituting (2.3) and (2.9) into (2.2), we have
where
and κi,W−2(t), ωi,l(t) are in (2.4) and (2.10), respectively.
In summary, the linear combination of ml approximates the right Caputo derivative tDλbm(tq). After calculation, it is found that each Ei,lq is proportional to η−λΓ(3−λ). After sorting (2.5),(2.6) and (2.11), the discrete Caputo derivative ηLλbmq be obtained as follows
where the values of all the coefficients ″E,F,G,H″ have been carefully calculated as follows
With the observe of the above coefficients, we will find that all the coefficients only depend on the constant of λ. With the help of (2.13), we have
Therefore, this sufficiently shows that the high-order numerical scheme corresponding to the above Eq (2.1) is
2.2. Error estimation
In order to analysis the local error of scheme (2.15), we bring in the following operator
where the values of all the coefficients ″E,F,G,H″ are defined by (2.14), which is replace mq with m(tq) in (2.13).
Theorem 2.1. Suppose m(t)∈C3[a,b],0<λ<1, let
we have
where Cm is a positive constant independent of η.
Proof. The following error estimates are estimated with the help of Taylor's theorem,
where ξW(ρ)∈[tW−2,tW], ξl(ρ)∈[tl−1,tl+1].
When q=W−1, according to (2.19), we get the following
where Cm is a positive constant independent of η.
When q=W−2, by (2.19), we have
When q≤W−3, according to (2.20), we obtain that
Next, we estimate each item at the right hand side of (2.23). By using (2.19), S1 is calculated as following
S2 and S3 are calculated as follows
Bringing (2.24) and (2.25) into (2.23), we can get
Combining (2.21), (2.22) with (2.26), Theorem 2.1 has been proved.
2.3. Stability analysis
We will now analyze the stability of the scheme (2.15), analogous to integer order differential equations, considering the problem
where θ>0 is a real number, first introducing the symbol,
Now rewrite the scheme (2.15) as follows, for q≤W−3,
When q=W−1,W−2, we have
Equations (2.29) and (2.30) are solved together, and to further simplify the representation, the coefficients are introduced, when q≤W−3,
Therefore, (2.29) is re-expressed as follows
Before analyzing the stability of (2.32), two auxiliary lemmas are given.
Lemma 2.2. For all 0<λ<1,q≤W−3, the coefficients in the problem (2.32) satisfy:
1) W∑s=q+1dqs=1;
2) dqs>0,s=W,W−1,…,q+4,q+3;
3) 0<dqq+1<43;
4) The symbol for dqq+2 can not be determined;
5) 14(dqq+1)2+dqq+2>0.
Proof. For a detailed proof, see Appendix A.
In Lemma 2.2, we find that the value of dqq+2 can be positive or negative. In order to prove that the high order scheme is absolutely stability for λ∈(0,1), let
then the Eq (2.32) can be re-expressed as follows
Now we denote
The equivalent form of (2.29) is
Lemma 2.3. For all 0<λ<1,q≤W−3, the coefficients of (2.36) satisfy
1) 0<τ<23;
2) ˉdqs>0,s=W,W−1,…,q+2;
3) τ+W−1∑s=q+2ˉdqs+ˉdqW≤1.
Proof. 1) According to 3) in Lemma 2.2 and (2.33), it is obviously provable.
2) When s=q+2,
According to 5) in Lemma 2.2, we have ˉdqq+2>0, thus
and by 2) in Lemma 2.2, τ>0, we have
3) Assume Qq=τ+W−1∑s=q+2ˉdqs+ˉdqW, for (2.34), we have
Further available,
According to 1), 2), 5) in Lemma 2.2 and (2.33), we get
Therefore, Qq=τ+W−1∑s=q+2ˉdqs+ˉdqW≤1. The proof of Lemma 2.3 is then completed.
Lemma 2.4. For all 0<λ<1, there is
Proof. For a detailed proof, see Appendix B.
Theorem 2.5. When g is defined by (2.27), the high-order numerical scheme to (2.1)is stable, and have
Proof. According to Lemma 2.4, we get
From (2.35) and ms=ˉms+τms+1=…=ˉms+τˉms+1+…+τW−1−sˉmW−1+τW−smW, bying using (2.38) and 1) in Lemma 2.3, we have
To sum up, Theorem 2.5 certificate is completed.
3.
High-order numerical schemes of FPDEs
3.1. High-order numerical schemes of 1D time FPDEs
Consider the following time FPDEs
where λ represents the order of a fractional derivative in regard to time t, m0(y) is a known function, and the relevant definition of tDλbm(y,t) is shown below
To construct a time-dependent high-order numerical scheme, [a,b] is divided into W equal parts. Marking η=b−aW,tq=qη,q=0,1,2,…,W, the interval [c,d] is divided into M equal parts, set Δy=d−cM,yi=c+iΔy,i=0,1,…,M. And let the numerical solution of (3.1) at (yi,tq) as mi,q. Next, construct the discrete scheme of (3.1) at time t, the derivation process is similar to the derivation of (2.13), then we have
where the values of the coefficients ″E,F,G,H″ are defined by (2.14).
Therefore, the semi-discrete scheme of (3.1) is
First analyze the truncation error of the above scheme, similar to Theorem 2.1, we have the following lemma.
Lemma 3.1. Suppose that m(y,t) is the solution of (3.1), and it has fourth-order continuous partial derivative in regard to time variable t, then
where Cm is a constant independent of η.
On the discretization of ∂2m(y,t)∂y2, for the fixed tq, use the idea of the central second-order difference scheme to discretize, and the method is as follows
Lemma 3.2. For a fixed time t, use the second-order central difference method to discretize ∂2m(y,t)∂y2, the scheme (3.6) is known that it has the second-order accuracy under suitable conditions.
Proof. Its detail proof can be found in reference [30].
We make use of the finite difference scheme (3.4) in the discretization of time and (3.6) in the discretization of space, and the high-order numerical scheme of the FPDE (3.1) is as follows
where mi,q represents the numerical solution of m(yi,tq), gi,q represents g(yi,tq), and i=1,2,…,M−1. The error of the discrete scheme is given below, here, we first bring in two defined the operators
Theorem 3.3. Assume m(y,t) is the solution of the Eq (3.1) and regarding the time variable t and space variable y both have 4th-order continuous partial derivatives, so
where Cm is a constant independent of η and Δy.
Proof. According to Lemmas 3.1, 3.2 and (3.8), we have already proved above, we are able to immediately gain
Theorem 3.3 is then proved.
3.2. High-order numerical schemes of 2D FPDEs
Consider the following FPDEs
where λ represents the order of the fractional derivative in regard to time t, m0(y,z) is a known function, and the relevant definition of tDλbm(y,z,t) is shown below
Similar to (3.7), let Δy=Δz=d−cM,zl=c+lΔz,l=0,1,…,M, we have
where mi,l,q represents the numerical solution of m(yi,zl,tq), gi,l,q represents g(yi,zl,tq),
where i,l=1,2,…,M−1. The error of the scheme is given below, we first define the two operators
where δ2y,δ2z are defined as following for i,l=1,2,…,M−1,
The proof is given below that similar to Theorem 3.3.
Theorem 3.4. Assume m(y,z,t) is the solution of the Eq (3.11) and has 4th-order continuous partial derivatives with respect to t,y,z, then
where Cm is a constant independent of η,Δy,Δz.
4.
Numerical results
In this section, three computational examples will be used to verify our conclusions and methods in the previous section.
Example 4.1. Case (1). the exact solution is smooth. For the problem (2.1), let the initial value m0=0, and
It can be verified that the exact solutions of the above two cases are all m(t)=(1−t)4. We take a=0, b=1, the step size is η=1W, W=2ˉn, where ˉn=4,…,10. The following two cases for (4.1) and (4.2) use several different values of λ, and we will calculate the convergence order with the help of ln(e2η/eη)/ln2, where eη=maxq|m(tq)−mq|.
Firstly, when the function g is defined by (4.1) which g is a linear case of m, it can be seen from Table 1 that when λ takes 0.3,0.5 and 0.7, their convergence orders tend to be 2.7,2.5 and 2.3, respectively. In Table 2, we take λ→0 or 1, that is λ=0.01 and 0.99, we find convergence orders close to 2.99 and 2.01, respectively. Therefore, from Tables 1 and 2, when λ take different values, even when it tends to 0 or 1, the convergence rate is still 3−λ,λ∈(0,1), and this is basically consistent with our theoretical expected results.
Secondly, when g is defined by (4.2) which g is the nonlinear case of m, it can be seen from Table 3 that when λ takes 0.2,0.5, and 0.8, the convergence order tends to 2.8,2.5, and 2.2, respectively. Similar to the above Table 2, λ=0.01 and 0.99 are also taken in Table 4, and the convergence order tends to 2.99 and 2.01, respectively. In short, the results in Tables 3 and 4 still verify that our order of convergence is 3−λ. Therefore, in case (1) of the Example 4.1, when 0<λ<1, by taking different values of λ, we find that the obtained convergence orders are all 3−λ, which again confirms the theoretical prediction in Theorem 2.1.
Case (2). The exact solution is non-smooth. Consider the problem (2.1), and m0=0 with an exact analytical solution:
It can be checked that the corresponding right hand side
The choices of the fractional order are now also taken as λ=0.3,0.5 and 0.7. Other settings are the same as previous case (1) in the Example 4.1.
The results are given in Table 5, from which we can obtain that when 0<λ<1, the convergence order is close to λ. The reason lies that the exact solution function has singularity at b=1 and is non-smooth on [0,1].
Case (3). We consider an exact analytical solution and the corresponding right hand side are (4.3) and (4.4), respectively.
In this case, we choose graded mesh is ti=1−(1−i/W)β with a grading parameter β>1 based on the idea of [27] for i=0,1,2,…,W. The choices of the fractional order are now also taken as λ=0.3,0.5,0.7 and W=8,16,32,64,128,256,512 and β=3.
The results are given in Table 6, from which we can obtain that when 0<λ<1, the convergence order is close to min(βλ,3−λ) which is result of the Lemma 8 in [27].
Example 4.2. To demonstrate the validity of the algorithm and observe the approximation of the numerical solution to the exact solution, we solve the Eq (3.1) by means of the scheme (3.7), and the corresponding right-hand member as following
The analytic solution was verified as m(y,t)=(1−t)4sin2πy.
We take a=0,b=1, the interval of the space be [0, 1], and set eη,Δy=maxi,q|mi,q−m(yi,tq)|. We start by looking at spatial accuracy. To prevent the effect of time-dependent error on spatial accuracy, we need to fix the time step to be small enough. Let W=10,000 in Table 7, by comparing the error of λ, Δy under different values and the order of convergence. When 0<λ<1, the spatial accuracy is second-order convergence, this result is accord with the theoretical analysis obtained in Theorem 3.3.
Next, we check the time-dependent convergence rate. In Table 8, we also list the value of eη,Δy and the corresponding order, when λ, η takes a series of different values, where Δy=O(η3−λ2) is taken. When λ takes 0.4,0.6, and 0.8, the convergence order tends to 2.6,2.4, and 2.2, respectively, this can show that the convergence order in the time is 3−λ.
Example 4.3. In order to further test the feasibility of the algorithm, we use the scheme (3.13) to solve the Eq (3.11), the maximum error is eη,Δy,Δz=maxi,l,q|mi,l,q−m(yi,zl,tq)|, and take the right side of the equation as follows
and its exact solution is m(y,z,t)=(1−t)4sin2πysin2πz.
To study the convergence rate of the space, in Table 9, let W=6000, the domain of space be [0,1]×[0,1]. By observing the error and convergence order with λ and Δy=Δz choosing different values, we find that its spatial accuracy is second-order convergence.
In Table 10, where Δy=Δz=O(η3−λ2) is taken, studying the convergence order in time. When λ takes 0.3,0.5, and 0.8, the convergence order tends to 2.7,2.5, and 2.2, respectively, the numerical results reveal the theoretical analysis is accord with the numerical results, the convergence order is 3−λ.
5.
Conclusions
In this paper, the high-order numerical scheme of the right Caputo FODE is constructed, and the local truncation error and stability are analyzed in detail based on the idea and methods of [1]. Secondly, the high-order scheme is used to solve the time FPDEs. Thirdly, three numerical examples are used to verify the validity of our conclusions that the optimal convergence rate of time is 3−λ,λ∈(0,1) with uniform accuracy. Due to the length limitation of the paper, we only give the local error estimate for FPDEs. The convergence analysis can be directly obtained by using the method of stability of FODE in this paper and the ideas of [4]. In the future, We will study higher order numerical schemes with low smoothness based on the good idea of [26,27,28] by using the non-uniform mesh and we expect that the constructed efficient high-order scheme can be applied to the fractional order optimal control problem and topology optimization nonlocal problem of composites plate based on the idea of [31,32,33].
Acknowledgments
The work was partially supported by National Natural Science Foundation of China (11961009 and 11901135), Guizhou Provincial Science and Technology Projects ([2020]1Y015), and Natural Science Research Project of Department of Education of Guizhou Province (QJJ2022015 and QJJ2022047). The authors thank the anonymous referees for their valuable suggestions to improve the quality of this work significantly.
Conflict of interest
The authors declare there is no conflicts of interest.
Appendix
A. The proof of Lemma 2.2
Proof. 1) can also be checked by a direct calculation using the definition of dqs and summing them for all s=q+1,…,W. For example, for q=W−3, by using (2.14) and (2.28), we have
The case of other values of q can be verified similarly.
2) Because
when q≤W−3, using Taylor's formula, we have
where
it is calculated that when r≥6, ax is a convergent alternating series with positive first term, so
when q=W−1,W−2, using a method similar to (A.1), it can be shown that, dqW−1>0,dqW−1>0,
In summary,
3) Calculated from (2.31)
thus
let φ(λ)=20−5λ+(3λ−18)21−λ, by calculation we get
therefore φ′(λ) is monotonically decreasing and 0<φ′(1)<φ′(λ)<φ′(0), so φ(λ) is monotonically increasing and φ(0)<φ(λ)<φ(1)=0, that is, there is
Let ˜ϕ(λ)=2−λ(2λ−12)−3λ+12, a tedious but routine calculation gives ˜ϕ(λ)>0, Therefore dqq+1>0, in short, 0<dqq+1<43.
4) Because
where 4−λ>0,λ∈(0,1), so the sign of dqq+2 is determined by ϕ1(λ). It can be known by calculation, ϕ′1(0)>0,ϕ′1(1)<0, ϕ1(λ) increase at first and then decrease, and ϕ1(0)=0,ϕ1(1)=−1<0. Therefore, when λ∈(0,1), dqq+2 have positive and negative values. Therefore, The symbol for dqq+2 can not be determined.
5) Because
where φ1(λ)=(3λ−12)(4−λ)+3(λ2−10λ+24)22−λ+4(12λ−λ2−32)31−λ+(6−λ)241−λ. According to Lagrange's mean value theorem, 41−λ>43+λ⋅31−λ, we get
where
so ∑+∞n=0fn is a convergent alternating series, that is, 0<∑+∞n=0fn<f0. Thus
where
It can be obtained by calculation
Because of λ∈(0,1), so λ3<λ2,λ6>0, then
By Calculating, we get ˉφ″3(λ)<0, so ˉφ′3(λ) is monotonically decreasing, so ˉφ′3(1)<ˉφ′3(λ)<ˉφ′3(0), a similar method can be used to find the first derivative of ˉφ′3(λ), because ˉφ′3(0)>0,ˉφ′3(1)<0, thus ˉφ′3(λ) changes from positive to negative, so ˉφ3(λ) first increases and then decreases, and ˉφ3(0)=0,ˉφ3(1)=191>0, so ˉφ3(λ)>0,φ1(λ)>φ2(λ)>ˉφ3(λ)>0. To sum up,
The proof of Lemma 2.2 is completed.
B. The proof of Lemma 2.4
Proof. Bringing (2.27) into scheme (2.13), when q=W−1,W−2,
and β0=ηλθ the relevant value of ″E″ is detailed in (2.14). It can be obtained by calculation,
where
So, when q=W−1,
After detailed calculation, it can be obtained,
and
With the above calculation, for all λ∈(0,1), when q=W−1, we have
Therefore, we obtain
When q=W−2, it can be proved that
holds by a similar method as q=W−1.
The following proves that when q=W−3, bring in (2.32) there is
Due to τ=12[3−21−λ(6−λ)4−λ],dW−3W−2=3−31−λ(λ+4)4−λ, and τ≠12dW−3W−2, according to (2.36) we get
where ˉdW−3W−2=dW−3W−2−τ,ˉdW−3W−1=τˉdW−3W−2+dW−3W−1,ˉdW−3W=τˉdW−3W−1+dW−3W.
Next, we will prove ˉdW−3W−2≥0,ˉdW−3W−1≥0,ˉdW−3W≥0. By carefully calculation,
Because dW−3W−1=14−λ31−λ(4λ+4)−3, then
By a tedious but routine calculation gives ˉdW−3W−1>0.
Because dW−3W=24−λ(4−λ2−λ2⋅32−λ)>0, so ˉdW−3W=τˉdW−3W−1+dW−3W>0.
Next calculate ˉdW−3W−2+ˉdW−3W−1+ˉdW−3W≤1. By carefully calculate, we have
According to 1) in Lemma 2.2, we get
In summary,
When q=W−3, multiply 2ˉmW−3 on both sides of (B.4) to have
For the left side of the above equation
By using (B.5), we have
By carefully calculation, we have
It is easy to check as follows
From conditions (B.2), (B.3), (B.8) and 3) in Lemma 2.3, (B.6) becomes following as
Therefore, when q=W−3, (2.37) holds.
When q=W−4, there is
According to Lemma 2.3, we have τ+W−1∑s=W−2ˉdqs+ˉdqW≤1,ˉdW−4s>0,s=W−2,W−1,W, multiply both sides of (B.10) by 2ˉmW−4 at the same time, we have
Due to 0<τ<23, using (B.2), (B.3) and (B.9), then
When q≤W−5, using a similar method above, multiply 2ˉmq on both sides of (2.36), and after sorting, we can get,
by mathematical induction, Lemma 2.3, we obtain,
In summary, the proof of Lemma 2.4 is completed.