Literature DB >> 27777859

A new multi-step technique with differential transform method for analytical solution of some nonlinear variable delay differential equations.

Brahim Benhammouda1, Hector Vazquez-Leal2.   

Abstract

This work presents an analytical solution of some nonlinear delay differential equations (DDEs) with variable delays. Such DDEs are difficult to treat numerically and cannot be solved by existing general purpose codes. A new method of steps combined with the differential transform method (DTM) is proposed as a powerful tool to solve these DDEs. This method reduces the DDEs to ordinary differential equations that are then solved by the DTM. Furthermore, we show that the solutions can be improved by Laplace-Padé resummation method. Two examples are presented to show the efficiency of the proposed technique. The main advantage of this technique is that it possesses a simple procedure based on a few straight forward steps and can be combined with any analytical method, other than the DTM, like the homotopy perturbation method.

Entities:  

Keywords:  Differential transform; Laplace–Padé method; Nonlinear variable delay differential equations

Year:  2016        PMID: 27777859      PMCID: PMC5052250          DOI: 10.1186/s40064-016-3386-8

Source DB:  PubMed          Journal:  Springerplus        ISSN: 2193-1801


Background

Differential equations are relevant tools to model a wide variety of physical phenomena across all areas of applied sciences and engineering. Analytical techniques are applied to find the exact solutions of some cases of differential equations; nevertheless, when the differential equations are nonlinear there are no general techniques of solutions. Therefore, numerical methods are important tools to study and understand the quantitative behavior of the nonlinear differential equations with unknown exact solutions. However, such methods can exhibit numerical instabilities, oscillations or false equilibrium states, among others (Gumel 2002, 2003; de Markus and Mickens 1999). This means that the numerical solution may not correspond to the real solution of the original problem. This situation is further aggravated for some types of differential equations like differential-algebraic equations, fractional differential equations and time delay differential equations (DDEs), among others (Ford and Wulf 2000; Engelborghs et al. 2000; Ascher and Petzold 1998; Campbell et al. 2008; Wanner and Hairer 1998). Time DDEs appear in propagation and transport phenomena, population dynamics, bioscience problems, neural network model, control problems, electrical networks containing lossless transmission lines and economical systems where decisions and effect are separated by some time intervals, among many other applications (Martín and García 2002a, b; Aiello and Freedman 1990; Gourley and Kuang 2004a, b; Kuang 1993; Li et al. 2006; Li and Jiang 2013; Zhang and Zhang 2013). Given the importance of this kind of equations, some numerical methods have been developed to solve them; among them one can mention: variable multi-step methods (Martín and García 2002a, b), Chebyshev polynomials for pantograph differential equation (Sedaghat et al. 2012), a spectral Galerkin method for nonlinear delay convection-diffusion reaction equations (Liu and Zhang 2015), power series method (Benhammouda et al. 2014a) and the differential transform method (DTM) (Karako and Bereketoglu 2009). In this work, we propose some case studies of nonlinear DDEs with variable time delays. For such case studies, there are no known numerical methods available. Therefore, we propose a multi-step method with a modified version of the DTM (Zhou 1986; Keskin et al. 2007; Benhammouda et al. 2014b; Benhammouda and Vazquez-Leal 2015; Biazar and Eslami 2010; Chen and Liu 1998; Ayaz 2004; Kangalgil and Ayaz 2009; Kanth and Aruna 2009; Arikoglu and Ozkol 2007; Chang and Chang 2008; Kanth and Aruna 2008; Lal and Ahlawat 2015; Odibat et al. 2010; El-Zahar 2013; Gökdoğan et al. 2012) in order to find the solutions of variable time DDEs. In the literature, there are some works reporting variable time delays. For instance, Taylor series method is applied to solve time-dependent stochastic DDEs (Milošević and Jovanović 2011). In addition, a study of asymptotic neutral type differential equations with variable time delay is given in Skvortsova (2015). Besides, a study of asymptotic behavior of first order differential equations with variable delay is presented in Dix (2005). In Ding et al. (2010), some new conditions for the boundness and stability by means of the contraction mapping principle are given for nonlinear scalar DDEs with variable delays. The issue of uniqueness of variable DDEs is investigated in Winston (1970), Eloe et al. (2005), Liu and Clements (2002) and Luo et al. (2013). Finally, a research for the existence of attractors for differential equations with a variable delay is presented in Graef and Qian (2000) and Caraballo et al. (2001). Nonetheless, in the present study, we propose different types of variable delays in terms of algebraic expressions of the time. What is more, given that the approximate solutions of the DTM are power series solutions, we propose to extend the domain of convergence by using a combined scheme of a multi-step technique and the Laplace–Padé resummation method (Vazquez-Leal and Guerrero 2014; Filobello-Nino et al. 2013; Jiao et al. 2002; Sweilam and Khader 2009; Momani et al. 2009; Khan and Faraz 2011; Momani and Ertürk 2008; Tsai and Chen 2010; Ebaid 2011). This paper is organized as follows: in “Differential transform method” section, we introduce the basic concept of the DTM. Then, in “Multi-step technique and DTM for nonlinear variable DDEs” section, the proposed multi-step technique with the use of the DTM to deal with variable time DDEs is presented. “Padé approximant” and “Laplace–Padé resummation method” sections are devoted to present the basic concepts of Padé and Laplace–Padé resummation methods. In “Case studies” section, the type of delay differential equations with variable delays under study is presented and solved using the proposed technique. Finally, discussion and conclusions are given in “Discussion” and “Conclusions” sections, respectively.

