999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Nonlinear Problems via a Convergence Accelerated Decomposition Method of Adomian

2021-04-27 10:28:46MustafaTurkyilmazoglu

Mustafa Turkyilmazoglu

1Department of Mathematics,Hacettepe University,Ankara,06532,Turkey

2Department of Medical Research,China Medical University Hospital,China Medical University,Taichung,40447,Taiwan

ABSTRACT The present paper is devoted to the convergence control and accelerating the traditional Decomposition Method of Adomian (ADM).By means of perturbing the initial or early terms of the Adomian iterates by adding a parameterized term,containing an embedded parameter,new modified ADM is constructed.The optimal value of this parameter is later determined via squared residual minimizing the error.The failure of the classical ADM is also prevented by a suitable value of the embedded parameter,particularly beneficial for the Duan–Rach modification of the ADM incorporating all the boundaries into the formulation.With the presented squared residual error analysis,there is noneed to check out the results against the numerical ones,as usually has to be done in the traditional ADM studies to convince the readers that the results are indeed converged to the realistic solutions.Physical examples selected from the recent application of ADM demonstrate the validity,accuracy and power of the presented novel approach in this paper.Hence,the highly nonlinear equations arising from engineering applications can be safely treated by the outlined method for which the classical ADM may fail or be slow to converge.

KEYWORDS Nonlinear equations;Adomian decomposition method;modification;convergence acceleration

1 Introduction

Researchers prefer an easily accessible and user friendly method requiring less computational labor while accurately approximating highly nonlinear equations resulting from mathematical modeling of real-life phenomena.The Adomian decomposition method (ADM) is one such popular technique capable of dealing with the prevailing nonlinearities by means of Adomian polynomials [1,2].A modification of the classical ADM is proposed within the current study based on the recent publications [3,4] successfully generating fast convergent ADM series solutions with as small Adomian polynomials as possible in the solution series.

A quick literature survey exhibits that ADM has been applied to many nonlinear equations [5].To classify some of the recent bibliography,algebraic equations were contained within the references [6,7].The ordinary differential equations were dealt within the citations [8–14].The articles [15–18] covered the efforts to partial differential equations.Mathematical analysis of the convergence of ADM to certain nonlinear equations was fulfilled in the publications [19–21].It can be successfully used to gain correct physical parameters domain [22].A traffic model was also very recently treated in [23] via the Adomian method.The publications by [24,25] present investigation of some nonlinear problems via different numerical approaches.

It is now well-known that an inadequate arrangement of the classical ADM series may lead to non-convergent solutions or solutions with a poor convergence rate.To avoid these shortcomings,a parameter is generally inserted at the leading term of the Adomian series and later it is subtracted at the first order term not to break down the equation structure.This procedure was pursued by the recent publications [11,17,18].However,how a proper value of the inserted parameter will be determined was not mentioned in these references.Instead,a randomly chosen value was assigned to it.A variety of modifications were also offered in the articles [26–31].A successful formulation of the ADM was made in the recent work of [3] which was named as the optimal ADM.Further applications of the homotopy analytic approximate method may be found in the literature [32–35].

The motivation of the current work is,benefiting from the idea in [11],to devise a method that greatly improves the mathematical property of classical ADM.Within this aim,a reorganization of the ADM series is proposed by altering the early terms so that they incorporate extra controllable terms.The reason of such a treatment is to get a rapidly converging ADM solutions with the least Adomian polynomials.In place of randomly selecting,an optimum value of the introduced parameter is later determined through error on the grounds of total residual.With this value at our disposal,there is no doubt that the ADM method is convergent to the true solution in a most rapid way,not demanding a verification of the ADM solutions by numerical ones.The failure of the classical ADM in the usual form or in the Duan–Rach formulation is also prevented by a suitable value of the embedded parameter.The present approach can also extend the region of convergence of the traditional method.Examples of physical value are provided to justify and validate the given procedure.

