大物上机实验报告

大学物理I上机实验报告

                         实验一

实验名称:同方向不同频率简谐振动的合成,“拍”

实验原理:设存在两列振动频率分别为w1和w2,振幅均为A。由于两列波频率不同,所以存在一定的周期使得两列波的旋转矢量在坐标系上重合,两列波的位相相同。反映在实际振动过程中为两列波的位移相同其位移的合成直接相加即可。为方便研究,在此处我们记两列波位矢重合的时刻为零点。两列波的数学表达式为

x1=Acos(w1t+Φ0)

x2=Acos(w2t+Φ0)

根据数学表达式合成可得到合振动的数学表达式为

x=x1+x2=

根据合振动的数学表达式我们可以看到对于合振动不是简谐振动,但仍是周期运动。通过观察表达式可以发现前一余弦表达式的角频率较后一余弦表达式周期的角频率小。当w2、w1均较大而其差值相差较小时更是如此。下面我将2A*cos((w2-w1)*t/2)一同视作振幅的表达式。并探讨w2、w1的大小关系与合振动图像之间的关系。

 实验现象

实验现象1:

w1、w2均较大时,实验数据:A=10,w1=2011,w2=2012,x=30。如图1

实验现象2:

w1、w2均较小时,实验数据:A=10,w1=11,w2=12,x=30。如图2。

实验现象3:

w1较大,w2较小时,实验数据:A=10,w1=2012,w2=12,x=30。如图3。

   

               图1                              图2                         图3

对图描述:

对于前面的振幅项,我们可以发现它的贡献在于对于“宏观振幅”的影响。从上述三幅图中可以看到“宏观振幅”呈现周期变化。3个图的振幅都一样,图1图线密集,图2较稀疏,图3因w的不同而产生不同的图线。

实验结果分析与讨论:

由上述三种情况及实验结果图像我们可以发现在同方向不同频率简谐波振动的合成实验中我们得到了实验所需要观察得到的“拍”的现象。

                                    实验二

实验名称:阻尼与受迫振动

实验原理:简谐振动是一种理想的运动模型,它是一种等幅振动,其机械能始终守恒,也不会以波的形式向外辐射能量引起能量损耗。在此实验中我们将考虑振动系统受到阻力影响的情况下系统的振动情况。在实验中,设系统受到的阻尼表达式为f=-γV,γ为正的比例常数。它的大小由物体的形状、大小、表面情况及周围介质的性质决定。V表示系统中谐振子的运动速度。易得下列动力学方程:

         (其中β=γ/2m,w^2=k/m)   。。。。。。。。。。。。。。。。 式2-1

由上述动力学方程及二阶线性微分方程求解可得系统的振动方程的具体形式将由β与w的大小关系决定。

在实际的振动系统中,由于阻尼的影响,为了获得稳定而又持续的振动,通常需要给阻尼振动系统一个周期性的外力,以弥补阻尼振动过程中损失的能量,这种周期性的外力称为策动力,在策动力作用下的振动称为受迫振动。

设策动力为Fcos(pt)。则式2-1可表达为

(其中β=γ/2m,w^2=k/m)   。。。。。。。。。。式2-2

故由线性非齐次微分方程的求解可知,方程的解由暂态解和稳态解构成。同时有稳态解的形式Fcos(pt)相同。

综上可解得

 

 

实验将考虑稳态时系统的振动情况,即x=a*cos(pt+θ)的相关情况。

实验现象:

阻尼实验:

受迫振动实验:

该实验分为3种情况;

情况1:f(F)〈f(固)(图一:阻尼较大情况下,图四:阻尼较小情况下)。

情况2:f(F)=f(固)(图二:阻尼较大情况下,图五:阻尼较小情况下)。

情况3:f(F)〉f(固)(图三:阻尼较大情况下,图六:阻尼较小情况下)。

图一                    图二                  图三

图四                   图五                   图六

对图描述:

在欠阻尼的图像中,系统是在震荡下恢复平衡的,随着阻尼系数的增大,震荡情况减弱而变为曲线下降,使得系统达到稳定时间加快;

在临界阻尼的图像中,振子的运动不再有周期性,但能较快回到平衡位置;

在过阻尼的图像中,振子的运动不再有周期性,且只能较缓慢地回到平衡位置。

在受迫振动实验中发现在三种不同的策动力频率下,系统恢复振动平衡后的振幅有变化。

  

实验结果分析及讨论:

在本次实验中,自己对阻尼振动,受迫振动进行了探究,在阻尼振动实验中发现在三种不同阻尼情况下系统的振动情况差异较大,更具体表现在了系统恢复平衡状态所需的时间上,按时间从大到小依次排列是欠阻尼,临界阻尼,过阻尼。

在受迫振动实验中发现在三种不同的策动力频率下,系统恢复振动平衡后的振幅有变化。

看来阻尼的存在大小和策动力的频率将会给振动系统带来很大的影响。

         

                                     

                         实验三

实验名称:相互垂直的简谐振动的合成