Differential transform method

For convenience of the reader, we will give a review of the differential transform method (DTM) (Zhou 1986; Keskin et al. 2007; Benhammouda et al. 2014b; Benhammouda and Vazquez-Leal 2015; Biazar and Eslami 2010; Chen and Liu 1998; Ayaz 2004; Kangalgil and Ayaz 2009; Kanth and Aruna 2009; Arikoglu and Ozkol 2007; Chang and Chang 2008; Kanth and Aruna 2008; Lal and Ahlawat 2015; Odibat et al. 2010; El-Zahar 2013; Gökdoğan et al. 2012). We will also describe the DTM to solve systems of ordinary differential equations.

Definition 1

(Zhou 1986; Keskin et al. 2007) If a function u(t) is analytical with respect to t in the domain of interest, thenis the transformed function of u(t).

Definition 2

(Zhou 1986; Keskin et al. 2007) The differential inverse transforms of the set is defined bySubstituting (1) into (2), we deduce thatFrom the above definitions, it is easy to see that the concept of the DTM is obtained from the power series expansion. To illustrate the application of the proposed DTM to solve systems of ordinary differential equations, we consider the nonlinear systemwhere is a nonlinear smooth function. System (4) is supplied with some initial conditionsThe DTM establishes that the solution of (4) can be written aswhere are unknowns to be determined by the DTM. Applying the DTM to the initial conditions (5) and system (4) respectively, we obtain the transformed initial conditionsand the recursion systemwhere is the differential transforms of . Using (7) and (8), the unknowns can be determined. Then, the differential inverse transformation of the set of values gives the approximate solutionwhere m is the approximation order of the solution. The exact solution of problem (4)–(5) is then given by (6). If U(k) and V(k) are the differential transforms of u(t) and v(t) respectively, then the main operations of the DTM are shown in Table 1.
Table 1

Main operations of the DTM

FunctionDifferential transform
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha u(t)\pm \beta v(t) $$\end{document}αu(t)±βv(t) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha U( k) \pm \beta V( k) $$\end{document}αU(k)±βV(k)
u(t)v(t) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\sum _{r=0}^{k}}}U\left( r\right) V\left( k-r\right) $$\end{document}r=0kUrVk-r
u(t)v(t) w(t) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\sum _{r=0}^{k}\sum _{l=0}^{r}}}U\left( l\right) V\left( r-l\right) W\left( k-r\right) $$\end{document}r=0kl=0rUlVr-lWk-r
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{d^{n}}{dt^{n}}[u(t) ]$$\end{document}dndtn[u(t)] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( k+1\right) \ldots \left( k+n\right) U\left( k+n\right) $$\end{document}k+1k+nUk+n
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^{\lambda t} $$\end{document}eλt \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{\lambda ^{k}e^{\lambda t_{0}}}{k!}$$\end{document}λkeλt0k!
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin \left( \omega t\right) $$\end{document}sinωt \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{\omega ^{k}}{k!}\sin \left( \omega t_{0}+\dfrac{\pi k}{2}\right) $$\end{document}ωkk!sinωt0+πk2
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos \left( \omega t\right) $$\end{document}cosωt \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{\omega ^{k}}{k!}\cos \left( \omega t_{0}+\dfrac{\pi k}{2}\right) $$\end{document}ωkk!cosωt0+πk2
Main operations of the DTM The process of the DTM can be described as: Apply the differential transform to the initial conditions (5). Apply the differential transform to the differential system (4) to obtain a recursion system for the unknowns Use the transformed initial conditions (7) and the recursion system (8) to determine the unknowns Use the differential inverse transform formula (9) to obtain an approximate solution for initial-value problem (4)–(5). The solutions series obtained from the DTM may have limited regions of convergence. Therefore, we propose to apply the Laplace–Padé resummation method (Vazquez-Leal and Guerrero 2014; Filobello-Nino et al. 2013; Jiao et al. 2002; Sweilam and Khader 2009; Momani et al. 2009; Khan and Faraz 2011; Momani and Ertürk 2008; Tsai and Chen 2010; Ebaid 2011) to the DTM truncated series to enlarge the convergence region as depicted in “Padé approximant” and “Laplace–Padé resummation method” sections.

