计算物理作业2

计算物理作业

第一题:

a.用最小二乘法拟合下面的组数据

01234567

7.827.937.987.597.927.917.807.71

寻求经验公式，并拟合以上数据。

答：

matlab程序如下:

n=7;%n表示拟合的精度，在此取7

x=0:l:7;

y=[7.827.937.987.597.927.917.807.71];

al=polyfit(x,y,n);

xI=0:0.1:7;

yl=polyval(al,xl);

p]ot(x,y；*',xl,yl；-r');%作出x-y的散点图和xl-yl的拟合曲线

程序运行之后：

al-0.00240.0610-0.60733.0190-7.75769.4799-4.08277.8200

所以该组数据的经验公式就是：

y=7.82-4.0827X+9.4799/-7.7576X3+3.0190%4-0.6073%5+0.0610%6-0.0024%7

用matlab拟合的曲线

8.2

蓝色的散点图是x-y图，红色的多项式曲线就是拟合后的曲线。

当n取6或者更小时，拟合效果并没有上面的好，如下n=6时的拟合曲线所示:

5

第二题：

复化梯形计算定积分:

/=Jsinxdx

」o

要求：递交算法说明过程，源程序及实际结果。

答：

复化梯形的迭代公式为：

n-1

JV(x)dx=h/2f(a)+22f(“+f(b)

aj=l,

在这里,a=0,b=7i,©)=f(b)=0。

算法如下：

x=zeros(l,100);

y=zcros(1,100);

%x、y是两个一维零矩阵，用来存储不同的n和与之对应的梯形公式的定积分％

t=0;j=l;

forn=l:l:100;

fori=l:n-l;

%每个n对应的2年二卜（牛）的值赋值给t%

t=t+2*sin(i*pi/n);

end;

tl=(pi/(2*n))*t;

y(Lj)=ti;

x(lj)=n；%每个n（存储在x矩阵）对应的定积分值存储在y矩阵％

j=i+i;

t=0;%n值递增，t归零，j递增来将不同n对应的值y矩阵的不同位置％

plot(x,y);%作图x-y%

2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

0102030405060708090100

图梯形

算法在计算精度n不同时的取值

可以从matlab中读出y矩阵中的不同元素,比方n=10时,y=1.9835；n=10时,y=1.9959o

n=2-100,

1.999817732515361.999821510044761.999825171343901.99982872113273

1.999832163893991.99983550388744

可以看到当n取较小值时，梯形公式计算定积分的值逐渐接近精确值2而且是迅速的上升;

在n比拟大时，值接近平缓，也无限接近理论值2。

第三题：

仅用Romberg算法计算定积分

ln

/=Jsinxdx

"o

要求算法说明，源程序，最后结果，并与理论值比拟。

答：

Romberg迭代公式为：

Rkj=Rk,j一1+―$一1_1—

在Matlab中设计实现积分功能的M函数，然后在matlab对话框中输入要计算的函数，给

出区间和精度即可。

新建function文件，如下：

function[y]=romberg(f,n,a,b)

z=zeros(n,n);

h=b-a;

z(l,l)=(h/2)*(f(a)+f(b));

fl=0;

fori=2:n

fork=l:2A(i-2)

fl=fl+f(a+(k-0.5)*h);

end

h=h/2;

z(i,l)=0.5*z(i-1,1)+0.5*h*fl;

forj=2:i

A

z(i,j)=z(izj-l)+(z：izj-l)-z(i-l,j-1))/(4(j-l)-l);

end

end

y=z(n,n);

z

end

在matlab的命令窗口输入以下程序:

»clear

»f=inline('sin(x)');

»I二romberg(f,10,0,pi)

z=

Columns1through7

0.0000000000

0.78541.047200000

1.34081.52591.55780000

1.65751.76311.77891.7824000

1.82551.88151.88941.89121.891600

1.91201.94081.94471.94561.94581.94590

1.95581.97041.97241.97281.97291.97291.9729

1.97781.98521.98621.98641.98651.98651.9865

1.98891.99261.99311.99321.99321.99321.9932

1.99451.99631.99651.99661.99661.99661.9966

Columns8through10

000

000

000

000

000

000

000

1.986500

1.99321.99320

1.99661.99661.9966

1.9966

»

第四题：

用改良的欧拉公式，求解常微分方程初值问题

2

di=y

y(O.l)=1

'0.1<x<0,4

答：

+/(%”〃)

