• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    A Unification of the Concepts of the Variational Iteration,Adomian Decomposition and Picard Iteration Methods;and a Local Variational Iteration Method

    2016-12-13 01:26:39XuechuanWangSatyaAtluri

    Xuechuan Wang,Satya N.Atluri

    A Unification of the Concepts of the Variational Iteration,Adomian Decomposition and Picard Iteration Methods;and a Local Variational Iteration Method

    Xuechuan Wang1,Satya N.Atluri2

    This paper compares the variational iteration method(VIM),the Adomian decomposition method(ADM)and the Picard iteration method(PIM)for solving a system of first order nonlinear ordinary differential equations(ODEs).A unification of the concepts underlying these three methods is attempted by considering a very general iterative algorithm for VIM.It is found that all the three methods can be regarded as special cases of using a very general matrix of Lagrange multipliers in the iterative algorithm of VIM.The global variational iteration method is briefly reviewed,and further recast into a Local VIM,which is much more convenient and capable of predicting long term complex dynamic responses of nonlinear systems even if they are chaotic.

    1 Introduction

    The rapidly increasing interest of researchers in nonlinear problems has prompted a variety of analytical,semi-analytical and numerical methods.Among them,the variational iteration method(VIM)[Inokuti,Sekine and Mura(1978);He(1999)],the Adomian decomposition method(ADM)[Adomian(1988)]and the Picard iteration method(PIM)[Fukushimam(1997)](often combined with other approximation techniques,such as the modified Chebyshev-Picard iteration(MCPI)method[Woollands,Younes and Junkins(2015)])have received wide attention.Unlike in classical perturbation methods,these three methods do not depend on the existence of a “small parameter”in the nonlinear problem,and they are not limited to cases of weak nonlinearities.The applications of the three methods are relatively straightforward in that there is no need to determine the so-called small parameter or make linearization transformations to the original equation.Further,the results obtained with these methods are valid for the whole parameter domain.Thus they are regularly utilized in areas where the traditional techniques are not applicable.In literature,some comparisons have been made between them on the efficiency and scope of application.The advantages and disadvantages of each method are illustrated in papers[He(1999);Rach(1987)].

    The VIM is capable of solving a large class of nonlinear problems,including both the systems expressed as either ordinary or partial differential equations.It is an analytical asymptotic approach[He(2006)],where the initial guess function is corrected step by step,and finally reaches the true solution.The correctional iterative formula of VIM is similar to that of the Newton-Raphson method,except that the former method is for functions while the latter is for fixed point of nonlinear algebraic equations(NAEs).From the view of Inokuti,Sekine and Mura(1978),we can see that the VIM is somehow an extension of the Newton method to the problem in function space.And the extension can also be made to problems expressed in algebraic,differential,integral,or finite-difference equations,or the combination of them.This idea is illustrated in detail and referred as the general use of Lagrange multipliers in Inokuti,Sekine and Mura(1978).However,for the case of general nonlinear terms,the derivation of exact generalized Lagrange multiplier is normally impossible.Considering that,the VIM of He[He(1999)]omits the nonlinearity in the adjoint equation satisfied by the generalized Lagrange multiplier.Other than that,the VIM of He[He(1999)]incorporates an artificial parameter in the approximation of the periodic solution.By eliminating the secular terms including the parameter,the accuracy of the result can be improved significantly.But this trick does not always work.As is illustrated in the section 3,it seems that only in some special cases can the approach show its efficiency.

    The ADM is a decomposition method.It treats the solution as the sum of a series of functions,and decomposes the nonlinear terms into the rearranged Taylor series expansion in terms of Adomian polynomials.With an initial guess given,the solution is corrected in each iteration step by adding an integral of an Adomian polynomial.As the iteration proceeds,the nonlinear equation is solved gradually.This method provides great convenience for the calculation of the correctional function,but the construction of Adomian polynomials is somewhat complex,especially for the higher order terms.The conditions for convergence of this method have been discussed in literature.It should be noted that the method is just locally convergent because it is based on a Taylor series expansion.Some modifications were made to it by recursively applying it over divided time segments,so that the global convergence can be achieved[Ghosh,Roy and Roy(2007)].

    Compared to VIM and ADM,Picard’s method takes a very simple form.How-ever,its application is very limited,since it needs to integrate the nonlinear terms in each iteration step,which could be very difficult to implement.Compared to Adomian’s method,the Picard iteration lacks the ease of computation and also the ability to solve a wide class of equations.But by combining with other computational techniques,it is still possible to make advances to PIM.For example,the MCPI method[Woolland,Younes and Junkins(2015)]combines the Chebyshev polynomials with the Picard’s method and is applied to the two body gravitational integration problem.It is shown that the integration process of Picard’s method becomes very simple and the method achieves high accuracy and efficiency.

    Despite all the detailed differences in VIM,ADMand PIM,it is the objective of this paper to show that they are all not unrelated.As will be illustrated in section 2,the iteration formula of the Picard’s method and the ADM can be completely deduced from the iterative formula of the VIM.This indicates that a common mathematical principle,which is the general use of Lagrange multipliers,underlies these methods,and a unification of the concepts underlying the three methods can be elucidated.It will help to explain the similarity of the three methods,and how we can combine the ideas in VIM,ADM and PIM to achieve possibly better analytical asymptotic methods.

    After that,a local variational iteration method(LVIM)is proposed in section 3.Firstly,the differential transform(DT)method[Jang,Chen and Liy(2000)]is introduced to help identify the generalized Lagrange multipliers.Then a piecewise solution is obtained by dividing the whole time domain into a finite number of small segments.The forcing term is approximated by a truncated Taylor series expansion in the LVIM,so that the calculation of integrals can be simplified.The advantages of the LVIM are as follows,

    (1)The initial guess function can take a very simple form.Unlike the conventional VIM,there is no need to pay much attention to the construction of an initial guess.A linear function could be good enough for the LVIM,as long as it meets the boundary condition of the problem.

    (2)Instead of restricting the nonlinear term from variation artificially,the nonlinear term can be included in the derivation of generalized Lagrange multipliers in the LVIM.So the correctional formula is more effective.

    (3)The introduction of DT method provides the approximation of generalized Lagrange multipliers in Taylor series expansion.It makes the derivation of generalized Lagrange multipliers much simpler,and also provides convenience to the integration of the correctional formula.

    (4)The LVIM has applicability to long term predictions of many kinds of complex nonlinear problems,including chaotic motions.

    (5)The accuracy and computational speed of LVIM is relatively very high,even compared to the numerical methods such as RK4.

    To convey the ideas in a simple and clear way,a system of first order differential equations is taken for illustration in this paper.However,the conclusions made in the following sections are also applicable to the higher order differential equations.It is well known that a higher order differential equation can always be transformed into a system of first order differential equations.For example,the equation

    can be rewritten as a system of equations like

    by introducing the variables x0,...xn?1.Therefore,the nonlinear ordinary differential equations can be expressed in a general form as

    whereWhere,here onwards,a bold symbol indicates a vector or a matrix,and for brevity,the differential operatoris denoted asLin the following sections.

    In the preceding equation,F(x,τ)is a nonlinear function of the state vectorxand the independent variable τ.

    2 PIM,ADM,and VIM

    The three methods,i.e.VIM,ADM and PIM,were proposed independently by different researchers at different times.ADM and PIM have been well known for a long time,while VIM has a history of no more than four decades.Although they appear to be unrelated to one another at first,a close examination leads to a unification of the three methods with the perspective provided by the concepts of VIM.In the following,the PIM and the ADM are briefly introduced.Then the VIM is described and some transformations are made to elucidate the relationships between all the three methods.

    2.1 PIM

    Consider an initial value problem(IVP)governed by a system of first-order differential equations

    It is equal to the associated integral equations

    The PIM solves this problem in a recursive way by constructing a series of approximating functions.The process works as

    (1)Give an initial guess of the solutionx0(τ)that satisfies the initial conditionx0(t0)=x(t0).

    (2)Substitute it into the recursive formula that holds

    for n≥0.

    2.2 ADM

    For convenience,we rewrite the dynamical system in the preceding IVP asLx=F(x,τ),whereLstands for the first-order derivative.LetL?1be the inverse operator ofL,namely the definite integrator.The equation can be recast as

    withx(t0)being the initial condition.

    The innovative part of ADM is that it approximates both the solution and the nonlinear part as sequences of functions and introduces the Adomian polynomials so that the original problem will be solved progressively.

    This enables us to rewritte the Eq.(3)as

    The solution can thus be derived recursively in the following way

    The Adomian polynomialsAnare generated by simply rearranging the Taylor series expansion ofF(x,τ)aboutˉx0=x(t0).They are expressed as

    2.3 VIM

    Consider the following general nonlinear system

    whereLis the first order differential operator andFis a nonlinear operator.The solution of this system can be approximated with an initial approximationx0(t)and the correctional formula as

    where λ (τ)is a matrix of Lagrange multipliers which are yet to be determined.Suppose Π [x(τ),λ (τ)]is a vector function ofx(τ)and λ (τ),where t0≤τ≤t.

    Let?x(τ)be the exact solution ofLx(τ)=F[x(τ),τ].Naturally,it satisfies the expression

    Now we want to make the function Π [x(t),λ (t)]stationary aboutxatx(t)=?x(t).Firstly,the variation of Π [x(t),λ (t)]is derived as

    IfFis not an explicit function of τ,the term[?F(x,τ)/?τ]δτ can be omitted in the preceding formula.

    Then we collect the terms including δx(τ)|τ=tand δx(τ),

    Note that the boundary value ofx(τ)at τ=t0is prescribed,that is to say δx(τ)|τ=t0=0.Thus the stationary condition for Π [x(τ),λ (τ)]is obtained as

    Noting that the exact solution?xis unknown,therefore the truly optimal λ (τ)is not available herein.As an alternative,λ (τ)is approximated by replacing?xwithxn.Ifxnis a neighbored function of?x,i.e.?x?xn=δ?x,the error caused will not exceed O2(δ?x).

    In the VIM of He(1999),sometimes the nonlinear term is considered as restricted from variation,in order to simplify the calculation of λ .As an asymptotic method for solving nonlinear problems,VIM provides plenty of freedom to the user.It does not restrict the selection of initial guess functions and even allows the existence of unknown parameters in it.Further,depending on the selection of a restricted variation,the matrix of generalized Lagrange multipliers λ can be determined in various ways.The less restricted the variational terms are,the more accurate will the Lagrange multipliers be.

    3 VIM and its relationship to PIM and ADM

    Instead of PIM,if VIM of He[He(1999)]is utilized to solve the IVP governed by Eq.(1),we will get the corresponding correctional formula with λ

    where?F(xn,τ)is considered as being restricted from variation.The matrix of Lagrange multipliers λ are determined from the following restricted stationary conditions:

    The Lagrange multipliers,therefore,can be identified as λ (τ)=diag[?1,?1,...].As a result,we obtain the following iteration formula

    Noting thatxn(t0)is equal to the initial valuex(t0),it is exactly the Picard iteration.In the paper of He(1999),a first-order differential equation is considered and the ADM is proved to be a specific version of VIM.This work is reviewed in the following.

    Consider the equation Lx+Rx+Nx=g(τ),in which L is a first-order differential operator,while R is a linear operator.For this equation,the correctional formula of VIM is

    Considering Rxn+Nxnas restricted from variation,λ can be easily obtained as?1 in this case.Suppose x0=ˉx0=x(t0)+L?1g(τ),substituting it into the formula gives

    Similarly the iteration process gives

    Ifˉx1is regarded as relatively small,expanding the nonlinear functional aboutˉx0and ignoring the higher-level small functions leads to N(ˉx0+ˉx1)=N(ˉx0)+ˉx1N′(ˉx0).Then we have

    Further,in the nth step of iteration,regardˉxn?1as the relatively small function and omit the terms O(ˉxm1ˉxm2...ˉxmk),where m1+m2+...+mk≥n.It leads to the correctional formula of ADM.

    It is worth to note that the above comparison can also be made in the case of multidimensional system.Herein,a general form of the nonlinear ordinary differential equations is taken for illustration.

    Lx=F(x,τ),

    The ADM gives an iteration formula as

    while the correctional formula of VIM is

    Set λ as diag[?1,?1,...]and follow the manipulations we made to the one dimensional case.Supposex0=ˉx0,from the correctional formula of VIM we have

    Accordingly,

    Regardˉx1as a relatively small amplitude vector of functions.The nonlinear term can be rewritten as

    Omitting the small terms O(ˉx21),it leads to

    As the iteration goes on,we have

    which becomes the ADM correctional formula by omitting the small terms

    From the statements above,it is clear that the ADM and PIM can both be regarded as variants of VIM.Although they are treated as different methods in literature,and were developed into various forms such as modified ADM and MCPI method,a common principle,namely the general use of Lagrange multipliers guides all of them.

    4 Local variational iteration method

    4.1 Limitations of Global VIM

    It has been shown that the VIM is effective in approximating the periodic motions of nonlinear systems,especially the limit cycle oscillation.To verify that,an unforced Duffing equation is solved with VIM herein.The governing equation is

    where c=0,k1=1,k2=1.

    In the approach proposed by He(1999),the Lagrange multiplier λ is derived with the nonlinear term being restricted from variation,and it is shown to be tenable in the unforced case with the initial conditions x1(t0)=A,x2(t0)=0 and initial guess x1(t)=Acosωt.The correctional formula is

    With the nonlinear term?x31restricted from variation,the matrix of Lagrange multipliers satisfies the condition:

    Therefore,the matrix of Lagrange multipliers λ (τ)is identified as

    The solution corrected only once with VIM can already give a result comparable with that of the RK4 method.The example of an unforced Duffing equation shows that the VIM could be very efficient in approximating the periodic solution of nonlinear systems by eliminating the secular term appearing in the correctional formula.The numerical results are plotted in Fig.1 with the initial condition x(0)=1,˙x(0)=0 and initial guess x1(t)=cosωt.By substituting the initial guess into the correctional formula of VIM,we can obtain a corrected solution,of which the secular term is(2?17ω2+9ω4)cos(t).To eliminate the secular term,the frequency is determined to be ω=1.32774.

    However,in a more general case withx(t0)=?A B?T,the performance of VIM is far inferior to the preceding one.The initial guess is selected as x(t)=Acosωt+(B/ω)sinωt.Let A=1 and B=1,the variational iteration formula gives the corrected solution,in which the terms of sint and cost are

    It is easily identified that the coefficients of sint and cost can not both be equal to zero.Here only the term of cost is eliminated,leading to the approximation ω≈1.4493.Substituting it into the corrected solution,we have x1(t)=1.00598cos(1.4493t)?0.00597962cos(4.34791t)?0.111705sin(t)+0.694112sin(1.4493t)+0.0243163sin(4.34791t).

    Figure 1:Comparison of VIM and RK4.solid line:RK4.Dash line:VIM

    This result is compared with that of RK4 in Fig.2.

    It is believed that if both terms of sint and cost vanish,the corrected solution will be more accurate.But with just one parameter ω in the initial guess,this will be an impossible mission.Thus we tried to bring in two parameters ω1, ω2in the initial guess,which takes the form x1(t)=Acosω1t+(B/ω2)sinω2t.If the resulted corrected solution can be obtained analytically,the accuracy will be much improved.But along with it,the calculation of variational iteration formula will become much more complex.We found that even with the help of Mathematica,this job is still too burdensome to be accomplished.

    Moreover,it is also indicated that a carefully selected initial guess is important to ensure the convergence of this method.If the approximate function takes the form as x1(t)=A+Bt rather than x1(t)=Acosωt,the result given by VIM will be more and more divergent from the true solution with the iteration.It is illustrated in the case A=1 and B=1.With VIM,the initial guess x1,0(t)=1+t is corrected for 2 times and the results are given as follows.

    The solution x1,1(t)is compared with that of RK4 in Fig.3.

    Figure 2:Solid line:RK4.Dash line:VIM

    Figure 3:Solid line:RK4.Dotted line:VIM

    The above numerical examples and analysis show that VIM is fastidious to the initial condition of the problem,and also too sensitive to the initial guess function,which is analogous to the Newton’s iteration method for solving nonlinear algebraic equations(NAEs).This method also fails to predict the long term responses of some complex nonlinear phenomenon such as chaos or quasi periodic motion,because these patterns of motion are non-periodic and can hardly be approximated by any analytical functions.

    4.2 Methodology of LVIM

    As is shown in Atluri (2005)and Dong,Alotaibi,Mohiuddine and Atluri(2014),local approximations are far more effective in solving initial value and boundary value problems.We use these ideas in developing a local variational iteration method.To remedy the drawbacks shown above,we propose to further improve the VIM.Here is the primary idea.Firstly,the nonlinear term should be kept in the derivation of Lagrange multipliers so that a better correction can be obtained for arbitrary approximation.And then the entire time domain is divided into small intervals and the LVIM is applied repeatedly in each time interval ti?1≤ τ≤ti.In each interval ti?1≤ τ≤ ti,we approximate locally using an arbitrary functionx(τ)=A+Bτ for instance;then the corrected solution will be a function ofAandB,i.e.the initial condition of each interval,thus the solution in the entire time domain can be obtained by the repetition ofAandBstep by step.

    With this in mind,the IVP of the unforced Duffing equation is solved with the initial condition as x1(t0)=1,x2(t0)=1.The correctional formula is

    where ti?1≤ t≤ ti.The stationary condition is

    Instead of neglecting the nonlinear term,we kept it in the derivation of λ .This makes the differential equations more difficult to be solved analytically.But an exact analytical solution is not necessary here.The differential transform(DT)method is introduced to obtain the approximated λ in the form of power series.With the nonlinear term included,the approximated λ is actually more precise than that in the VIM of He(1999).

    Supposing that the initial guess of x1takes the form as x1,0(τ)=Bτ+A,an approximated λ can be obtained using the DT method.Herein only the first four DT terms of λ (τ)is kept.

    where ti?1≤ t≤tiand ti?1≤ τ≤ti.Using the correctional formula for one time,we have

    where ti?1≤ t≤ ti.

    In Fig.4(a),the preceding solution is compared with the numerical results obtained by ODE45,for which the relative and absolute accuracy are both set as 10?15.Fig.4(b)shows the computational error of the LVIM with respect to ODE45.The time step size is set as Δt=0.01 and the simulation is carried in t=[0,100].From Fig.4(b),it can be seen that the error of the LVIM is less than 10?3even for t=100.

    Figure 4:Comparison of the LVIM and the ODE45.

    To further inspect the performance of the LVIM in predicting complicated responses,the forced Duffing equation is investigated.

    It is shown to be a chaotic system by RK4 and other reliable methods.Unsurprisingly,the VIM of He[He(1999)]fails to solve it because the solution is too complicated to be computable.An approximation is made to the forced term of the system.Since it is involved in the integration,a simple form of it will accelerate the calculation.For that,the forced term is expanded into power series as well.The stationary condition for λ is

    where ti?1≤ τ≤tiand ti?1≤t≤ti.

    Suppose the initial guess in each interval ti?1≤ τ≤ tiis x0(τ)=A+Bτ.With DT method,we can get the fourth-order approximated λ(τ,t)in terms of A and B,which is

    Then substituting it into the correction formula,with the forced term being replaced by power series,will lead to the corrected function x1(t)

    where ti?1≤ t≤ ti.It has a compacted form with A,B,ti?1(denoted as C in the corrected function)as the coefficients.When applied to multiple time intervals,it only needs to recompute A,Band t0,which correspond to the final position,velocity and time instant of the last interval.The LVIM is very fast,precise and efficient on solving nonlinear problems.Fig.5 shows the chaotic motions of a forced Duffing oscillator predicted by both LVIM and RK4 methods.The time step size of LVIM is set as 0.2,while that of RK4 needs to be set as 0.02 in order to achieve the same accuracy.The numerical simulation carried out with MATLB shows the LVIM is almost 10 times faster than RK4 in this case.The computation time of the LVIM and the RK4 are 0.013471s and 0.135445s,respectively.

    Similar to the modified ADM,the variational iteration method is further improved herein by bringing in some local approximation techniques,so as to be practical in predicting long term motion and complex dynamical responses.The LVIM makes the derivation of Lagrange multipliers λ and the calculation of the variational iteration formula much simpler,yet still provide reliable solution in each sub time domain.

    Figure 5:Comparison of the LVIM and the RK4 method for the forced Duffing oscillator.Dots:LVIM.Solid line:RK4

    5 Conclusions

    The Picard iteration method,the Adomian decomposition method and the variational iteration method are unified in this paper.From the view point provided by the concept of generalized Lagrange multipliers,it is found that the three methods are rooted in the same mathematical approach of generalized Lagrange multipliers.By some simple transformations,the Picard iteration formula and Adomian decomposition formula can be derived from the variational iteration formula.

    However,the Global VIM is very laborious in practice and is unable to predict long term responses of nonlinear systems,although it is able to provide an analytical solution.Being aware of that,we introduced a Local VIM that simplified the derivation of the matrix of generalized Lagrange multipliers λ with the DT method.By dividing the entire domain into small intervals and updating the initial conditions in each sub-interval,the Local VIM method achieves very high accuracy and very high computational speeds.

    Acknowledgement:The support of Texas Tech University,as well as the support of the Chinese Scholarship Council to the first author are gratefully acknowledged.

    Adomian,G.(1988):A review of the decomposition method in applied mathematics.Journal of Mathematical Analysis and Applications,vol.135,pp.501–544.

    Atluri,S.N.(2005):Methods of Computer Modeling in Engineering&The Sciences,Volume I.Tech Science Press,Forsyth.

    Dong,L.;Alotaibi,A.;Mohiuddine,S.A.;Atluri,S.N.(2014):Computational methods in engineering:a variety of primal&mixed methods,with global&local interpolations,for well-posed or Ill-posed BCs.Computer Modeling in Engineering&Sciences,vol.99,no.1,pp.1–85.

    Fukushimam,T.(1997):Picard iteration method,Chebyshev polynomial approximation,and global numerical integration of dynamical motions.The Astronomical Journal,vol.113,no.5,pp.1909–1914.

    Ghosh,S.;Roy,A.;Roy,D.(2007):An adaption of Adomian decomposition for numerical-analytic integration of strongly nonlinear and chaotic oscillators.Computational Methods in Applied Mechanics and Engineering,vol.196,pp.1133–1153.

    He,J.(1999):Variational iteration method:a kind of non-linear analytical technique:some examples.International Journal of Non-linear Mechanics,vol.34,pp.699–708.

    He,J.(2006):Some asymptotic methods for strongly nonlinear equations.International Journal of Modern Physics B,vol.20,no.10,pp.1141–1199.

    Inokuti,M.;Sekine,H.;Mura,T.(1978):General use of the Lagrange multiplier in nonlinear mathematical physics,in:S.Nemat-Nasser(Ed.),Variational Method in the Mechanics of Solids.Pergamon Press,Oxford,pp.156–162.

    Jang,M.;Chen,C.;Liy,Y.(2000):On solving the initial-value problems using the differential transformation method.Applied Mathematics and Computation,vol.115,pp.145–160.

    Rach,R.(1987):On the Adomian decomposition method and comparisons with Picard iteration.Journal of Mathematical Analysis and Applications,vol.128,pp.480–483.

    Woollands,R.M.;Younes,A.B.;Junkins,J.L.(2015):New solutions for the perturbed lambert problem using regularization and picard iteration.Journal of Guidance,Control,and Dynamics,vol.38,no.9,pp.1548–1562.

    1Texas Tech University;Center for Advanced Research in the Engineering Sciences;Visiting Scholar from Northwestern Polytechnical University,China.

    2Texas Tech University;Center for Advanced Research in the Engineering Sciences.

    亚洲少妇的诱惑av| 亚洲精品一区蜜桃| 欧美日韩av久久| 免费一级毛片在线播放高清视频 | 在线观看国产h片| 精品一区二区三区四区五区乱码 | 色婷婷久久久亚洲欧美| 水蜜桃什么品种好| 国产成人免费无遮挡视频| www.熟女人妻精品国产| 国产一区二区 视频在线| 高清不卡的av网站| 精品卡一卡二卡四卡免费| xxxhd国产人妻xxx| 日本av手机在线免费观看| 黑人欧美特级aaaaaa片| 久久 成人 亚洲| 亚洲男人天堂网一区| 久久亚洲国产成人精品v| www.999成人在线观看| 十八禁网站网址无遮挡| 啦啦啦 在线观看视频| 欧美黑人欧美精品刺激| 亚洲美女黄色视频免费看| 91麻豆av在线| 日本av手机在线免费观看| 一边摸一边做爽爽视频免费| 亚洲精品一二三| 香蕉丝袜av| 国产一区二区在线观看av| 啦啦啦中文免费视频观看日本| 日韩免费高清中文字幕av| 超碰成人久久| 成年av动漫网址| 1024香蕉在线观看| 精品人妻熟女毛片av久久网站| 免费在线观看完整版高清| 亚洲中文av在线| 香蕉国产在线看| 制服人妻中文乱码| 亚洲av综合色区一区| 欧美老熟妇乱子伦牲交| 久久女婷五月综合色啪小说| 欧美激情高清一区二区三区| 超色免费av| 婷婷色综合大香蕉| 亚洲 欧美一区二区三区| 下体分泌物呈黄色| 亚洲欧美激情在线| 欧美日韩黄片免| 久久国产精品男人的天堂亚洲| 看免费成人av毛片| 老司机影院毛片| 黄频高清免费视频| 美女大奶头黄色视频| 超色免费av| 国产成人一区二区三区免费视频网站 | 成人影院久久| 晚上一个人看的免费电影| 欧美黄色片欧美黄色片| av福利片在线| 日日爽夜夜爽网站| 丁香六月天网| 久久亚洲国产成人精品v| 国产不卡av网站在线观看| 国产亚洲欧美精品永久| 男女国产视频网站| 91精品伊人久久大香线蕉| 亚洲成人免费电影在线观看 | 人人妻人人爽人人添夜夜欢视频| 免费看av在线观看网站| 免费一级毛片在线播放高清视频 | 午夜激情av网站| avwww免费| 天天操日日干夜夜撸| av又黄又爽大尺度在线免费看| av电影中文网址| 女人高潮潮喷娇喘18禁视频| 男女下面插进去视频免费观看| 中文精品一卡2卡3卡4更新| av不卡在线播放| 欧美av亚洲av综合av国产av| 久久精品成人免费网站| 久久青草综合色| 午夜视频精品福利| 国产一区二区激情短视频 | 99国产精品99久久久久| 国产成人系列免费观看| 欧美激情高清一区二区三区| 伦理电影免费视频| 女人被躁到高潮嗷嗷叫费观| 精品久久久久久久毛片微露脸 | 夫妻午夜视频| 久久精品久久久久久久性| 亚洲成av片中文字幕在线观看| 青草久久国产| 中文字幕高清在线视频| 亚洲成人免费av在线播放| 国产有黄有色有爽视频| www.精华液| 又紧又爽又黄一区二区| 久久久久精品人妻al黑| 亚洲免费av在线视频| 婷婷色综合www| 国产精品麻豆人妻色哟哟久久| 肉色欧美久久久久久久蜜桃| 激情五月婷婷亚洲| 国产爽快片一区二区三区| 亚洲人成77777在线视频| 免费不卡黄色视频| 亚洲欧美一区二区三区国产| www.av在线官网国产| 人人妻人人澡人人看| 91九色精品人成在线观看| 午夜福利视频在线观看免费| 国产成人精品久久二区二区91| 国产男女内射视频| 久久精品久久久久久久性| 大香蕉久久网| 啦啦啦在线观看免费高清www| 亚洲国产精品成人久久小说| 欧美日本中文国产一区发布| 国产一区二区在线观看av| 夜夜骑夜夜射夜夜干| 少妇精品久久久久久久| 欧美日韩成人在线一区二区| 精品国产国语对白av| 丝袜喷水一区| 超色免费av| xxxhd国产人妻xxx| 涩涩av久久男人的天堂| 亚洲欧美精品自产自拍| 50天的宝宝边吃奶边哭怎么回事| 在线观看人妻少妇| 日韩大码丰满熟妇| 国产免费一区二区三区四区乱码| 久久 成人 亚洲| av线在线观看网站| 日韩人妻精品一区2区三区| 青春草亚洲视频在线观看| 每晚都被弄得嗷嗷叫到高潮| 国产有黄有色有爽视频| 亚洲国产精品一区三区| 丁香六月天网| 91国产中文字幕| 精品少妇内射三级| 80岁老熟妇乱子伦牲交| 少妇人妻 视频| 最近手机中文字幕大全| 深夜精品福利| 首页视频小说图片口味搜索 | 久9热在线精品视频| 午夜日韩欧美国产| 另类亚洲欧美激情| 亚洲国产av新网站| 亚洲专区国产一区二区| 久久精品人人爽人人爽视色| 多毛熟女@视频| 超碰97精品在线观看| 女人久久www免费人成看片| 女人精品久久久久毛片| 丰满饥渴人妻一区二区三| 五月开心婷婷网| 18禁裸乳无遮挡动漫免费视频| 国产精品99久久99久久久不卡| 欧美激情 高清一区二区三区| 久久精品久久久久久噜噜老黄| 我要看黄色一级片免费的| 80岁老熟妇乱子伦牲交| 少妇人妻久久综合中文| 国产视频一区二区在线看| 19禁男女啪啪无遮挡网站| 女人精品久久久久毛片| 天堂俺去俺来也www色官网| 五月开心婷婷网| 精品欧美一区二区三区在线| 美女视频免费永久观看网站| 久久99一区二区三区| 久久午夜综合久久蜜桃| 国产精品麻豆人妻色哟哟久久| 亚洲一区二区三区欧美精品| 久久精品国产综合久久久| 日韩视频在线欧美| 国产成人av教育| 男人爽女人下面视频在线观看| 亚洲国产精品999| 国产精品秋霞免费鲁丝片| 欧美精品啪啪一区二区三区 | 国产男女超爽视频在线观看| 制服人妻中文乱码| 亚洲av在线观看美女高潮| 国产福利在线免费观看视频| 中文字幕av电影在线播放| 国产精品麻豆人妻色哟哟久久| 亚洲人成电影观看| 久久精品国产综合久久久| 久久久精品免费免费高清| 国产熟女午夜一区二区三区| 欧美精品亚洲一区二区| 国产高清不卡午夜福利| 国产视频首页在线观看| videos熟女内射| 亚洲精品自拍成人| 美国免费a级毛片| 亚洲人成网站在线观看播放| 国产精品秋霞免费鲁丝片| 国产精品一区二区精品视频观看| 久久国产精品大桥未久av| 一级毛片电影观看| www.自偷自拍.com| 国产黄频视频在线观看| 国语对白做爰xxxⅹ性视频网站| 国产精品久久久人人做人人爽| 日韩伦理黄色片| 美女国产高潮福利片在线看| 成人亚洲欧美一区二区av| √禁漫天堂资源中文www| 中文字幕亚洲精品专区| 18在线观看网站| 国产免费一区二区三区四区乱码| 国产欧美日韩精品亚洲av| 精品久久蜜臀av无| 亚洲 国产 在线| 美女大奶头黄色视频| 少妇猛男粗大的猛烈进出视频| 日韩精品免费视频一区二区三区| 精品一区二区三区四区五区乱码 | 亚洲欧美一区二区三区黑人| 黄色毛片三级朝国网站| 18禁国产床啪视频网站| 亚洲国产欧美在线一区| 国产深夜福利视频在线观看| 最新的欧美精品一区二区| netflix在线观看网站| 久久天堂一区二区三区四区| 777米奇影视久久| 纯流量卡能插随身wifi吗| cao死你这个sao货| 午夜福利影视在线免费观看| 成年人免费黄色播放视频| 啦啦啦在线观看免费高清www| 久久久久久久大尺度免费视频| 丝袜美腿诱惑在线| 国产av国产精品国产| 日韩 欧美 亚洲 中文字幕| 日本午夜av视频| 手机成人av网站| 大香蕉久久网| 欧美日韩视频高清一区二区三区二| 国产一区有黄有色的免费视频| 久久九九热精品免费| 五月天丁香电影| 国产精品.久久久| 久久久久精品国产欧美久久久 | 亚洲精品久久成人aⅴ小说| 少妇人妻 视频| 一边摸一边做爽爽视频免费| 亚洲av欧美aⅴ国产| 丝袜美足系列| 看免费成人av毛片| 欧美日韩视频高清一区二区三区二| 日韩视频在线欧美| 在线av久久热| 亚洲av成人不卡在线观看播放网 | 久久影院123| 久久久国产欧美日韩av| 女警被强在线播放| 黄色 视频免费看| 99国产精品一区二区蜜桃av | 国产亚洲精品久久久久5区| 啦啦啦中文免费视频观看日本| 久久久精品免费免费高清| 七月丁香在线播放| 成人国产一区最新在线观看 | 中文字幕亚洲精品专区| 国产一区二区三区av在线| 亚洲精品日本国产第一区| 亚洲成色77777| av天堂久久9| 国产真人三级小视频在线观看| 精品国产乱码久久久久久小说| 97精品久久久久久久久久精品| 免费女性裸体啪啪无遮挡网站| xxxhd国产人妻xxx| 国产一卡二卡三卡精品| 激情五月婷婷亚洲| 两个人免费观看高清视频| 一本—道久久a久久精品蜜桃钙片| 欧美亚洲日本最大视频资源| 成人手机av| a级毛片黄视频| 最近中文字幕2019免费版| 日本av免费视频播放| 水蜜桃什么品种好| 交换朋友夫妻互换小说| 国产精品久久久久成人av| 国产精品九九99| 9色porny在线观看| 亚洲第一青青草原| 久久ye,这里只有精品| 波多野结衣av一区二区av| 在线观看免费午夜福利视频| 国产一级毛片在线| 国产精品一国产av| 久久久精品区二区三区| 免费看十八禁软件| 亚洲专区国产一区二区| 亚洲av日韩在线播放| 亚洲国产欧美网| 国产欧美日韩精品亚洲av| 中文字幕最新亚洲高清| 久久天堂一区二区三区四区| 国产片内射在线| 国产精品麻豆人妻色哟哟久久| 亚洲国产av新网站| 黄色 视频免费看| 18在线观看网站| 下体分泌物呈黄色| 性色av乱码一区二区三区2| 国产成人免费无遮挡视频| 免费看不卡的av| 中文精品一卡2卡3卡4更新| av在线app专区| av福利片在线| 超碰97精品在线观看| 亚洲一区二区三区欧美精品| 亚洲欧美日韩另类电影网站| 人成视频在线观看免费观看| 在现免费观看毛片| 精品久久久久久久毛片微露脸 | 亚洲欧洲日产国产| 日韩 欧美 亚洲 中文字幕| 在线av久久热| 日韩制服丝袜自拍偷拍| 久久久亚洲精品成人影院| 丝袜美足系列| 夫妻性生交免费视频一级片| 成人影院久久| 涩涩av久久男人的天堂| 视频区图区小说| 欧美日韩黄片免| 十八禁人妻一区二区| 9热在线视频观看99| 免费女性裸体啪啪无遮挡网站| 王馨瑶露胸无遮挡在线观看| 夫妻性生交免费视频一级片| 亚洲国产最新在线播放| 午夜福利在线免费观看网站| 黄色毛片三级朝国网站| 看十八女毛片水多多多| 精品卡一卡二卡四卡免费| 岛国毛片在线播放| 国产1区2区3区精品| 亚洲,欧美,日韩| 欧美另类一区| 久久久欧美国产精品| 欧美中文综合在线视频| 两个人免费观看高清视频| 美女扒开内裤让男人捅视频| 人妻 亚洲 视频| 日本猛色少妇xxxxx猛交久久| 巨乳人妻的诱惑在线观看| 久久久久久亚洲精品国产蜜桃av| 涩涩av久久男人的天堂| 成年女人毛片免费观看观看9 | 婷婷色综合大香蕉| 国产爽快片一区二区三区| 性高湖久久久久久久久免费观看| 在线观看免费视频网站a站| 少妇猛男粗大的猛烈进出视频| 两个人看的免费小视频| 国产精品国产三级专区第一集| 在线观看免费日韩欧美大片| 日韩大码丰满熟妇| 美女脱内裤让男人舔精品视频| 99国产精品一区二区蜜桃av | 狂野欧美激情性bbbbbb| 亚洲av成人不卡在线观看播放网 | 超碰成人久久| 日本欧美国产在线视频| 飞空精品影院首页| 在线精品无人区一区二区三| 亚洲国产精品一区二区三区在线| 男人操女人黄网站| 午夜激情久久久久久久| av国产精品久久久久影院| 国产精品一区二区精品视频观看| 精品人妻1区二区| 色94色欧美一区二区| 麻豆国产av国片精品| 亚洲国产精品一区三区| 又紧又爽又黄一区二区| 国产真人三级小视频在线观看| 天天添夜夜摸| 天天躁夜夜躁狠狠久久av| 国产精品秋霞免费鲁丝片| 亚洲精品久久午夜乱码| 亚洲欧美精品综合一区二区三区| 亚洲伊人色综图| 亚洲av日韩精品久久久久久密 | 99精品久久久久人妻精品| 国产高清videossex| 精品福利观看| 国产淫语在线视频| 成年美女黄网站色视频大全免费| 亚洲伊人久久精品综合| 黄色怎么调成土黄色| 亚洲专区国产一区二区| 两个人免费观看高清视频| 亚洲成色77777| 国产麻豆69| 丁香六月欧美| av在线app专区| 精品免费久久久久久久清纯 | 国产精品一区二区在线观看99| 少妇人妻 视频| 手机成人av网站| 国产欧美亚洲国产| 日本欧美视频一区| 亚洲一区二区三区欧美精品| 日日摸夜夜添夜夜爱| 国产精品久久久久久精品古装| 黄片播放在线免费| 麻豆乱淫一区二区| 一区福利在线观看| 国产女主播在线喷水免费视频网站| 两性夫妻黄色片| 啦啦啦在线观看免费高清www| 在现免费观看毛片| 亚洲国产毛片av蜜桃av| 免费女性裸体啪啪无遮挡网站| 国产麻豆69| 国产老妇伦熟女老妇高清| 亚洲第一av免费看| 欧美精品人与动牲交sv欧美| 亚洲成国产人片在线观看| 99久久综合免费| 国产亚洲av片在线观看秒播厂| 成人亚洲欧美一区二区av| 日韩一卡2卡3卡4卡2021年| 9191精品国产免费久久| av天堂在线播放| 欧美日本中文国产一区发布| 如日韩欧美国产精品一区二区三区| 50天的宝宝边吃奶边哭怎么回事| 亚洲av综合色区一区| 99re6热这里在线精品视频| 老司机午夜十八禁免费视频| 午夜福利影视在线免费观看| 无遮挡黄片免费观看| 中文字幕另类日韩欧美亚洲嫩草| 爱豆传媒免费全集在线观看| 一本久久精品| 十八禁高潮呻吟视频| 中文字幕另类日韩欧美亚洲嫩草| 操美女的视频在线观看| 久久精品国产综合久久久| 欧美+亚洲+日韩+国产| 亚洲五月婷婷丁香| 飞空精品影院首页| 黄色怎么调成土黄色| 亚洲欧美日韩高清在线视频 | 黄片小视频在线播放| 国产福利在线免费观看视频| 黄色a级毛片大全视频| 999精品在线视频| 美女高潮到喷水免费观看| 最新在线观看一区二区三区 | 亚洲欧美一区二区三区久久| 欧美日韩黄片免| 精品国产一区二区三区久久久樱花| 国产精品一国产av| 人妻人人澡人人爽人人| 青春草视频在线免费观看| 日本午夜av视频| 亚洲精品第二区| 波野结衣二区三区在线| 国产高清videossex| 亚洲成人免费av在线播放| av网站免费在线观看视频| 国产一卡二卡三卡精品| 最黄视频免费看| 老司机亚洲免费影院| 久久av网站| 亚洲精品乱久久久久久| 精品人妻1区二区| 日本色播在线视频| 热99国产精品久久久久久7| 一本一本久久a久久精品综合妖精| 亚洲综合色网址| 精品亚洲成a人片在线观看| 亚洲国产欧美一区二区综合| 亚洲色图综合在线观看| 亚洲国产看品久久| 男人操女人黄网站| 天天躁日日躁夜夜躁夜夜| 一区福利在线观看| 国产精品免费大片| 国产野战对白在线观看| 亚洲伊人久久精品综合| 午夜免费成人在线视频| 亚洲天堂av无毛| 一级a爱视频在线免费观看| 99精国产麻豆久久婷婷| 人人澡人人妻人| av欧美777| 18禁裸乳无遮挡动漫免费视频| 国产欧美日韩一区二区三区在线| 亚洲一区二区三区欧美精品| 最近手机中文字幕大全| 久久鲁丝午夜福利片| 精品一区在线观看国产| 日韩精品免费视频一区二区三区| 七月丁香在线播放| 晚上一个人看的免费电影| 欧美亚洲 丝袜 人妻 在线| 777久久人妻少妇嫩草av网站| 精品国产一区二区三区久久久樱花| 精品少妇久久久久久888优播| 午夜影院在线不卡| 青草久久国产| 在线 av 中文字幕| 19禁男女啪啪无遮挡网站| 国产激情久久老熟女| www.自偷自拍.com| 18在线观看网站| 欧美人与性动交α欧美软件| www.av在线官网国产| 日韩av免费高清视频| 999久久久国产精品视频| 免费高清在线观看视频在线观看| 欧美乱码精品一区二区三区| 久久av网站| av网站在线播放免费| www.自偷自拍.com| 欧美成狂野欧美在线观看| 国产老妇伦熟女老妇高清| 夫妻性生交免费视频一级片| 两个人看的免费小视频| tube8黄色片| 少妇人妻久久综合中文| 精品国产超薄肉色丝袜足j| 久久99一区二区三区| 国产成人a∨麻豆精品| 一本一本久久a久久精品综合妖精| 亚洲av男天堂| 黑人巨大精品欧美一区二区蜜桃| 后天国语完整版免费观看| 欧美成人午夜精品| 电影成人av| 午夜福利影视在线免费观看| 欧美激情极品国产一区二区三区| 黑人猛操日本美女一级片| 国产成人系列免费观看| 国产成人免费观看mmmm| 久久久久久人人人人人| 久久性视频一级片| 国产精品一区二区精品视频观看| 久久人妻熟女aⅴ| 欧美人与性动交α欧美软件| 久久精品国产亚洲av高清一级| 亚洲国产欧美网| 亚洲精品美女久久久久99蜜臀 | 午夜两性在线视频| 国产女主播在线喷水免费视频网站| 国产免费一区二区三区四区乱码| 中文字幕人妻丝袜一区二区| 国产人伦9x9x在线观看| 日本欧美视频一区| 操美女的视频在线观看| 乱人伦中国视频| 亚洲欧洲国产日韩| 99久久人妻综合| 人人妻人人澡人人看| 亚洲精品国产av蜜桃| 别揉我奶头~嗯~啊~动态视频 | 高清欧美精品videossex| 成人免费观看视频高清| 黑人欧美特级aaaaaa片| h视频一区二区三区| 亚洲国产欧美在线一区| 精品人妻一区二区三区麻豆| 久久精品久久久久久噜噜老黄| 国产一级毛片在线| 麻豆国产av国片精品| 韩国高清视频一区二区三区| 成人手机av| 久久精品久久久久久久性| 嫁个100分男人电影在线观看 | 99国产精品一区二区蜜桃av | 亚洲欧美清纯卡通| 夫妻性生交免费视频一级片| 91精品伊人久久大香线蕉| 国产精品九九99| 男男h啪啪无遮挡| 99精国产麻豆久久婷婷| 国产精品久久久久成人av| 中文字幕制服av| 日本色播在线视频| 好男人电影高清在线观看| 波多野结衣一区麻豆| 少妇粗大呻吟视频| 亚洲一卡2卡3卡4卡5卡精品中文| 美女视频免费永久观看网站| 后天国语完整版免费观看| 中文字幕亚洲精品专区| 国产欧美日韩综合在线一区二区| 悠悠久久av| 日韩一卡2卡3卡4卡2021年| 国产亚洲精品第一综合不卡| 爱豆传媒免费全集在线观看| 亚洲欧美色中文字幕在线| 精品久久蜜臀av无| 巨乳人妻的诱惑在线观看| 可以免费在线观看a视频的电影网站|