计算物理作业 2

1、计算物理作业第一题：a. 用最小二乘法拟合下面的一组数据012345677.827.937.987.597.927.917.807.71寻求经验公式，并拟合以上数据。答：matlab程序如下：n=7; % n表示拟合的精度，在此取7x=0:1:7;y=7.82 7.93 7.98 7.59 7.92 7.91 7.80 7.71;a1=polyfit(x,y,n);x1=0:0.1:7;y1=polyval(a1,x1);plot(x,y,*,x1,y1,-r); %作出x-y的散点图和x1-y1的拟合曲线程序运行之后：a1 -0.00240.0610-0.60733.0190-7.75769

2、.4799-4.08277.8200所以该组数据的经验公式就是：y=7.82-4.0827x+9.4799x2-7.7576x3+3.0190x4-0.6073x5+0.0610x6-0.0024x7用matlab拟合的曲线蓝色的散点图是x-y图，红色的多项式曲线就是拟合后的曲线。当n取6或者更小时，拟合效果并没有上面的好，如下n=6时的拟合曲线所示：b在某次实验中需要观察水分的渗透速度，测得时间t与水重量w的数据t1248163264w4.224.023.854.593.443.022.59已知t与w关系 w=Ats，试用最小二乘法确定A、S。答：先对式子两边取对数，化为一阶，然后使用上题的