Multi-step technique and DTM for nonlinear variable DDEs

The type of nonlinear variable delay differential equations (DDEs) which are considered here are given by where the initial condition is given and and n is a positive integer. For the nonlinear lag function the delay function is and is such that . The solution of (10)–(11) is assumed to exist, unique and analytical. To solve (10)–(11), we start by finding the first interval of approximation. This is achieved by solving the inequalities and which lead to the interval , where Therefore, the first sub-problem to solve is given by Since for then substituting (13) into (12) reduces problem (12)–(13) to the following initial-value problem where (15) is obtained from (11). To solve (14)–(15), the differential transform is applied to it to get the recursion where is the differential transform of . Using (16)–(17), the unknowns , for can be determined. Then, the differential inverse transformation of the set of values gives the approximate solutionNow if , then the above process is terminated and is the approximate solution. Otherwise, if , then the inequalities and are solved to get . Then, the solution is extended by solving the following sub-problem Since for then substituting (20) into (19) reduces problem (19)–(20) to the following initial-value problem where (22) is obtained from (20). To solve (21)–(22), the differential transform is applied to it to get the recursion where is the differential transforms of . Using (23)–(24), the unknowns can be determined. Then, the differential inverse transformation of the set of values gives the approximate solutionFinally, continuing this process, one can extend the domain of the solution to the desired interval.

Padé approximant

Given an analytical function u(t) with Maclaurin’s expansionThe Padé approximant to u(t) of order [L, M] which we denote by is defined by (Baker and Graves-Morris 1996)where we considered , and the numerator and denominator have no common factors. The numerator and the denominator in (27) are constructed so that u(t) and and their derivatives agree at up to . That isFrom (28), we haveFrom (29), we get the following algebraic linear systemsandFrom (30), we calculate first all the coefficients . Then, the coefficients are determined from (31). Note that for a fixed value of , the error (28) is smallest when the numerator and denominator of (27) have the same degree or when the numerator has degree one higher than the denominator.

Laplace–Padé resummation method

Several approximate methods provide power series solutions (polynomial). Nevertheless, sometimes, this type of solutions lack large domains of convergence. Therefore, Laplace–Padé resummation method (Vazquez-Leal and Guerrero 2014; Filobello-Nino et al. 2013; Jiao et al. 2002; Sweilam and Khader 2009; Momani et al. 2009; Khan and Faraz 2011; Momani and Ertürk 2008; Tsai and Chen 2010; Ebaid 2011) is used in literature to enlarge the domain of convergence of solutions or inclusive to find the exact solutions. The Laplace–Padé method can be summarized as follows: First, Laplace transformation is applied to power series (9). Next, s is substituted by 1/t in the resulting equation. After that, the transformed series is convert into a meromorphic function by forming its Padé approximant of order [N/M]. N and M are arbitrarily chosen, but they should be smaller than the order of the power series. In this step, the Padé approximant extends the domain of the truncated series solution to obtain better accuracy and convergence. Then, t is substituted by 1/s. Finally, by using the inverse Laplace s transformation, the exact or an approximate solution is obtained.

Case studies