2 Traditional Decomposition Method of Adomian

The usual steps of traditional ADM can be inferred from the aforementioned citations.The methodology in brief is such that under an invertible linear operatorLand a forcing functionf,it is desired to approximate the functionuhaving the nonlinearityN(u)and satisfying the general nonlinear equation

with the initial and/or boundary restrictions

Having inverted (1) under the restrictions (2) generally leads to

wheregis due to the conditions in (2).Ifuis a single scalar parameter like (1) representing an algebraic equation,then there is no suchgin (3),whereas,in the case of a variableu,L?1denotes an integral operator giving rise togin (3).Then (3) is a mixed Volterra-Fredholm type equation so-called as the Duan–Rach formulation in the recent literature,see for instance [13,14].

The subsequent series decompositions ofuandN(u)

in whichAn’s are the classical Adomian polynomials,are later substituted into (3).The solutionuof (1) is finally generated from the recurrence relation

As a result,by means of the relations from (5),an approximate series solution of orderMis obtained as

which serves for practical purposes.

In general,the procedure in (5) yields convergent ADM series solutions,see for instance [19–21].If not,to achieve convergent solutions or for computational conveniences some modifications in the termsuiin (5) are implemented as in the articles [11,12],without a proper mathematical evidence and support.

3 A Modified Decomposition Method of Adomian

To overcome the divergence of classical ADM or to speed up the convergence rate of the ADM series,the leading order termu0in (5) (which is in compliance with the previous implementations,in for instance [11]) or some of the early terms,callue,0 ≤e≤Me,in (5) will be modified so that the modification will account for the change in the ADM by incorporating parameterized terms(h)withhan embedded unknown parameter [3].To exemplify,consider the modified version of (5) in the manner

The following conditions for parameterized terms in the new algorithm (7) should be added

so that it can be reduced to the traditional ADM forh=0.

It is remarked that there is no a unique way of selecting the(h)terms,but they must be as simple as possible in line with the simplicity inspiration of the ADM method itself.For instance,constants depending onhor simple power functions (or simple integrable functions)in combination withhcan be used,refer to the illustrations belove for more details.We should remark that in the absence of tilde terms (orh) (7) dully conforms with (5).Having inserted such tilde terms in the modified ADM (7),the question now arises,how to determine the correct or proper value ofh? So far,in ADM applications of this kind,only random or trial values are preassigned,see for instance [11] amongst others.On the other hand,we may either observe the least change in the plot of some physical quantities,the method is so-called as the constanth-level curves [36] giving rise to rough estimates forh,or we outline the following rigorous algorithm to determine the best or optimalh.

Algorithm.Consider the squared residual error corresponding to (1) defined by

where eitherL2(Ω)orL(Ω)norms are employed.Moreover,uin (8) is owing to the modified ADM from (7).The optimum value ofhis such that it minimizes the squared residual error (8)at the approximation levelM.

As a consequence,the above Algorithm will generate the best value ofhwhich will ensure the convergence of ADM series solution (7) in a fastest rate of convergence.The minimization task of(8) may be fulfilled by means of contemporary softwares,such as MAPPLE or MATHEMATICA.

4 Applications

Potential applications of the introduced ADM in (7) are given here.To control the error,we use the norm

with the exactueand ADM solutionu.

4.1 An Algebraic Equation

As stated by Adomian [1] the classical ADM method (5) fails to result in a convergent solution of

for the solutionu=?0.73205080757.On the other hand,when the new ADM is built via

Fig.1 displaysh-level curves at selected truncation ordersM.The intervalh∈(?0.35,?0.25)is observed to yield convergent ADM solutions,which excludes the traditional ADM withh=0 [1].

Through the residual minimization at the approximation levelM=8,h=?0.2679492 is obtained as the optimum.The history and why this value is the best for the convergence control,as compared to the failure of classical ADM can be visualized in Tab.1.