3、一阶拟合的程序，取n=1t=1 2 4 8 16 32 64;w=4.22 4.02 3.85 4.59 3.44 3.02 2.59;x=log(t);y=log(w);a1=polyfit(x,y,1);A=exp(a1(2);S=a1(1);x1=1:0.1:64;y1=A*x1.S;plot(t,w,*,x1,y1,-r);程序运行结果：a1 -0.1107 1.5153因此，A=e1.5153=4.5509 S=-0.1107拟合曲线：第二题：复化梯形计算定积分：I=0sinxdx要求：递交算法说明过程，源程序及实际结果。答：复化梯形的迭代公式为：abf(x)dx=h/2fa+2j=

4、1n-1fxj+f(b)在这里，a=0，b=，fa=fb=0。算法如下：x=zeros(1,100);y=zeros(1,100); %x、y是两个一维零矩阵，用来存储不同的n和与之对应的梯形公式的定积分%t=0;j=1;for n=1:1:100;for i=1:n-1;t=t+2*sin(i*pi/n); %每个n对应的2j=1n-1fxj的值赋值给t%end;t1=(pi/(2*n)*t;y(1,j)=t1;x(1,j)=n; %每个n（存储在x矩阵）对应的定积分值存储在y矩阵%j=j+1;t=0; %n值递增，t归零，j递增来将不同n对应的值y矩阵的不同位置%end;plot(x,y)

5、; %作图x-y% 图 梯形算法在计算精度n不同时的取值可以从matlab中读出y矩阵中的不同元素，比如n=10时，y=1.9835；n=10时，y=1.9959。n=2100,1.570796326794901.813799364234221.896118897937041.933765598092811.954097233313711.966316678765891.974231601945551.979650811216481.983523537509451.986386986581661.988563776584321.990257175347771.991600427355071.9

6、92683831530771.993570343772341.994304944309461.994920463583451.995441318320191.995885972708721.996268598739461.996600220269271.996889516466771.997143395803951.997367412545631.997566073264041.997743065358541.997901429465681.998043690970561.998171961343651.998288016964171.998393360970141.9984892721876

7、01.998576844134181.998657016333451.998730599624851.998798296749971.998860719196551.998918401057861.998971810497071.999021359278081.999067410726891.999110286411871.999150271773441.999187620887651.999222560512861.999255293540041.999286001945131.999314849324061.999341983076261.999367536291511.999391629

8、385201.999414371519651.999435861842901.999456190571241.999475439937671.999493685024891.999510994498601.999527431254561.999543052990811.999557912714681.999572059193141.999585537353381.999598388640041.999610651334151.999622360838611.999633549933971.999644249008121.999654486262881.999664287899991.99967

9、3678288961.999682680118711.999691314534761.999699601263501.999707558724961.999715204135311.999722553600021.999729622198761.999736424062841.999742972445841.999749279788281.999755357776721.999761217397971.999766868988701.999772322281151.999777586445001.999782670125971.999787581481331.999792328212631.9

10、99796917595941.999801356509761.999805651460761.999809808607661.999813833783361.999817732515361.999821510044761.999825171343901.999828721132731.999832163893991.99983550388744可以看到当n取较小值时，梯形公式计算定积分的值逐渐接近精确值2而且是迅速的上升；在n比较大时，值接近平缓，也无限接近理论值2。第三题：使用Romberg算法计算定积分I=0sinxdx要求算法说明，源程序，最后结果，并与理论值比较。答：Romberg迭代

11、公式为：Rk,j=Rk,j-1+Rk,j-1Rk-1,j-14j-1-1在Matlab中设计实现积分功能的M函数，然后在matlab对话框中输入要计算的函数，给出区间和精度即可。新建function文件，如下：function y = romberg( f,n,a,b )z=zeros(n,n);h=b-a;z(1,1)=(h/2)*(f(a)+f(b);f1=0;for i=2:n for k=1:2(i-2) f1=f1+f(a+(k-0.5)*h); end h=h/2; z(i,1)=0.5*z(i-1,1)+0.5*h*f1; for j=2:i z(i,j)=z(i,j-1)+(z

12、(i,j-1)-z(i-1,j-1)/(4(j-1)-1); endendy=z(n,n);z end在matlab的命令窗口输入以下程序： clear f=inline(sin(x); I=romberg(f,10,0,pi)z = Columns 1 through 7 0.0000 0 0 0 0 0 0 0.7854 1.0472 0 0 0 0 0 1.3408 1.5259 1.5578 0 0 0 0 1.6575 1.7631 1.7789 1.7824 0 0 0 1.8255 1.8815 1.8894 1.8912 1.8916 0 0 1.9120 1.9408 1.9

13、447 1.9456 1.9458 1.9459 0 1.9558 1.9704 1.9724 1.9728 1.9729 1.9729 1.9729 1.9778 1.9852 1.9862 1.9864 1.9865 1.9865 1.9865 1.9889 1.9926 1.9931 1.9932 1.9932 1.9932 1.9932 1.9945 1.9963 1.9965 1.9966 1.9966 1.9966 1.9966 Columns 8 through 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.9865 0 0 1.9

14、932 1.9932 0 1.9966 1.9966 1.9966I = 1.9966第四题：用改进的欧拉公式，求解常微分方程初值问题dydx=y2y0.1=10.1x0.4答：程序：x=0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4;%取h=0.02，代入梯形公式%f=zeros(1,16);%f为1*16维矩阵，用来存储y（x）的值%ya=1;f(1,1)=ya;for i=2:16; yb=ya+0.02*ya2; ya=ya+0.01*(ya2+yb2); f(1,i)= ya;en

15、d plot(x,f,r+); %作出x-f的散点图，用红色+点表示%hold on;z=dsolve(Dy=y2,y(0.1)=1,x); %解出y=f（x）的方程%ezplot(z,0.1,0.4); %作出z曲线%ylabel(y);legend(it caculate,it theory);hold off; %将散点图和曲线图组合到一个图表输出%由图可知计算值能很好的逼近理论值第五题：用四阶龙格-库塔法求解微分方程dydx=xyy2.0=12.0x2.6答：程序：x=2 2.05 2.1 2.15 2.20 2.25 2.3 2.35 2.4 2.45 2.5 2.55 2.6;%步

16、长=0.05%f=zeros(1,13);%1*13矩阵，用以存储y值%f(1,1)=1;for i=2:13; k1=x(1,i-1)/f(1,i-1); k2=(x(1,i-1)+0.025)/(f(1,i-1)+0.025*k1); k3=(x(1,i-1)+0.025)/(f(1,i-1)+0.025*k2); k4=(x(1,i-1)+0.05)/(f(1,i-1)+0.05*k3); f(1,i)=f(1,i-1)+(0.05/6)*(k1+2*k2+2*k3+k4);endplot(x,f,r+); % 作出x-f的散点图，用红色+点表示%hold on;z=dsolve(Dy=

17、x/y,y(2)=1,x); %解出y=f（x）的方程%ezplot(z,2,3); %作出z曲线%ylabel(f);legend(it caculate,it theory);hold off; %将散点图和曲线图组合到一个图表输出%第六题：计算3重积分：积分区间：平面x=0、平面y=0、平面z=0、平面x+y+z=1所围区域 要求：要有源程序及结果及不同投点数的比较答：设D=(x,y,z)是在V 中均匀分布的随机变量，则其中，E(f)为函数f在指定区间内的数学期望。证明：由于D是在V中均匀分布的随机变量, 因此其分布密度函数是p(D)= 因此即证这里我们在x=0,x=1,y=0,y=1,

18、z=0,z=1围成的空间V0内投点N，设在V中点的个数为N1，则V=N1/N*V0，同时选取N*3或3N*1维随机数组进行投点。程序：j=1;k=100; %k是可以变化%Na=zeros(1,k);Ib=zeros(1,k);for N=100:100:100*k %N值有k个进行循环%N1=0;sum_f=0;X=rand(N,1);Y=rand(N,1);Z=rand(N,1); %产生N个（0,1）的随机数%V=0;for i=1:N if X(i,1)+Y(i,1)+Z(i,1)=1 N1=N1+1; sum_f=sum_f+120*X(i,1)*Z(i,1); endendV=N1

19、/N;E_f=sum_f/N1;I=E_f*V; %I是每个N对应的I值%Na(1,j)=N;Ib(1,j)=I;j=j+1;end %j是用来控制I存储在Ib的位置%plot(Na,Ib,*);k是可以随意调节的，下面对比k=100,1000的散点图k=100k=1000计算值随投点数增加有逼近理论值1的趋势第七题：用牛顿差值公式计算y=x，x=2.15时的值，x2,2.4答：当节点等距分布时: 牛顿前插公式（一般当 x 靠近 x0 时用）牛顿后插公式（一般当 x 靠近 xn 时用）c+程序如下：#include #include #define N 5 /插值节点数目void main()

20、double tmp=1;int i,j;double xN; /插值节点坐标double yN;double gN;double u,newton; /需要计算的插值点cout输入插值点坐标：endl;for(i=0;ixi;cinyi;cout输入要计算的插值点：u;for(i=0;iN;i+)gi=yi;for(i=0;ii;j-)gj=(gj-gj-1)/(xj-xj-i-1);newton=g0;tmp=1;for(i=1;i=N;i+)tmp*=(u-xi-1);newton=newton+tmp*gi;cout所得结果为：newton=T0) U(i,s+1)=x(i,1); else U(i,s+1)=T0; end endendt=0:k;x=0:h:2;x=x(1:49);x=x;figure,mesh(t,x,U),view(-20,10)xlabel(t(ns),rotation,2)ylabel(x(mm),rota