In this section, we will demonstrate the effectiveness and accuracy of the multi-step technique proposed in “Multi-step technique and DTM for nonlinear variable DDEs” section with the differential transform method (DTM) (Zhou 1986; Keskin et al. 2007; Benhammouda et al. 2014b; Benhammouda and Vazquez-Leal 2015; Biazar and Eslami 2010; Chen and Liu 1998; Ayaz 2004; Kangalgil and Ayaz 2009; Kanth and Aruna 2009; Arikoglu and Ozkol 2007; Chang and Chang 2008; Kanth and Aruna 2008; Lal and Ahlawat 2015; Odibat et al. 2010; El-Zahar 2013; Gökdoğan et al. 2012) through two nonlinear variable delay differential equations (DDEs).

Example 1

Consider the following nonlinear variable delay differential equation where and For this problem the lag function is and the variable delay is . Following the procedure described in “Multi-step technique and DTM for nonlinear variable DDEs” section, problem (32)–(33) is solved step by step. The first interval of approximation is determined by solving the inequalities and which lead to the interval where Therefore, the first sub-problem to solve is given by Since for then substituting (35) into (34) reduces problem (34)–(35) into the following initial-value problem where condition (37) is obtained from (33). To solve (36)–(37), the differential transform is applied to it to obtain the following recursion From recursion (38)–(39), the following values are obtained From these values, an eighth-order approximate solution is constructedwhere Note that one can use more terms in the above series solutions. Nevertheless, sometimes, this may not increase the accuracy of the solution for large intervals. Therefore, we use the Laplace–Padé resummation method presented in “Laplace–Padé resummation method” section to enlarge the domain of convergence of solutions or to find the exact solutions as follows. Applying Laplace transform to yieldsFor simplicity let , thenAll of the -Padé approximants of (44) with and and yieldNow since , we obtain in terms of s as followsFinally, applying the inverse Laplace transform to (46) gives the following approximate solutionwhich is the exact solution of delay equation (34)–(35). In a similar manner, the second interval of approximation is determined. The conditions and yield the interval where Therefore, we consider the sub-problem: Since for then substituting (49) into (48) reduces problem (48)–(49) to the following initial-value problem where condition (51) is obtained from (49). To solve (50)–(51), the differential transform is applied to it to obtain the following recursion From recursion (52)–(53), the following values are obtained From these values, an eighth-order approximate solution is constructedwhere Applying Laplace transform to yieldsFor simplicity let , thenAll of the [L/M] -Padé approximants of (58) with L and M and yieldNow since , we obtain in terms of s as followsFinally, applying the inverse Laplace transform to (60) gives the following approximate solutionwhich is the exact solution of delay equation (48)–(49). Finally, combining (47) and (61), the exact solution for the nonlinear variable delay problem (32)–(33) is obtained.

Example 2

Consider the following nonlinear variable delay differential equation where and For this problem the lag function is and the delay is As for example 1, problem (62)–(63) is solved step by step. The first interval of approximation is determined by solving the inequalities and which lead to the interval where Thus, the first sub-problem to solve is given by: Since for then substituting (65) into (64) reduces problem (64)–(65) to the following initial-value problem where conditions (67) are obtained from (63). To solve (66)–(67), the differential transform is applied to it to obtain the following recursionwhere and and where represents the differential transform of at From recursion (68), the following values are obtained From these values, an eighth-order approximate solution is constructedwhere To enlarge the domain of convergence of the approximate solution or to find the exact solution, the Laplace–Padé resummation method is used as in the previous example. Applying Laplace transform to yieldsFor simplicity let , thenAll of the -Padé approximants of (73) with and and yieldNow since , we obtain in terms of s as followsFinally, applying the inverse Laplace transform to (75) gives the following approximate solutionThus, the exact solution of the nonlinear variable delay equation (64)–(65) is . In a similar manner, the second interval of approximation is determined. This interval is obtained by solving the inequalities and which lead to the interval , where Thus, the second sub-problem to solve is given by: Since for then substituting (78) into (77) reduces problem (77)–(78) to the following initial-value problem where conditions (80) are obtained from (78). To solve (79)–(80), the differential transform is applied to it to obtain the following recursionwhere and and where represents the differential transform of at From recursion (81), the following values are obtained From these values, an eighth-order approximate solution is constructedwhere Applying Laplace transform to yieldsFor simplicity let , thenAll of the -Padé approximants of (86) with and and yieldNow since , we obtain in terms of s as followsFinally, applying the inverse Laplace transform to (88) gives the following approximate solutionThus, the exact solution of the nonlinear variable delay equation (77)–(78) is Finally, combining (76) and (89), the exact solution of the nonlinear variable delay problem (62)–(63) is obtained.

Discussion