Figure 1:Convergence control parameter h regarding (11)

Table 1:Convergence of modified ADM (12) for Eq.(11) with different h.Paranthesis is for the absolute error (10)

The convergent solution of (11) with the new modified ADM (12) at the approximation levelM=8 is found to be

Steve went home from school, thoughtful, that afternoon. Walking into the house, he took one look around. Both parents were passed out, in various stages of undress, and the stench was overpowering! He, quickly, gathered up his camping gear, a jar of peanut butter, a loaf of bread, a bottle of water, and this time...his schoolbooks. Grim faced and determined7, he headed for the woods.

for which the optimumhis tabulated in Tab.1.

4.2 Equation Involving Integral

Consider the equation given in [30]

The modified ADM method (7) for the current integral problem is adopted as

We findh∈(?1500,1500)for the convergence interval in Fig.2,refer also to [3].A comparison and convergence accelerating feature of the present ADM can be visualized from the Tab.2.It is noticed from Tab.2 that even the 4th-order modified ADM (16) is able to produce the solution accurate of order 10?10,which is adequate in practical purposes.

Figure 2:Convergence control parameter h regarding (15)

Table 2:Errors in (15) regarding (14).Parenthesis denotes the optimum values of h

4.3 A Fin with Porosity Feature

As taken from [14],a porous fin can be modelled via

In (17),temperature along the fin isu,andsandβare physical parameters.

The modified ADM algorithm (7) here is

At the selected valuess=5 andβ=1,the history of convergence is tabulated in Tab.3.It is seen how useful the modified ADM over the classical one by more than doubling the convergence rate at the same number of iterations.The CPU times evaluated by MATHEMATICA through the residual are shown in Tab.3.The list proves the improved efficiency of the current method over the classical ADM.

Table 3:Residual errors of classical ADM and modified ADM (18) for Eq.(17) with s=5 and β=1.Parenthesis denote the optimum values of h

Table 4:Numerical and modified ADM results regarding (17) [3]

Instead of the modification of ADM in (18),we may use the Duan–Rach formulation involving no unknown parameters within it except the embedded parameterh

Choosings=β=1,Fig.3 showsh?level curves,indicating that the prescription ofhmust be made within the range(0.39,0.49)to ensure the convergence of the modified ADM 19.The exact values ofu(0)andu′(1)are respectively,0.522738093570 and 1.13937891581.

Figure 3:Convergence control parameter h regarding (17)

To demonstrate the power of the modified ADM (19),Tab.5 shows the squared residual error(9) from both the novel and classical ADM.It is unfortunate to observe that the Duan–Rach formulation (19) withh=0 fails to converge,however,the optimum embedded parameterhinsures that the modified ADM is convergent for the present physical problem,even if the convergence is not as fast as the modified formulation in (18).

Table 5:Values offor (17) from the classical and modified ADM methods

Table 5:Values offor (17) from the classical and modified ADM methods

M √Res(0) √Res(h) h 4 0.91×101 9.379592752×10?3 0.5717083134 8 5.90×102 8.485113566×10?5 0.4466738418 12 1.90×105 1.712138085×10?5 0.4205569783 16 9.50×107 8.888688579×10?6 0.4108495377 20 5.68×1010 3.872408211×10?6 0.4061899784

4.4 Gelfand Equation

The Gelfand equation [5] involves exponential nonlinearity [8]

withu′(0)=?0.4636325917 [3].

In line with the publications [5,8] whenh=0,the modified ADM is

u0(η)=0,

whereAn(h,η)are the polynomials of Adomian foreu(η),encompassing all boundary conditions in (20) within the Duan–Rach approach,and we have

Fig.4 shows the predicted convergence control parameters.WithM=12,an optimum value for the embedding parameterhis found to be?0.01274 from the Algorithm in (9).We find that the residual error is 1.6835117258×10?15withu′(0)=?0.4636325917 from the present approach.On the other hand,and error of 4.813740934×10?9withu′(0)=?0.4636325899 is calculated from the traditional ADM.Our value is given by

