计算物理ComputationalPhysics计算物理 (5).pdf

Computational physiCsChaosOrdinary differential equations Initial-value problems The Euler methods Predictor-corrector methods The Runge-Kutta method Chaotic dynamics Boundary-value problems The shooting method Linear equations Eigenvalue problemsthe motion of a classical particle is described by Newtons equationThe motion of a quantum particle is described by the Schrodinger equationThe dynamics and statics of bulk materials such as fluids and solids are all described by differential equations.Most problems in physics and engineering appear in the form of differential equations.22dtrdmdtvdmamfVmti222For example In general,we can classify ordinary differential equations into three major categories:initial-value problemstime-dependent equations with given initial conditionsboundary-value problemsdifferential equations with specified boundary conditionseigenvalue problemssolutions for selected parameters(eigenvalues)in the equationsInitial-value problems Typically,initial-value problems involve dynamical systems.For example,the motion of the moon,earth,and sun,the dynamics of a rocket,or the propagation of ocean waves.A dynamical system can be described by a set of first-order differential equations:),(),(),(),(21tygtygtygtyglthe generalized velocity vectorthe generalized position vector),(tygdtyd),(21lyyyyExample A particle moving in one dimension under an elastic force Define y1=x;y2=v;Then we obtain:If the initial position y1(0)=x(0)and the initial velocity y2(0)=v(0)are given,we can solve the problem numerically.xkdtvdmamf,1221ymkdtdyydtdyThe Euler methodThe accuracy of this algorithm is relatively low.At the end of the calculation after a total of n steps,the error accumulated in the calculation is on the order of nO(t2)O(t).iiiiiiiiiiittOgyytygttyydtdy12111)(),(We can formally rewrite the above equation as an integralwhich is the exact solution if the integral can be obtained exactly.dttygyyjiittiji),(Because we can not obtain the integral exactly in general,we have to approximate it.The accuracy in the approximation of the integral determines the accuracy of the solution.If we take the simplest case of j=1 and approximate g(y,t)=gi in the integral,we recover the Euler algorithm.Code example 4.1.Euler.cpp(1.3.Intro.cpp)Predictor-corrector method Use the solution from the Euler method as the starting point.Use a numerical quadrature to carry out the integration.For example,if we choose j=1 and use the trapezoid rule for the integral.)()(2311OggyyiiiiCode example The harmonic oscillation.Euler method:poor accuracy with t=0.02p.Predictor-corrector method:much better?/Predict the next position and velocity xi+1=xi+vi*dt;vi+1=vi-xi*dt;/Correct the new position and velocity xi+1=xi+(vi+vi+1)*dt/2;vi+1=vi-(xi+xi+1)*dt/2;Code example 4.2.Predictor-Corrector.cpp Another way to improve an algorithm is by increasing the number of mesh points j.Thus we can apply a better quadrature to the integral.For example,take j=2 and then use the linear interpolation scheme to approximate g(y,t)in the integral from gi and gi+1:)()(),(211Ogttgtttygiiiidttygyyjiittiji),(Now if we carry out the integration with g(y,t)given from this equation,we obtain a new algorithmwhich has an accuracy one order higher than that of the Euler algorithm.However,we need the values of the first two points in order to start this algorithm,because gi+1=g(yi+1,ti+1).)(2312Ogyyiii We can make the accuracy even higher by using a better quadrature.For example,we can take j=2 in above equation and apply the Simpson rule to the integral.Then we haveThis implicit algorithm can be used as the corrector if the previous algorithm is used as the predictor.)()4(35122OgggyyiiiiiA car jump over the yellow river 1997年，香港回归前夕，柯受良驾驶跑车成功飞越了黄河天堑壶口瀑布，长度达55米。飞越当天刮着大风，第一次飞越没有成功，但第二次成功了，其中有过很多危险的动作，但他都安全度过了，因此获得了“亚洲第一飞人”的称号。1953-2003 Let us take a simple model of a car jump over a gap as an example.The air resistance on a moving object is roughly given by ,where A is cross section of the moving object,r is the density of the air,and c is a coefficient that accounts for all the other factors.So the motion of the system is described by the equation setrcAfr.,ymgfmfadtvdvdtrd f is the total force on the car of a total mass m.Here y is the unit vector pointing upward.Assuming that we have the first point given,that is,r0 and v0 at t=0.Code example4.3.FlyingCar.cppThe RungeKutta method Formally,we can expand y(t+)in terms of the quantities at t with the Taylor expansion:A particle moving in one dimension under an elastic force We know the initial condition x(0),v(0).xkdtvdmamf x(t)=x(0)+x(0)t+x(0)t2/2+.v(t)=v(0)+v(0)t+v(0)t2/2+.x=v;x=v=-kx/m v=-kx/m=-kv/m The same process for higher orders xn and vn:x=v;v=-kv/m=k2x/m2;x=v;v=k2v/m2Code example 4th-order RungeKutta algorithm for the harmonic oscillator4.4.RungeKutta.cppChaotic dynamics nonlinear item nonlinear physics chaosAn LC circuit An LC circuit,also called a resonant circuit,tank circuit,or tuned circuit,is an electric circuit consisting of an inductor,represented by the letter L,and a capacitor,represented by the letter C,connected together.The circuit can act as an electrical resonator,an electrical analogue of a tuning fork,storing energy oscillating at the circuits resonant frequency.Equations of LC circuithttps:/en.wikipedia.org/wiki/LC_circuithttps:/ Use the 4th order Runge-Kutta method to solve a LC circuit with resistance&excitation.Boundary-value problems The solution of the Poisson equation with a given charge distribution and known boundary values of the electrostatic potential.Wave equations with given boundary conditions.The stationary Schrodinger equation with a given potential and boundary conditions.One-dimensional example Where u is a function of x,u and u are the 1st and 2nd derivatives of u with respect to x;f(u,u;x)is a function of u,u,and x.Either u or u is given at each boundary point.We can always choose a coordinate system so that the boundaries of the system are at x=0 and x=1 without losing any generality if the system is finite.);,(xuufu For example,if the actual boundaries are at x=x1 and x=x2 for a given problem,we can always bring them back to x=0 and x=1 by moving and scaling with a transformation:For problems in one dimension,we can have a total of four possible types of boundary conditions:(1)u(0)=u0 and u(1)=u1;(2)u(0)=u0 and u(1)=v1;(3)u(0)=v0 and u(1)=u1;(4)u(0)=v0 and u(1)=v1.(2)is the same as(3)by reversing the direction.)/()(121xxxxx The boundary-value problem is more difficult to solve than the similar initial-value problem with the differential equation.For example,if we want to solve an initial-value problem and the initial conditions u(0)=u0 and u(0)=v0,the solution will follow the algorithms discussed earlier.However,for the boundary-value problem,we know only u(0)or u(0),which is not sufficient to start an algorithm for the initial-value problem without some further work.The equation describing the stationary solution of elastic waves is .If both ends(x=0 and x=1)of the rod are fixed,the boundary conditions are u(0)=u(1)=0.If one end(x=0)is fixed and the other end(x=1)is free,the boundary conditions are u(0)=0 and u(1)=0.Example:longitudinal vibrations along an elastic rod For example,if both ends of the rod are fixed,the eigenfunctions are the possible solutions of the differential equation.Here the eigenvalues are given bywith l=1,2,.,.22)(plklThe shooting method The key here is to make the problem look like an initial-value problem by introducing an adjustable parameter;the solution is then obtained by varying the parameter.For example,given u(0)and u(1),we can guess a value of u(0)=a,where a is the parameter to be adjusted.The shooting method For a specific a,the value of the function u(1),resulting from the integration with u(0)=a to x=1,would not be the same as u1.The idea of the shooting method is to use one of the root search algorithms to find the appropriate that ensures f()=u(1)-u(1)=0 within a given tolerance.The shooting methodGuess a value for au(1)-ua(1)Calculate ua(1)Correct aconvergenotconvergeadjust a With given boundary conditions u(0)=0 and u(1)=1,We can define new variables y1=u and y2=u;Assume that this equation set has the initial values y1(0)=0 and y2(0)=.Here is a parameter to be adjusted in order to have f()=u(1)-1=0.We can combine the secant method for the root search and the 4th-order RungeKutta method for initial-value problems to solve the above equation set.Example)1(4,12221ydxdyydxdyp Many eigenvalue or boundary-value problems are in the form of linear equations,such as Assume that the boundary conditions are u(0)=u0 and u(1)=u1.If all d(x),q(x),and s(x)are smooth,we can solve the equation with the shooting method as shown above.However,an extensive search for the parameter a from f()=u(1)u1=0 is unnecessary in this case,because of the principle of superposition of linear equations:any linear combination of the solutions is also a solution of the equation.Linear equations We need only two trial solutions u0(x)and u1(x),where 0 and 1 are two different parameters.The correct solution of the equation is given by where a and b are determined from u(0)=u0 and u(1)=u1.Note that u0(0)=u1(0)=u(0)=u0.So we haveCode example 4.5.Boundary.cpp