Variable delays differential equations (DDEs) with arbitrary types of nonlinear functions for the time delay is an open area of research that require new numerical and analytical methods in order to deal with their solution. Therefore, from the examples above, it is worthwhile to remark that the multi-step technique proposed in this work combined with modified differential transform method (DTM) (Zhou 1986; Keskin et al. 2007; Benhammouda et al. 2014b; Benhammouda and Vazquez-Leal 2015; Biazar and Eslami 2010; Chen and Liu 1998; Ayaz 2004; Kangalgil and Ayaz 2009; Kanth and Aruna 2009; Arikoglu and Ozkol 2007; Chang and Chang 2008; Kanth and Aruna 2008; Lal and Ahlawat 2015; Odibat et al. 2010; El-Zahar 2013; Gökdoğan et al. 2012) based on Laplace–Padé resummation method (Vazquez-Leal and Guerrero 2014; Filobello-Nino et al. 2013; Jiao et al. 2002; Sweilam and Khader 2009; Momani et al. 2009; Khan and Faraz 2011; Momani and Ertürk 2008; Tsai and Chen 2010; Ebaid 2011) was able to obtain the exact solutions of nonlinear DDEs with variable delays. It is important to highlight that this multi-step technique was able to obtain the exact solutions for both case studies within the given intervals. The straight forward procedure was able to deal with different algebraic time delays as quadratic and cubic term highlighting the malleability of the technique presented to solve nonlinear DDEs with variable delays. As far as the knowledge of authors goes, no numerical or analytical approximation methods to solve the type of case studies of this work have been reported in the literature. Hence, we propose as a proof of concept, to solve problems with known exact solutions. Therefore, this multi-step technique in combination with the DTM and Laplace–Padé resummation was able to obtain the exact solutions within the given intervals. However, it is still pending for future work in order to treat DDEs with unknown exact solutions. For such hard to solve problems, we will measure the error of approximation using the mean square residual (MSR) error as in Benhammouda and Vazquez-Leal (2015). Finally, further research is required to extend this proposed methodology to other delay functions and systems of DDEs, delay differential-algebraic equations and delay partial differential equations with variable time delays.

Conclusions

This article deals with the solution of nonlinear DDEs with variable delays using a new multi-step method and the DTM (Zhou 1986; Keskin et al. 2007; Benhammouda et al. 2014b; Benhammouda and Vazquez-Leal 2015; Biazar and Eslami 2010; Chen and Liu 1998; Ayaz 2004; Kangalgil and Ayaz 2009; Kanth and Aruna 2009; Arikoglu and Ozkol 2007; Chang and Chang 2008; Kanth and Aruna 2008; Lal and Ahlawat 2015; Odibat et al. 2010; El-Zahar 2013; Gökdoğan et al. 2012). This technique was tested on two nonlinear problems. The results obtained show that the technique can be applied to solve these types of equations efficiently obtaining the exact solution. On the one hand, it is important to highlight that this kind of variable delay present series issues for traditional numerical and analytical methods, and on the other, the DTM in combination with Laplace–Padé resummation method (Vazquez-Leal and Guerrero 2014; Filobello-Nino et al. 2013; Jiao et al. 2002; Sweilam and Khader 2009; Momani et al. 2009; Khan and Faraz 2011; Momani and Ertürk 2008; Tsai and Chen 2010; Ebaid 2011) was able to obtain the exact solutions of nonlinear DDEs with variable delays. Future work is necessary to involve into the application of the proposed methodology for the solution of nonlinear DDEs with variable delays and unknown exact solutions. Other types of problems, like delay differential-algebraic equations and delay partial differential equations with variable time delays, will also be considered.
  3 in total

Review 1.  A stage structured predator-prey model and its dependence on maturation delay and death rate.

Authors:  Stephen A Gourley; Yang Kuang
Journal:  J Math Biol       Date:  2004-05-31       Impact factor: 2.259

2.  A time-delay model of single-species growth with stage structure.

Authors:  W G Aiello; H I Freedman
Journal:  Math Biosci       Date:  1990-10       Impact factor: 2.144

3.  Modeling the glucose-insulin regulatory system and ultradian insulin secretory oscillations with two explicit time delays.

Authors:  Jiaxu Li; Yang Kuang; Clinton C Mason
Journal:  J Theor Biol       Date:  2006-05-18       Impact factor: 2.691

  3 in total

北京卡尤迪生物科技股份有限公司 © 2022-2023.