The success of the present modified ADM (21) is thus obvious.

Figure 4:Convergence control parameter h regarding (20)

4.5 Electrostatic Cantilever Micro-Electromechanical System

The beam-type electrostatic actuators for the nonlinear cantilever micro-electro mechanical systems are modelled by the fourth-order boundary value problem from [11]

To comply with the Duan–Rach Adomian decomposition method in [11],the present modified ADM is

where the Adomian polynomialsAn(h,η)are due to the negative-power nonlinearities in (23),see[11] for more details.

For the fixed parametersK=3,αK=0.2,β=0.5 andγ=0.25,Fig.5 produced from (24)helps us guess the proper values of embedded parameterh.Small values ofhare seen to be adequate to get convergent ADM solutions.

Figure 5:Convergence control parameter h regarding (23)

Fig.6 demonstrates different approximation levelsM,and it signifies toh=0.1045421730 as the optimumhwhenM=10.With this optimum value of the embedding parameter,the squared residual error for the current problem is=3.2727623195×10?9,whereas the classical residual error withh=0 is=8.7520447833×10?5.This implies that much more Adomian series terms are required for the classical ADM to reach the accuracy of the modified ADM here.Moreover,from [11] the physical values areu′′(0)=?0.56764 andu′′′(0)=1.07585 evaluated with 12 Adomian polynomials.The present values areu′′(0)=?0.5676684138 andu′′′(0)=1.0758898885 evaluated with only 10 Adomian polynomials,correct to 9 decimal places.Hence,the advantage of the present modification is clear for the present highly nonlinear physical problem.We should remark that with the present Algorithm,the validation of the ADM results against the numerical ones is no longer a prerequisite.Also,there is no need to evaluate error remainder functions,nor to consult to Pade-approximates to increase the accuracy,both of which as implemented in [11].

Figure 6:Error regarding (23) for various M

The convergence accelerating feature of the present modified ADM (24) as compared to the classical ADM is better visualized from the Tab.6.Table also shows the comparable CPU times.

Table 6:Convergence history of modified ADM vs.classical ADM for Eq.(23)

4.6 Electrostatic Cantilever Nano-Electromechanical System

Nonlinear model for the electrostatic double cantilever nano-electromechanical system in the case of Casimir force (K=4) is given by [11]

We adopt the subsequent modified ADM,that conforms to the classical ADM (h=0) given in [15]

where the Adomian polynomialsAn(h,η)are due to the negative-power nonlinearities in (25),see[11] for more details.

For the specific parametersαK=1,β=1.5 andγ=0.5,Fig.7 displays the constanth-level curves drawn at the approximation levelM=8.Similar to the previous example,very small values ofhare seen to be adequate to get convergent ADM.

In order to evaluate the performance of modified ADM over the classical one,Tab.7 shows the unknown physical quantitiesu′′(0)=u′′(1)andu′′′(1)at several truncation ordersM.The faster convergence rate of the present ADM is apparent.

To illustrate,the analytical formula computed via the present algorithm (26) atM=4 for the value ofu′′(0)is given by

which is of almost nine degree of accuracy as seen from Tab.7.

Table 7:Values of u′′(0) and u′′′(1) for (25) from the classical and modified ADM methods

4.7 Lane–Emden Equation

We consider the Lane–Emden type boundary value problem from [9]

that models the oxygen diffusion in a spherical cell with Michaelis–Menten uptake kinetics.We take into account the subsequent constants to comply with the literature [9]

k=4,r=0.76129,s=0.03119.

The modified ADM that is offered for the present problem is then

which conforms with the classical ADM of [9] in the limith→0 and the Adomian polynomialsAn(h,η)are to account for the nonlinearity.The interest is to determine the physical value ofu(0)=a.