实验原理:

  一个质点同时参与两个相互垂直方向的简谐振动,一个沿x轴方向,另一个沿y轴方向,其振动频率相同,振动方程分别为

                             x=A1cos(wt+φ1)

                             y=A2cos(wt+φ2)

易得——合振动轨迹方程:

             x2/A12+y2/A22-2xycos(φ21)/A1A2=sin221)

    所以,质点合振动的轨迹一般为一椭圆。

实验现象:

     椭圆轨道不会超出以2A1和2A2为边长的矩形范围。当2个分振动的振幅A1和A2给定时,椭圆的形状由△φ=φ21决定。

 

对图描述:

  当△φ=kπ(k=0,1,2.........)时,合振动为一条经过一三象限或者二四象限的直线,仍然为简谐振动;当△φ≠kπ时,合振动为一个椭圆。

  按图示顺序,合振动轨迹由直线变为椭圆,再由椭圆变为直线。

         

实验结果分析及讨论:

   在实验中,两个互相垂直的简谐振动的位相差决定着合振动的图形,其图形一般为椭圆。若2个分振动的位相差取其它值,则猜想合振动的轨迹将为形状和方位各不相同的椭圆,质点的运动方向则可能为逆时针或顺时针。

 

第二篇:数值代数上机实验报告

试验项目名称:平方根法与改进平方根法

实验内容:先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序,然后用你编写的程

序求解对称正定方程组

Ax=b,其中,

A=[10

1 10 1

1 10 1

1 10]100*100

b随机生成,比较计算结果,评论方法优劣。

实验要求:平方根法与改进的平方根的解法步骤;存储单元,变量名称说明;系数矩阵

与右端项的生成;结果分析。

实验报告 姓名:罗胜利 班级:信息与计算科学0802 学号:实验一、平方根法与改进平方根法

先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序,然后用你编写的程序求解对称正定方程组AX=b,其中 系数矩阵为40阶Hilbert矩阵,即系数矩阵A的第i行第j列元素为=,

向量b的第i个分量为=.

平方根法函数程序如下:

function [x,b]=pingfanggenfa(A,b)

n=size(A);

n=n(1);

x=A^-1*b; %矩阵求解 disp('Matlab自带解即为x');

for k=1:n

A(k,k)=sqrt(A(k,k));

A(k+1:n,k)=A(k+1:n,k)/A(k,k);

for j=k+1:n;

A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);

end

end

for j=1:n-1

b(j)=b(j)/A(j,j);

b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);

end

b(n)=b(n)/A(n,n);

A=A';

for j=n:-1:2

b(j)=b(j)/A(j,j);

b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);

end

b(1)=b(1)/A(1,1);

disp('平方根法的解即为b');

end

function [x]=ave(A,b,n)

L=zeros(n,n);

D=diag(n,0);

S=L*D;

for i=1:n

L(i,i)=1; %Cholesky分解 %前代法 %回代法 %用改进平方根法求解Ax=b %L为n*n矩阵 %D为n*n的主对角矩阵 %L的主对角元素均为1

end

for i=1:n

for j=1:n %验证A是否为对称正定矩阵

if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时,输出wrong disp('wrong');break;end

end

end

D(1,1)=A(1,1);

for i=2:n

for j=1:i-1

S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');

L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);

end

D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');

end

y=zeros(n,1);

x=zeros(n,1);

for i=1:n

y(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i);

end

for i=n:-1:1

x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %将A分解使得A=LDLT % x,y为n*1阶矩阵 %通过LDy=b解得y的值 %通过LTx=y解得x的值

end

改进平方根法函数程序如下:

function b=gaijinpinfanggenfa(A,b)

n=size(A);

n=n(1);

v=zeros(n,1);

for j=1:n

for i=1:j-1

v(i)=A(j,i)*A(i,i);

end

A(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);

A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j); end B=diag(A);

D=zeros(n);

for i=1:n

D(i,i)=B(i);

A(i,i)=1;

End

A=tril(A);

for j=1:n-1

b(j)=b(j)/A(j,j);

b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);

end

b(n)=b(n)/A(n,n); A=D*(A');

for j=n:-1:2

b(j)=b(j)/A(j,j);

b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);

end

b(1)=b(1)/A(1,1); disp('改进平方根法解得的解即为b');

end

调用函数解题:

clear;clc;

n=input('请输入矩阵维数:');b=zeros(n,1); A=zeros(n);

for i=1:n

for j=1:n

A(i,j)=1/(i+j-1); %LDL'分解 %得到L和D %前代法 %回代法

b(i)=b(i)+1/(i+j-1);

end

end %生成hilbert矩阵

[x,b]=pingfanggenfa(A,b)

b=gaijinpinfanggenfa(A,b)

运行结果:

请输入矩阵维数:40

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4

In qiujie at 10

Matlab自带解即为x

平方根法的解即为b

x =

1.6035

8.9685

0.8562

1.0195

0.9375

-50.2500

-3.0000

-16.0000

24.0000

-49.5000

-30.0000

