计算物理ComputationalPhysics计算物理 (4).pdf
Isaac NewtonGottfried Leibniz Computational physiCsNumerical calculusxpositionvvelocityaaccelerationderivative derivative integralintegralNumerical calculus Numerical differentiation Numerical integration Roots of an equation Extremes of a functionNumerical differentiationTaylor expansion:To calculate f,f,f.single-variable:)(!)()(!2)()()()()(0)(0020000 xfnxxxfxxxfxxxfxfnnmultivariable:),(!2)(2),(!2)(),(!2)(),()(),()(),(),(00000020002000000000yxfyyxxyxfyyyxfxxyxfyyyxfxxyxfyxfxyyyxxyxyxffxy/2Now if we divide the space into discrete points xi with evenly spaced intervals xi+1-xi=h and label the function at the lattice points as fi=f(xi),we obtain the simplest expression for the first-order derivative.x)f(xx)f(x)(xfii0 xilimThe first-order derivative of a single-variable function f(x)around a point xi is defined from the limit.the two-point formula for the first-order derivative)(1hOhfffiiiAn improved choice:The accuracy is improved from O(h)to O(h2).3/0206/2/6/2/311321321iiiiiiiiiiiiiifhfhfffhfhfhfffhfhfhffThe three-point formula for the first-order derivative)(O2211hhfffiii)(O2311hfhffiii)()88(121)2(8)1(42112hOffffhfiiiii(1)(2)A five-point formula can be derived by including the expansions of fi+2 and fi-2 around xi.)(O325)3(311hfhfhffiiii)(O3845)3(322hfhfhffiiiiMore pointsHigher accuracySmaller hnumber of pointsinaccuracy2O(h)3O(h2)5O(h4)SummarySimilarly,we can combine the expansions of fi2 and fi1 around xi and fi to cancel the fi,f(3)i,f(4)i,and f(5)i terms;then we havethe five-point formula for the second-order derivative)(O163016121421122hfffffhfiiiiii ）（The three-point formula for the second-order derivative 6/2/6/2/321321iiiiiiiiiifhfhfhfffhfhfhff)(O2)(O222114211hhffffhfhfffiiiiiiii Example Given f(x)=sin(x),lets calculate f(x)&f(x).Divide the region from 0 to p/2 to 100 equal-length intervals with 101 points i*p/200(i=0,1,2.100).For boundary points(i=0,1&99,100),we can use Lagrange interpolation to extrapolate the derivatives.Code example:3.1.Differentiation.cpp Three-point formula for fThree-point formula for fNumerical calculus Numerical differentiation Numerical integration Roots of an equation Extremes of a functionNumerical integrationDivide the region a,b into n slices,evenly spaced with an interval h.If we label the data points as xi with i=0,1,.,n,we can write the entire integral as a summation of integrals,with each over an individual slice.In general,we want to obtain the numerical value of an integral,defined in the region a,b,.)(dxxfSba 101)()(nixxbaiidxxfdxxf)()(2/)()(21011hOffhShffxxfxfniiiiiiiThe above quadrature is commonly referred as the trapezoidal rule,which has an overall accuracy up to O(h2).The simplest quadrature is obtained if we approximate f(x)in the region xi,xi+1 linearly,that is,Trapezoidal rule We can obtain a quadrature with a higher accuracy by working on two slices together.If we apply the Lagrange interpolation to the function f(x)in the region xi-1,xi+1,we have)()()()()()()()(311111111111111hOfxxxxxxxxfxxxxxxxxfxxxxxxxxxfiiiiiiiiiiiiiiiiiiiii)()4(3412/022122hOfffhSnjjjjExample Given f=sin(x),integrate f from 0 to p/2.The analytic function:-cos(x)The exact value:cos(0)-cos(p/2)=1.We can use trapezoidal rule to see how the numerical value converges to 1.Code example 3.2.Integration.cppHomework Improved integration with the three-point Lagrange interpolation implemented.Comparison with the trapezoidal rule method.)()4(3412/022122hOfffhSnjjjjNumerical calculus Numerical differentiation Numerical integration Roots of an equation Extremes of a functionRoots of an equationIf we need to find a root for f(x)=a,then how?define g(x)=f(x)-a,and find a root for g(x)=0.In physics,we often encounter situations in which we need to find the possible value of x that ensures the equation f(x)=0,where f(x)can either be an explicit or an implicit function of x.If such a value exists,we call it a root or zero of the equation.Bisection method If we know that there is a root xr in the region a,b for f(x)=0,we can use the bisection method to find it within a required accuracy.Most intuitive method.f(a)f(b)0 x0=(a+b)/2|f(x0)|d?output!yesnof(a)f(x0)0 f(x0)f(b)0f(a)f(x0)0 or f(x0)f(b)0?b=x0 a=x0Code Example f(x)=sin(x)=0.5;x is within 0 to p/2.Analytically,we know the root is p/6.Numerically,the procedure is:since sin(0)-0.5*sin(p/2)-0.50 and sin(0)-0.5*sin(p/4)-0.50;then the root must be within(0,p/4).Then we calculate the value at p/8.3.3.Bisection.cppThe Newton method0)()()()(xfxxxfxfrrwhere x can be viewed as a trial value for the root of xi at the ith step and the approximate value of the next step xi+1 can be derived.This method is based on linear approximation of a smooth function around its root.We can formally expand the function f(xr)=0 in the neighborhood of the root xr through the Taylor expansion.(i=0,1,.)iiiiiiffxxxx/1Code example 3.4.NewtonRoot.cppExample:f(x)=sin(x)=0.5;g(x)=f(x)-0.5=sin(x)-0.5ixigigi00-0.5110.5Possible bugs If the function is not monotonous If fi=0 or very small at some points Works well when the function is monotonous,especially with moderate f.Secant method-discrete Newton methodIn many cases,especially when f(x)has an implicit dependence on x,an analytic expression for the first-order derivative needed in the Newton method may not exist or may be very difficult to obtain.We have to find an alternative scheme to achieve a similar algorithm.One way to do this is to replace the analytic f(x)with the two-point formula for the first-order derivative,which gives)/()(111iiiiiiifffxxxxCode example3.5.Secant.cppExample:f(x)=sin(x)=0.5;g(x)=f(x)-0.5=sin(x)-0.5ixigi00-0.51p/20.52p/4)/()(111iiiiiiigggxxxx4)5.05.0/(5.0)02(22pppxNumerical calculus Numerical differentiation Numerical integration Roots of an equation Extremes of a functionExtremes of a function An associated problem to find the root of an equation is finding the maxima and/or minima of a function.Examples of such situations in physics occur when considering the equilibrium position of an object,the potential surface of a field,and the optimized structures of molecules and small clusters.We know that an extreme of g(x)occurs at the point with which is a minimum(maximum)if f(x)=g(x)is greater(less)than zero.So all the root-search schemes discussed so far can be generalized here to search for the extremes of a single-variable function.0)()(dxxdgxfExample)exp(4)(0002rrVrerVpwhere e is the charge of a proton,0 is the electric permittivity of vacuum,and V0 and r0 are parameters of this effective interaction.The first term comes from the Coulomb interaction between the two ions,but the second term is the result of the electron distribution in the system.The(ionic)bond length of the diatomic moleculeNaClThe force:At equilibrium,the force between the two ions is zero.Therefore,we search for the root of f(x)=-dV(x)/dx=0.)exp(4)()(000202rrrVredrrdVrfpcode example parameters for NaCl e2/4p0=14.4 AeV V0=1.09 x 103 eV r0=0.33 A r starts from 1 A 3.6.NaCl.cpp/iiiffxIn the example program above,the search process is forced to move along the direction of descending the function g(x)when looking for a minimum.In other words,for xi+1=xi+xi,the increment xi has the sign opposite to g(xi).Based on this observation,an update scheme can be formulated:with a being a positive,small,and adjustable parameter.For the minimum,f(or g)must be positive.iiifagax This scheme can be generalized to the multivariable case asNote that step xi here is scaled by|g(xi)|and is forced to move toward the direction of the steepest descent.This is why this method is known as the steepest-descent method.where x=(x1,x 2,.,x l)and g(x)=(g/x1,g/x 2,.,g/xl).|)(|/)(1iiiiiixgxgaxxxxHomework Search for the minimum of the function g(x,y)=sin(x+y)+cos(x+2*y)in the whole space