We present Tab.8 to demonstrate the performance of the modified ADM (29) versus the classical ADM.The expected practical accuracy is met at lower Adomian series approximations via the modified method.

Table 8:Values of u(0)=a and for (28) from the classical and modified ADM methods

Table 8:Values of u(0)=a and for (28) from the classical and modified ADM methods

M √Res(0) √Res(h) a (h=0) a h 2 5.05920434×10?5 1.82258949×10?5 0.896870943665 0.896877048532?0.0289267765 4 2.24898539×10?7 2.18219055×10?8 0.896877026028 0.896877046709?0.0321993310 6 1.12229262×10?9 2.83931728×10?11 0.896877046624 0.896877046709?0.0335012655

4.8 The Fluid Flow of Jeffery–Hamel

The Jeffery–Hamel fluid flow problem is modelled via [22]

withα,ReandHaare the physical parameters [13].

Following the successful Duan–Rach ADM formulation of the problem (30) in [17],we propose the following modified version

whereAn(h,η)=?(4?Ha)α2u′n(η)?are the Adomian polynomials.

For the diverging channel,considering the specific parametersα=50,Re=50 andHa=1000 to be in line with [13],the constanth-level curves are depicted in Fig.8 at the approximation orderM=8.Accordingly,we expect the convergency of the modified ADM (31) in the vicinity of zero.

The performance of modified ADM (31) is next measured by computing the centerline velocityu(0.5)(numerical value is 0.764064240111) at different approximation levelsMas shown in Tab.9.It is observed that 10 digits of accuracy is quickly reached by the present ADM,whereas the classical ADM falls behind.Hence,even though it was not clearly mentioned in [13] (see Tab.1 therein),the accuracy of order 10?8as obtained via the classical ADM demands at least 15–20 Adomian polynomials,whereas only 6 Adomian polynomials are sufficient to gain the same accuracy with the present modification.

Figure 8:Convergence control parameter h regarding (30)

Table 9:Values of u(0.5) and for (30) from the classical and modified ADM methods

Table 9:Values of u(0.5) and for (30) from the classical and modified ADM methods

M √Res(0) √Res(h) u(0.5) (h=0) u(0.5) h 2 6.92050730×10?2 1.13520536×10?2 0.763699786377 0.764097626541?0.0180187189 4 1.94819917×10?3 1.87302552×10?5 0.764051346961 0.764064279681?0.0179816107 6 7.53651632×10?5 3.95453415×10?8 0.764063696627 0.764064239975?0.0178807933 8 3.34939742×10?6 9.72265613×10?11 0.764064214850 0.764064240111?0.0178880579

4.9 Squeezing Two Parallel Plates

The flow squeezed between two parallel plates are modelled by the nonlinear equations [22]

see [14] for the flow parameters.

In accordance with the Duan–Rach ADM formulation of the physical problem (32) in [14],we set the modified ADM in the form

where

are the Adomian polynomials.

To make a comparison with the classical ADM in [14],we set the parametersS=1,φ=0.02 and concentrate on Cu-Water nanofluid withρs=8933 andρf=997.1.The corresponding constanth-level curves are plotted in Fig.9.It appears that the convergence of the modified ADM(32) is guaranteed ifhis selected in the overlapping interval(?1,1).

Figure 9:Convergence control parameter h regarding (32)

The effects of iterative numberMon the skin frictionu′′(1)are next demonstrated in Tab.10.The advantage of the modified ADM (33) with optimum values of embedding parameterhis clearly observed against the classical ADM withh=0.It appears that the accuracy of classical ADM of O(10?8)as displayed in Tab.2 of [18] clearly requires more Adomian series terms (seeM=8 in Tab.10) as compared to the less terms needed in the modified ADM here.

Table 10:Values of u′′(1) and for (32) from the classical and modified ADM methods

Table 10:Values of u′′(1) and for (32) from the classical and modified ADM methods

