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

《计算物理ComputationalPhysics计算物理 (3).pdf》由会员分享，可在线阅读，更多相关《计算物理ComputationalPhysics计算物理 (3).pdf（32页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、Computational physiCs Approximation of functions Linear interpolation Lagrange interpolation Newton interpolation Linear system method Least-squares approximation Millikan experimentWhat is interpolation?Interpolation is needed when we want to infer some local information from a set of incomplete or

2、 discrete data.trajectory of a golf ball/football/missile/.Interpolation between two points Linear interpolation y=y(x)More points:Direct connections between two nearest neighbor points.010010 xxxxyyyyMore smooth interpolation Lagrange interpolation Lets start from the simplest case:1100)()()(yxAyxA

3、xy0)(1)(1000 xAxA1)(0)(1101xAxA(x0,1)(x1,0)1010 xxxxA(x1,1)(x0,0)0101xxxxALagrange interpolation-Three pointsy1y2x0 x1XYOy=f(x)x2y0221100)()()()(yxAyxAyxAxy)()()()()()(120210201210212010210 xxxxxxxxAxxxxxxxxAxxxxxxxxA111x0 x1x2x0 x1x2x0 x1x2)(0 xA)(1xA)(2xALagrange interpolation-a general formula In

4、 general,for n+1 points:jnjjyxAxy)()(0njiiijijxxxxxA0)(The principlex=xjAj(x)=1x=xi(i!=j)Aj(x)=0How to write a code?A subroutine Inputs:(xj,yj),and x Output:y(x)Algorithm in the black box:To calculate the coefficient Aj(x)Then we can obtain:jnjjyxAxy)()(0double interpolate(const double x,const doubl

5、e y,const int n,const double xx)/n+1:total points/x:array of xi/y:array of yi/xx:x in y(x)double yx=0;/y=y(x)for(int j=0;j=n;j+)/Aj*yjdouble aj=1;/Ajfor(int i=0;i a1:y1-y0/(x1-x0)N2(x2)=y2 -a2:y2-N1(x2)/(x2-x0)(x2-x1)N3(x3)=y3 -a3:y3-N2(x3)/(x3-x0)(x3-x1)(x3-x2).)()()()()()()(210323102120101xxxxxxaN

6、xNxxxxaNxNxxayxN)()()()()(2103102010 xxxxxxaxxxxaxxayxyHomework I Newton interpolation of 10 equal spacing points of cos(x)within 0,p.One more method-self-madeThe fact is:y(x)=f0+f1x+f2x2+f3x3+f4x4.+fnxna polynomial function for n+1 points.Then we get a system of linear equations involving the same

7、set of variables fk(k=0n).The coefficients are xik(k=0n).y(x)=f0+xf1+x2f2+x3f3+x4f4.+xnfn y0=y(x0)=a00f0+a01f1+a02f2+a03f3+.+a0nfn.yn=y(xn)=an0f0+an1f1+an2f2+an3f3+.+annfn here aij=xij.Totally n+1 equations with n+1 variables.n0n0nnn3n2n1n00n03020100yyffaaaaaaaaaa A*f=B;A is a matrix constructed by

8、elements aij.B is a column vector of yi.By solving the linear system,we can obtain the values of fis,then any y(x)can be calculated.Three points example Three points:(x0,y0),(x1,y1)(x2,y2)The interpolation function is y(x)=f0+f1x+f2x2.The linear equations are:y0=f0+f1x0+f2x02 y1=f0+f1x1+f2x12 y2=f0+

9、f1x2+f2x22210210222211200111yyyfffxxxxxxvoid interpolate(const double x,const double y,const int n,double f)/n:total points/x:array of x/y:array of y/f:array of f double ann;/Matrix Afor(int i=0;in;i+)/Aj*yj double xj=1;const double xi=xi;for(int j=0;j Overall approximation or fitting A typical exam

10、ple is a polynomial fit to a set of experimental data with error bars.Process What do we have?A fitting function p(x)vs the data(xi,yi).p(x)is close to but may not pass through(xi,yi).The differences between p(xi)and yi are the error bars.The aim is to reduce the error bars to a minimum level by adj

11、usting the coefficients of p(x).Example:p(x)is a mth-order polynomialnm,or the error bars can be zerom+1 variableskkmkmxaxp0)(for discrete data202)(iimnikyxpaTo minimize the error bars.We get m+1 linear equations.Then solve the linear system for ak.The simplest example:linear fitting02lkaa21002101)(

12、)(iinikyxaaaxaaxpBy calculating the partial derivative311002100)1(cacaccacan02lkaaiiniiniiniiniyxcycxcxc030220100,where.)1()1(,)1(120320112021300cnccnccacncccccaJohann Carl Friedrich Gauss A German mathematician who contributed significantly to many fields,including number theory,algebra,statistics,

13、analysis,differential geometry,geodesy,geophysics,mechanics,electrostatics,astronomy,matrix theory,and optics.Referred to as the Prince of Mathematicians and greatest mathematician since antiquity.1777-1855Robert Andrews MillikanIn 1910,Millikan published his famous work on the oil drop experiment i

14、n Science.Based on the measurements of the charges carried by all the oil drops,Millikan concluded that the charge carried by any object is a multiple(with a sign)of a fundamental charge,the charge of an electron(for negative charges)or the charge of a proton(for positive charges).1868 1953Hermann v

15、on Helmholtz-Albert Abraham Michelson(1st)-Robert Andrews Millikan(2nd)-Chung-Yao Chao(赵忠尧)President of the California Institute of TechnologyPresident of the American Physical Societyfor his work on the elementary charge of electricity and on the photoelectric effect(1923)Input dataThe Millikan exp

16、erimentFitting equation:qk=k*Qe+Q0Output:Qe&Q0Millikans original oil-drop apparatus,1909-1910void Millikan(const double k,const double q,const int n)/n:total points/k:array of k/q:array of q double c4=0,0,0,0;/coefficientsfor(int i=0;in;i+)c0+=ki;c1+=ki*ki;c2+=qi;c3+=ki*qi;const double q0=(c0*c3-c1*

17、c2)/(c0*c0-n*c1);const double qe=(c0*c2-n*c3)/(c0*c0-n*c1);coutThe elementary charge is:tqeendl;Code example.)1()1(,)1(120320112021300cnccnccacncccccaCode example 2.2.Millikan.cppNonlinear systems-I For example:the physical behavior is in the form of p(x)=a*expb*x,like nuclear decay.If you fit data

18、directly,the huge contrast between large x and small x regions will make the fitting inaccurate.An alternation is to fit lny and lnp(x).lnp(x)=lna+b*x,which is a linear function.Nonlinear systems-II For example:the physical behavior is in the form of p(x)=a/(x-b),which diverges at x=b.An alternation is to fit 1/p(x).1/p(x)=x/a-b/a,which is a linear function.Homework II Least-squares approximation for nuclear decay.To find the half-life of an unknown nucleus.time0-11-22-33-44-55-6fission300200150906050You can do it as a project!