h(-、

%+i=%+司/(/，%)+/(芍+i，y“+i)

程序：

x=[0.10.120.14().16().18().20.220.240.26().28().3().320.340.36().380.4];

%取心0.02,代入梯形公式％

kzeros(l,16);%f为1*16维矩阵，用来存储y(x)的值％

ya=1；

f(l,l)=ya;

fori=2:16;

yb=ya+0.02*yaA2;

ya=ya+().()1*(yaA2+ybA2);

f(l,i)=ya；

end

plot(x,f；r+');%作出x-f的散点图，用红色+点表示％

holdon;

z=dsolve('Dy=yA2'；y(0.1)=1'；x');%解出y=f(x)的方程％

ezplot(z,[().1,0.4]);%作出z曲线％

ylabel('y');

legend('\itcaculate','\ittheory);

holdoff;%将散点图和曲线图组合到一个图表输出％

11/10)

由图可知计算值能很好的逼近理论值

第五题:

dyx

石=歹

用四阶龙格-库塔法求解微分方程y(2.o)=1

、2.0<x<2.6

答:

y*=+-(kl+2k2+2k3+44)

yi6

k1=hf(ti,yi)

左2=步&+0.5九、+0.50)

k3=hf(ti+0.5h,yi+0.5k2)

k4=hf(ti+用,％+左3)

程序：

x=[22.052.12.152.202.252.32.352.42.452.52.552.6];

%步长=().()5%

kzeros(l,13);%l*13矩阵，用以存储丫值％

f(U)=l；

fori=2:13;

k2=(x(l,i-l)+0.025)/(f(l,i-l)+0.025*k1);

k3=(x(I,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(l,i)=f(l,i-l)+(0.05/6)*(kl+2*k2+2*k3+k4);

end

plol(x,f,T4);%作出x-f的散点图，用红色+点表示％

holdon;

z=dsolve('Dy二x/y','y(2)=l','x');%解出y=f(x)的方程％

ezplot(z，l2.3]):%作出z曲线％

ylabel('f);

legend('\itcaculate','\ittheory)；

holdoff;%将散点图和曲线图组合到一个图表输出％

第六题：

干算3重积分：I=[\\n^zdxdydz

积分区间：平面x=0、平面y=0、平面片0、平面x+y+z=l所围区域

要求：要有源程序及结果及不同投点数的比拟

答：

设D二(x,y,z)是在V中均匀分布的随机变量，那么

,z)dP=夕/(x,y,z)|(x,y,z)^VyV

其中，E⑴为函数f在指定区间内的数学期望。

证明：由于D是在V中均匀分布的随机变量，因此其分布密度函数是

--DeV

PfD)二

L0otherwise

因此

£(，(0)=JV(0P(0〃=飘/(〃)〃=圳,f(x,y,z)dxdydz

即证

这里我们在x=O,x=l,y=O,y=l,z=O,z=1围成的空间Vo内投点N,设在V中点的个数为Ni,那么

V=Ni/N*Vo,同时选取N*3或3xN*l维随机数组进行投点。

程序：

j=h

k=100;%k是可以变化%

Na=zcros(1,k);

Ib=zcros(l,k);

forN=100:100:100*k%N值有k个进行循环％

Nl=0;

sum_f=0;

X=rand(N,l);

Y=rand(N,l);

Z=rand(N,l);%产生N个(0,1)的随机数％

V=0;

fori=l:N

ifX(i,l)+Y(i,l)+Z(i,l)<=l

N1=N1+1;

sum_f=sum_f+120*X(i,l)*Z(i,l);

end

end

V=N1/N;

E_f=sum_f/N1;

I=E_f*V;%I是每个N对应的1值％

Na(l,j)=N;

j=j+l；

end%j是用来控制I存储在lb的位置％

plotCNaJb,'*1);

k是可以随意调节的，下面比照k=100,1000的散点图

k=100

1.25

1.2

1.15

1.05

1

0.95

0.9

0.85

0.8

0.75

10

x104

k=1000

廿算值随投点数增加有逼近理论值1的趋势

第七题：

用牛顿差值公式计算y=/,x=2.15时的值，XG[2,2.4]

答：

工(X)=/(AQ)+九期,石［(Y一%)+・・・+/［玛,・・・,与］仁一%)・・・(Y-大工J

x=x^+ih(i=0,…，n)

当节点等距分布时:t

牛顿前插公式(一般当X靠近x()时用)

设x=x°+〃h则N.(x)=N.(x°+th)=E；"/(Xo)

f(〃+i“GA=OIAJ

A〃(x)=)…1T)…-“，^e(x0,xM)

牛顿后插公式(一般当X靠近x„时用)

将节点顺序倒置：

N8)=/(x〃)+/[为居1]住-工)+・・・+力丹・・・，引住一£)・・・仅一内)

设x=x〃+〃7,则M(x)=N〃(x〃+〃[)=£(-以：▽"(£)

k=QA

C++程序如下：

#include<math.h>

#includc<iostrcam.h>

#defineN5〃插值节点数目

voidmain()

(

doubletmp=1;

inlij;

doublex[N];〃插值节点坐标

doubley[NJ;

doubleg[N];

doubleu,newton;〃需要计算的插值点

coul<<"输入插值点坐标："«endl;

for(i=0;i<N;i++)

cin»y[i];

)

couiv〈”输入要计算的插值点：H«endl;

cin»u;

for(i=0;i<N;i++)

g；i]=yliJ；

for(i=0;i<N;i++)

(

forU=N;j>i;j-)

(

gU]=(gLj]-g[j-1])/(x[j]-x[j-i-1]);

)

)