M √Res(0) √Res(h) u′′(1) (h=0) u′′(1) h 2 2.67740258×10?1 6.68835495×10?2?3.6764267174?3.6917310584?0.0991340307 4 1.61057990×10?2 7.92481375×10?4?3.6952425904?3.6965216730?0.1097576248 6 1.09939458×10?3 6.61056361×10?6?3.6964553942?3.6965506329?0.1111585030 8 8.10446170×10?5 7.71479681×10?8?3.6965431988?3.6965505168?0.1112323962 10 6.28790630×10?6 1.09803715×10?9?3.6965499276?3.6965505117?0.1112270087

The following fourth-order modified ADM series solution for the skin friction may serve good to the purpose of engineering applications if not high accuracy is required

4.10 Nonlinear Oscillator Problem

Let us consider the nonlinear oscillator Duffing problem (see [36] (Chapter 5) and [4])

which involves a cubic nonlinearity.

The improved ADM can be given via

with the Adomian termsAn(t)in (35).

The classical Adomian method withh=0 in (36) is not convergent,whereas withh=0.68981924,the residual error becomes=0.02589439,and the convergency is satisfied as revealed in Fig.10 for the domain of definition.

Figure 10:Convergence control parameter h regarding (35)

4.11 Diffusion Equation

Let us consider the nonlinear diffusion equation,see [18] and [3]

for which [18] presents an exact solution.

The form of modified ADM for the partial differential equation (37) is

where

An(h,x,t)

are the Adomian polynomials corresponding to the function?u?u3+(u2ux)x.

Fig.11 shows the constanth-level curves at the approximation levelM=10,indicating a very large range of embedding parameterh.

Figure 11:Convergence control parameter h regarding (37)

Actually at this truncation of the modified ADM series,it is obtained

Defining the squared residual error for (37) as

Tab.11 tabulates how the modified ADM has smaller residual errors.

Table 11:The residual errors for (37) from the classical and modified ADM methods

Table 11:The residual errors for (37) from the classical and modified ADM methods

M √Res(0) √Res(h) h 2 3.067520559×10?1 1.126425884×10?2 0.78155339805 4 1.905329025×10?2 2.920732957×10?4 0.84725848564 6 5.284431922×10?4 4.431860700×10?6 0.88282078470 8 8.251967252×10?6 4.360598101×10?8 0.90500266099 10 8.249541273×10?8 2.997208214×10?10 0.92013994524

4.12 Burger’s Equation

The final example is known as Burger’s equation [3]

with an exact solutionu(x,t)=.

Figure 12:Convergence control parameter h regarding (41)

The form of modified ADM for the partial differential equation (40) is

where

An(h,x,t)are the Adomian polynomials corresponding to the functionuxx?uux.

The traditional ADM withh=0 in the domain 0 ≤x≤1,0 ≤t≤1 is noticed to fail to produce a convergent series as also clear in Fig.12,becauseh=0 is not lying in the overlapping region.

With the definition

Tab.12 justifies the success of the present modified ADM over the classical divergent one,both in terms of accuracy and computational cost.

Table 12:The residual errorsfor (41) from the classical and modified ADM methods

Table 12:The residual errorsfor (41) from the classical and modified ADM methods

M √Res(0) CPU √Res(h) (hopt) CPU 5 1.75×102 0.012 4.674547981×10?2 (2.2068296567) 0.039 10 1.67×105 0.053 4.082713247×10?3 (2.3147042279) 0.078 15 1.39×108 0.101 6.755704561×10?4 (2.3159080855) 0.145 20 1.24×1011 0.242 2.741368603×10?5 (2.3820449852) 0.289

5 Concluding Remarks

The aim of the present work is to present superiority over the well-known Adomian decomposition method (ADM) often employed in the recent literature to analytically approximate solutions to highly nonlinear algebraic and differential equations of some real physical motions.Within this aim,a reformulation of the ADM is targeted to prevent first the failure and then convergence acceleration of the classical Adomian polynomials.