39.0000

22.0000

-64.0000

-12.0000

2.0000

10.2500

-10.5000

-1.0000

-10.8750

83.0000

46.0000

-98.0000

12.0000

-69.0000

68.0000

21.0000

17.0000 -50.7188 -8.7500 -8.0000 112.0000 6.0000 -68.7500 22.0000 44.0000 -28.0000 8.0000 -44.0000 12.0000

b =

1.0e+007 *

0.0000 -0.0000 0.0001 -0.0004 -0.0014 0.0424 -0.2980

1.1419 -2.7335

4.2539 -4.3018

2.7733 -1.1989 0.5406 -0.3688 0.3285 -0.4438 0.4621 -0.2513 0.0565 0.0000 -0.0051 0.0071 -0.0027 -0.0031

0.0036

-0.0019

0.0009

0.0002

-0.0002

-0.0006

0.0004

0.0001

-0.0002

0.0001

0.0000

-0.0000

0.0000

-0.0000

-0.0000

改进平方根法解得的解即为b

b =

1.0e+024 *

0.0000

-0.0000

0.0001

-0.0012

0.0139

-0.0954

0.4208

-1.2101

2.0624

-1.0394

-3.3343

6.2567

-0.2463

-7.4594

2.8030

3.6990

0.7277

-1.7484

-0.4854

-3.6010

0.2532

5.1862

-2.1299

1.4410

0.8738

-4.5654

1.0422

4.0920

-2.7764

-2.2148

-0.8953

0.3665

4.8967

1.0416

0.1281

-4.3387

-1.1902

-2.8334

8.4610

-3.6008

实验二、利用QR分解解线性方程组:

利用QR分解解线性方程组Ax=b,其中 A=[16 4 8 4;

4 10 8 4;

8 8 12 10;

4 4 10 12];

b=[32 26 38 30];

求解程序如下:

定义house函数:

function [v,B]=house(x)

n=length(x);

y=norm(x,inf);

x=x/y;

Q=x(2:n)'*x(2:n);

v(1)=1;v(2:n)=x(2:n);

if n==1

Q=0;

B=0;

else

a=sqrt(x(1)^2+Q);

if x(1)<=0

v(1)=x(1)-a;

else

v(1)=-Q/(x(1)+a);

end

B=2*v(1)^2/(Q+v(1)^2);

v=v/v(1);

end

end

进行QR分解:

clear;clc;

A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];

b=[32 26 38 30];

b=b';

x=size(A);

m=x(1);

n=x(2);

d=zeros(n,1);

for j=1:n

[v,B]=house(A(j:m,j));

A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n); d(j)=B;

if j<m

A(j+1:m,j)=v(2:m-j+1);

end

end %QR分解 R=triu(A); %得到R D=A;

I=eye(m,n);

Q=I;

for i=1:n

D(i,i)=1;

end

H=tril(D);

M=H';

for i=1:n

N=I-d(i)*H(1:m,i)*M(i,1:m);

Q=Q*N;

end %得到Q

b=(Q')*b; %Q是正交阵 for j=n:-1:2

b(j)=b(j)/R(j,j);

b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);

end

b(1)=b(1)/R(1,1); %回带法

运行结果如下:

R =

18.7617 9.8072 15.7769 11.0864

0 9.9909 9.3358 7.5341

0 0 5.9945 9.8013

0 0 0 -0.5126

Q =

0.8528 -0.4368 -0.2297 -0.1709

0.2132 0.7916 -0.4594 -0.3417

0.4264 0.3822 0.2844 0.7689

0.2132 0.1911 0.8095 -0.5126

b=

1.00000000000000

1.00000000000001

0.999999999999988

1.00000000000001

实验三、Newton下山法解非线性方程组:

3x-cos(yz)-=0,

-81+sinz+1.06=0,

exp(-xy)+20z+=0;

要求满足数值解=满足或.

定义所求方程组的函数:Newtonfun.m

function F = Newtonfun(X)

F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;

F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;

F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;

End

向量求导:Xiangliangqiudao.m

function J=xiangliangqiudao()

syms x y z

X=[x,y,z];

F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3];

J=jacobian(F,[x y z]);

End

代值函数:Jacobi.m

function F=Jacobi(x)

F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];

End

方程组求解:

format long; %数据表示为双精度型 X1=[0,0,0]';

eps=10^(-8);

k=1;i=1;

X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);

while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps)

if norm(Newtonfun(X2),2)<norm(Newtonfun(X1),2) %判断先后两次迭代的大小 X1=X2;

B=inv(Jacobi(X2));

C=Newtonfun(X2);

X2=X2-B*C;

i=i+1;

else

v=1/(2^k); %引入下山因子

X1=X2;

B=inv(Jacobi(X2));

C=Newtonfun(X2);

X2=X2-v*B*C;

k=k+1;

end

end

j=i+k-1 %迭代次数

X=X2 %输出结果

运行结果如下: j =

5

X =

0.500000000000000 -0.000000000000000 -0.523598775598299

相关推荐