newton=g[()l;

tmp=l;

for(i=l;i<=N;i++)

(

tmp*=(u-x[i-l]);

ncwton=ncwton+tmp*g[i];

)

cout«”所得结果为:"«ne\vton«endl;

)

程序运行结果:

/­D:\IDDOWNLOAD\Debug\newton.exe"

然后，输入差值点坐标，

21.4142132.11.4491382.21.4832402.31.5165752.41.549193

按回车

按回车，得到答案：1.46629,标准值：1.466287829

第八题：

1064〃m的脉冲激光作用于金属铝板，（厚度d=2mm）。求温度的分布以及演化的寸程。

根本要求：

1建立模型（方程和边界条件）

2查找所需的物理参数（考虑边界的热流，吸收系数）

3用有限差分法求解［用隐式的方法求解，一阶精度即可）

4使用。rigin软件给出某些时刻温度场分布、某些空间点温度变化曲线。

激光源作为边界条件处理，有定解问题:

S〃(x,t)hu(x,t)

-----------=Lf

dtaP

—k-----------=a/(力

dx

〃(x,0)=0

Q<x<2

经查阅相关文献，金属铝的密度夕=27(X)Ag/加3,铝材料对1064cm激光的吸收系数

。=0.052,比热c=78Q27+0.488T,热传导系数(其中T是温度，单位是K),方程中

D=K/cp

程序：

%激光加热问题

closeall

clearall

cic

k=30;

x=2;%材料长度mm

h=0.04;%空间步长mm

r=x/h;%空间点数

tt=l;%时间步长ns

absorb=0.052;%吸收系数

Ta=300;%室温

T0=300;%物体初始温度

IO=2Oe5;%激光中心处功率密度W

to=io；%脉冲宽度ns

density=0.27;%密度kg/mmA2

C=780.27+0.488*T0;%比热

ki=249.5-0.08*T0;%热导率

v=kr/(density*C);

a=(tt*v)/hA2;

A=zeros(r-l,r-l);%构造A矩阵

fori=2:l:r-l

A(iJ)=l+2*a;

end

%上边界条件

A(l,l)=l+a;

fori=l:r-2

A(i,i+l)=-a;

A(i+lJ)=-a;

end

U=zeros(r-l,k+l);%构造U

fori=l:(r-l)

U(i,l)=T0;%初值

end

fors=l:k

B=zeros(r-l,l)；%构造B矩阵

tx=s*tt;

F=(tx*exp(-tx/tO))/tO;

q=(h*IO*absorb*F)/kr;

B(l/l)=U(l,s)+a*q;%上边界条件

B(r-l/l)=U(19/s)+a*T0;%下边界条件

fori=2:r-2

B(i4)=U(bS);

end

x=A\B;%构造X矩阵

fori=l:r-l

if(x(i,l)>=T0)

U(i,s+l)=x(U);

else

U(i,s+1)=TO;

end

end

end

t=O:k;

x=0:h:2;

x=x(l:49);

x=x';

figure/mesh(t,x,U),view(-20,10)

xlabel('t(ns)','rotation',2)

ylabel('x(mm)7rotation'z-20)

zlabel('T(K)')

title。温度场分布)

计算结果：

温度场分布

65cli

600-

550-

500^

g

450-

30

t(ns)

Mmm)

Kns)

第九题：

激光脉冲宽度为10ns,时间分布为/(r)=Lexp(-L),波长为1064nm,能量密度为

，010

10mJ/mm：作用于铝板(厚度为2mm，长10mm)的中心区域产生的温度场分布及演化过

程。根本要求：

1建立模型；

2查找物理参数；

3有限差分方法或有限元软件(COMSOL求解)；

4给出温度场分布和温度变化曲线。

2.模型建立:

儿何模型如下:

材料参数如下:

属性名称值单位

V热膨胀系数alpha8e-6[l/K]1/K

V常压热容Cp900[J/(kg*K)]J/(kg-K)

，密度rho39OO[kg/mA3]kg/m3

V导热系数k27[W/(m*K)]W/(m-K)

y杨氏模量E300e9[Pa]