To accomplish the objective,the classical ADM is modified by inserting some simple parameterized terms into the early iterates involving an embedded parameter to control and pacing the convergence of the generated ADM series.In order to determine the best suitable value or the optimum value of this parameter,squared residual minimizing of the governing equation is proposed.This enables us to overcome the divergence of the classical ADM,and more importantly,there is no need to check out the results against the numerical ones,as usually has to be done in traditional ADM studies,since the optimum embedded parameter obtained is an insurance for ADM series convergence in a most rapid manner.

Physical examples selected from the recent application of ADM demonstrate the validity,accuracy and power of the present approach in terms of generating the convergent solution within the least number of iterations.In particular,the Duan-Rach modification of the ADM incorporating all the boundaries mostly used in the recent ADM applications takes great benefit of the present proposal,otherwise there is always the inevitable danger that it may lead to non physical solutions.The present approach successfully extends the convergence interval of the studied problem.In conclusion,the present formulation of ADM offers a promising tool to treat more strongly nonlinear equations/systems of real life phenomena.

Funding Statement:The author received no specific funding for this study.

Conflicts of Interest:The author declares that he has no conflicts of interest to report regarding the present study.

主站蜘蛛池模板: 青青青伊人色综合久久| 又大又硬又爽免费视频| 无码中文AⅤ在线观看| 91在线免费公开视频| 精品国产一二三区| 美女被操91视频| 欧美成人午夜在线全部免费| 精品国产欧美精品v| 亚洲A∨无码精品午夜在线观看| 无码人妻热线精品视频| 美女一级免费毛片| 国产免费久久精品99re不卡 | 日韩a在线观看免费观看| 国产精品网址你懂的| 国产精品第一区| 亚洲欧美日韩另类在线一| 国产剧情国内精品原创| 天天视频在线91频| 全裸无码专区| 一级毛片无毒不卡直接观看| 欧美一级高清片欧美国产欧美| 国产91精品调教在线播放| 无码AV高清毛片中国一级毛片| 超碰色了色| 99久久精品免费观看国产| 日a本亚洲中文在线观看| 囯产av无码片毛片一级| 亚洲午夜天堂| 一级毛片在线播放免费观看| 色综合久久无码网| 欧美日韩一区二区三区在线视频| 国产97视频在线观看| 国产精品视频猛进猛出| 国产超碰在线观看| 91在线激情在线观看| 无码日韩视频| 欧美在线一二区| 操国产美女| 色久综合在线| 中文无码影院| 亚洲网综合| 国产精品伦视频观看免费| 无码网站免费观看| 99精品伊人久久久大香线蕉| 黄网站欧美内射| 又粗又硬又大又爽免费视频播放| 国产SUV精品一区二区6| 久久婷婷色综合老司机| 色悠久久综合| 天天综合亚洲| 日韩欧美综合在线制服| 福利姬国产精品一区在线| AV熟女乱| 亚洲日韩每日更新| 免费在线成人网| 青青草国产一区二区三区| 亚洲天堂网在线观看视频| 亚洲综合18p| 久久青草精品一区二区三区| 无码人妻免费| 人妻精品久久久无码区色视| 一区二区三区高清视频国产女人| 国产午夜看片| 亚洲男人的天堂网| av在线人妻熟妇| 2021国产v亚洲v天堂无码| 亚洲国产日韩在线观看| 天天躁日日躁狠狠躁中文字幕| 色网站在线免费观看| 日本国产一区在线观看| 欧美性猛交一区二区三区| 国产青榴视频| 在线视频精品一区| 国产日本一区二区三区| 亚洲视频免费在线看| 99精品热视频这里只有精品7| 中国精品自拍| 人人看人人鲁狠狠高清| 免费一级大毛片a一观看不卡| 蜜臀av性久久久久蜜臀aⅴ麻豆| www.youjizz.com久久| 亚洲欧美日韩精品专区|