北邮matlab物理专业考试题!!

一、Matlab 基本命令，循环排序，函数，矩阵等问题 1、下列命令的意思：

det(A) — 求行列式 eig(A)—返回矩阵的特征值和特征向量 fliplr(A)—左右翻转矩阵A mean(A)—求A 的平均值

std(A)—求A 的均方差 angle — 求相位角

ifft — 快速傅里叶变换 sort(A)— 按升序排列

2、a =[3 2 5 6 1 4]。按下列"选择排序法"编程排序.

（1）找出全部N 个数据中最小的一个，把它和数列的第1个数字交换位置。

（2）找出剩下N-1个数据中最小的一个，把它和数列的第2个数字交换位置。

（3）..........

（4）一直做到只剩下一个数据为止。

clear;clc

a=[3 2 5 6 1 4];

for i=1:5

min=a(i);n=i;

for j=i+1:6

if a(j)

min=a(j);

n=j;

end

end

a(n)=a(i);a(i)=min;

end

a

3、分别用“函数文件”和“用函数句柄@表示函数”2种方法，编程求解：

(1)当3, 2.1x θ==，求函数22()cos ()f x x θ=+ 的值; (2)将计算结果存在cc.txt 文件中。 句柄函数

clear;clc

theta=2.1;

f=@(x)(cos(x^2))^2+theta;

A=f(3)

save cc.txt A -ascii

function fun=f(x,theta)

fun=(cos(x^2))^2+theta;

f(3,2.1)

4、给定1，2，...,7，编程计算：

（1）若任意顺序排列，一共有多少种不同的排列

a=[1,2,3,4,5,6,7];

A=perms(a);

length(A)

（2）若将其连成一个环，考虑旋转对称性，还剩下多少种排列？把它们都找出来存在文档p7.txt 文件中

A=perms(2:6);

[r,~]=size(A);

m=0;

for i=1:r-m

for j=i+1:r-m

if A(i,:)==fliplr(A(j,:))

A(j,:)=[];

m=m+1;

break ;

end

end

end

[r,~]=size(A);

I=zeros(r,1)+1;

A=[I,A];

save p7.txt A -ascii

（3）在Matlab 中调用p7.txt 文件，矩阵p7的某一行可写为a1，a2...a7，由此组成7阶方阵A 如下：

1111222233221111, n n n n n n n n n n a a a a a a a a a a a a A a a a a a a a a --------????

--????--=??????--??--???? n =8

该矩阵对称，行和为0，特征值的规律为1230n λλλλ=>≥≥≥ 。求出p7.txt 所有排列的λ2和λ7，并分别找出λ2最大值和λ7最小值所对应的排列；以及找出

λ7/λ2最大值和λ7/λ2最小值所对应的排列。

clear,clc

load p7.txt

[r,c]=size(p7);

for i=1:r

R=[p7(i,:),p7(i,:)];

A=zeros(c,2*c);

for j=1:c

A(j,[j j+1 j+2])=[R(c+j-1) -R(c+j-1)-R(j) R(j)];

B=A(:,2:c+1);

B(:,1)=A(:,c+2)+B(:,1);

B(:,c)=A(:,1)+B(:,c);

end

E=sort(eig(B),'descend');

E2(i)=E(2);E7(i)=E(7);

[r2max,r2maxid]=max(E2);

[r7min,r7minid]=min(E2);

[r72max,r72maxid]=max(E7./E2);

[r72min,r72minid]=min(E7./E2);

end

disp(['r2max=',num2str(r2max),'对应排列:',num2str(p7(r2maxid,:))])

disp(['r7min=',num2str(r7min),'对应排列:',num2str(p7(r7minid,:))])

disp(['r72max=',num2str(r72max),'对应排列:',num2str(p7(r72maxid,:))]) disp(['r72min=',num2str(r72min),'对应排列:',num2str(p7(r72minid,:))])

二、混沌映射

1、logistic映射：x(n+1) = a*x(n)*(1-x(n)) 。初值x(1)=0.3。

(1) 仅取n = [9900 10000]时（抛除n=9900之前的暂态），用subplot画出下面四个参数下的时间序列：a=2.5；a =3.2 ; a =3.8 ; a =4.0，并在图上分别标出横纵座标和图名。

clear;clc

a=[2.5,3.2,3.8,4];

str{1}='a=2.5μ??ì??ó3é?';

str{2}='a=3.2μ??ì??ó3é?';

str{3}='a=3.8μ??ì??ó3é?';

str{4}='a=4.0μ??ì??ó3é?';

for i=1:4

x=0.3;

A=[];

for n=1:9889

x=a(i)*x*(1-x);

end

for n=9900:10000

x=a(i)*x*(1-x);

A=[A;x];

end

subplot(2,2,i)

plot(9900:10000,A)

xlabel('t'); ylabel('x'); title(str{i})

end

(2) 当a =4时，先抛除暂态5000后，再用hist命令画出10000个x 值的分布图（称为序列A）. clear;clc

a=4;A=[];

x=0.3;

for n=1:5000

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

end

for n=5001:15000

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

A=[A;x];

end

hist(A,[0:0.0001:1])

(3) 求出序列A 的自相关系数；再随机取另一个初值，并产生另一个相同长度的序列（序列

B ），求A 和B 的互相关系数。

clear;clc

a=4;A=[];

x=0.3;

for n=1:5000

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

end

for n=5001:15000

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

A=[A;x];

end

a=4;B=[];

x=rand;

for n=1:5000

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

end

for n=5001:15000

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

B=[B;x];

end

R=cov(A)

R12=cov(A,B)

2、CML 同步

（1）方程：

， n 为迭代时间，i 为空间，i =1,..., L ，L =100。f(x)=1-ax 2,a=1.9,

1.0,7.021==γγ。初始条件为任意的[0 1]的随机数。周期边界条件:x 0n = x L n , x L +1n = x 1n

（2）方程：

为与上面参数一样的系统，如果初始条件不同（均为任意的随机数），验证这两系统是不能同步的。

clear;clc

c=1.9;a=0.7;b=0.1;L=100;

f=@(x)1-c*x.^2;

x=zeros(101,101);

y=zeros(101,101);

x(:,1)=rand(101,1);x(1,:)=rand(1,101);

y(:,1)=rand(101,1);y(1,:)=rand(1,101);

for n=1:100

for i=2:L

x(1,n)=x(L+1,n);

x(L+1,n)=x(1,n);

x(i,n+1)=(1-a-b)*f(x(i,n))+a*f(x(i-1,n))+b*f(x(i+1,n));

y(1,n)=y(L+1,n);

y(L+1,n)=y(1,n);

y(i,n+1)=(1-a-b)*f(y(i,n))+a*f(y(i-1,n))+b*f(y(i+1,n));

end

end

z=x-y;

[I,N]=meshgrid(1:101,1:101);

plot3(N,I,z)

plot(1:101,z(100,:))

（3）假定11

n n x y =，即用方程（1）驱动方程（2），这时方程（2）的边界条件为：

)()()1(,111111

-++-==L n L n L n n n y f y f y x y γγ。经过一段时间后[取10000]，两个时空混沌系统

会一一对应地同步。用时间n 作为横纵坐标，纵坐标为i y x i d i n i n /100)(+-=（i =1...100)。

a=1.9;b=0.7;c=0.1;

f=@(x) 1-a*x.^2;

L=100;N=15000;

x=rand(N,L+2);y=rand(N,L+2);

for n=1:N-1

y(n,1)=x(n,1);

y(n+1,L+1)=(1-b)*f(y(n,L+1))+b*f(y(n,L));

for i=2:L+1

x(n,1)=x(n,L+1);

x(n,L+2)=x(n,2);

x(n+1,i)=(1-b-c)*f(x(n,i))+b*f(x(n,i-1))+c*f(x(n,i+1));

y(n+1,i)=(1-b-c)*f(y(n,i))+b*f(y(n,i-1))+c*f(y(n,i+1));

end

end

d=zeros(size(x));

for i=2:L+1

d(:,i)=x(:,i)-y(:,i)+100/(i-1);

end

plot(1:N,d(1:N,:),'.')

（4）假设)()()1(),()()1(,1.01111111111-+-++-=+-===L n L n L n L n L n L n n n x f x f x y f y f y x y γγγγ，

此时两个系统经过一点时间后，在刚开始的一些空间格点会同步，但是到后面的空间格点就不再同步，请验证。

（5）假设

)()()1(),()()1(,]01[1111111111-+-++-=+-===L n L n L n L n L n L n n n x f x f x y f y f y x y γγγγ之间随机数，此时两个系统经过一段时间后，能够实现全部同步，请验证。

（6）逐渐增加2γ，当2γ大于何值时，问题（3）中的两个系统就不再同步了？

4)

a=1.9;r1=0.7;r2=0.1;

f=@(x) 1-a*x.^2;

L=100;N=1e+4;

x=rand(N,L+2);y=rand(N,L+2);

for n=1:N-1

x(n,1)=0.1;y(n,1)=0.1;

x(n+1,L+1)=(1-r1)*f(x(n,L+1))+r1*f(x(n,L));

y(n+1,L+1)=(1-r1)*f(y(n,L+1))+r1*f(y(n,L));

for i=2:L+1

x(n+1,i)=(1-r1-r2)*f(x(n,i))+r1*f(x(n,i-1))+r2*f(x(n,i+1)); y(n+1,i)=(1-r1-r2)*f(y(n,i))+r1*f(y(n,i-1))+r2*f(y(n,i+1)); end

end

d=zeros(size(x));

for i=2:L+1

d(:,i)=x(:,i)-y(:,i)+100/(i-1);

end

plot(9900:1e+4,d(9900:1e+4,30:100))

5)

a=1.9;r1=0.7;r2=0.1;

f=@(x) 1-a*x.^2;

L=100;N=1e+4;

x=rand(N,L+2);y=rand(N,L+2);

x0=rand();y0=rand();

for n=1:N-1

x(n,1)=x0;y(n,1)=y0;

x(n+1,L+1)=(1-r1)*f(x(n,L+1))+r1*f(x(n,L));

y(n+1,L+1)=(1-r1)*f(y(n,L+1))+r1*f(y(n,L));

for i=2:L+1

x(n+1,i)=(1-r1-r2)*f(x(n,i))+r1*f(x(n,i-1))+r2*f(x(n,i+1));

y(n+1,i)=(1-r1-r2)*f(y(n,i))+r1*f(y(n,i-1))+r2*f(y(n,i+1));

end

end

d=zeros(size(x));

for i=2:L+1

d(:,i)=x(:,i)-y(:,i)+100/(i-1);

end

plot(9900:1e+4,d(9900:1e+4,:))

6)

a=1.9;r1=0.7;r2=0.1;

f=@(x) 1-a*x.^2;

L=100;N=1e+4;k=0;

while 1

r2=r2+0.1;

x=rand(N,L+2);y=rand(N,L+2);

for n=1:N-1

y(n,1)=x(n,1);

y(n+1,L+1)=(1-r1)*f(y(n,L+1))+r1*f(y(n,L));

for i=2:L+1

x(n+1,i)=(1-r1-r2)*f(x(n,i))+r1*f(x(n,i-1))+r2*f(x(n,i+1));

y(n+1,i)=(1-r1-r2)*f(y(n,i))+r1*f(y(n,i-1))+r2*f(y(n,i+1));

end

end

d=zeros(size(x));

for i=2:L+1

d(:,i)=x(:,i)-y(:,i)+100/(i-1);

end

if std(d(9900:1e+4,10))>1e+2

cr=r2;

plot(9900:1e+4,d(9900:1e+4,2:L+1))

disp(cr);

break;

end

disp(['μ±?°r2=',num2str(r2)]);

end

三、庞加莱截面图是采用频闪法得到的一组“点”线。对于正、余弦等周期运动，频闪法可以很快得到她的周期。例如单摆运动，在摆球过零时用光电们记录，就是一种频闪记录方法，我们在玻尔共振实验中用的频闪也是这样判断周期和相位差的。在混沌系统中，因为有很多的周期（或周期为无限长），这是频闪就会得到一系列的点；但周期1态只应观察到一个值，

周期2就该观察到2个值，以此类推。对于chua 系统，具体做法是：判断x 变量为最大时，记录该时刻的y 值，来画不同参数时的分岔图，这个图是一种庞加莱截面图的集结。如果x y 均作周期1运动，只能观测到一个点（每次记录的y 值相同）；周期2时，就可记录到2个y 值，等等。

根据以上描述，编写一段程序求得chua 系统的分岔图。chua 电路方程：

(())

dx y x f x dt dy x y z dt dz y dt αβ?=--???=-+???=-??

, 其中，()0.5()(11)f x bx a b x x =+-+--, )25,14(,0.9,7/5,7/8∈=-=-=βαb a ，初值：0000,0.01,0.01x y z ===,或自选.

（1） 用自编RK4法求解Chua's 电路，分别画出6.15,16,18,24=β,15.2,14的相图(x-y) clear,clc

a=-8/7;b=-5/7;alf=9;beta=[24,18,16,15.6,15.2,14]

f=@(y)b*y+0.5*(a-b)*(abs(y+1)-abs(y-1));

for k=1:6

ydot=@(t,y)[alf*(y(2)-y(1)-f(y(1))),y(1)-y(2)+y(3),-beta(k)*y(2)]; y0=[0 0.01 0.01];

h=0.01;A=y0;t=0;

for i=1:140/h

k1=h*ydot(t,y0);

k2=h*ydot(t+h/2,y0+k1/2);

k3=h*ydot(t+h/2,y0+k2/2);

k4=h*ydot(t+h,y0+k3);

yy=y0+(k1+2*k2+2*k3+k4)/6;

A=[A;yy];

y0=yy;

end

subplot(2,3,k)

plot(A(10000:14000,1),A(10000:14000,2))

end

β=，用Matlab中的ode45方法，时间步长取0.01，计算chua's （2）当14.0:0.1:25.0

-），横坐标为β。

方程的分岔图（yβ

function ydot=chua(t,y)

global alf beta a b

f=@(y)b*y+0.5*(a-b)*(abs(y+1)-abs(y-1));

ydot=[alf*(y(2)-y(1)-f(y(1)));

y(1)-y(2)+y(3);

-beta*y(2)];

end

clear,clc

global alf beta a b

a=-8/7;b=-5/7;alf=9;

l=14:0.1:25;h=0.01;

y0=[0,0.01,0.01];

L=length(l);

for i=1:L

beta=l(i);

[T,Y]=ode45(@chua,0:0.01:100,y0);

A=[0 0];

for j=50/h:length(T)-1

if Y(j,1)>Y(j-1,1)&&Y(j,1)>Y(j+1,1)

A=[A;beta,Y(j,2)];

end

end

plot(A(2:end,1),A(2:end,2),'.');

axis([14 25 -.6 .6])

hold on

end

四、如下图所示：矩形的盒子被一档板等体积的分成两部分，气体分子现在左边的区域由一点（共800个分子）扩散开来，达到稳定后（设t=10000）把档板抽去，气体分子扩散到右边区域。请编程模拟这一过程，在2维平面画出t =1，10，50，100，500，1000，5000，9500，10000和t =10001, 10010, 10050, 10100, 10500, 11000, 15000,19500,20000时的分子扩散状态图，并计算熵随时间的变化。

clear all;clc

ns=2e+4;np=800;

rand('state',sum(100*clock));

ss=2;

x=zeros(ns,np);y=zeros(ns,np);

for s=2:ns

for p=1:np

if s<1e+4

bound=[-100 100;-100 100];

else

bound=[-100 300;-100 100];

end

theta=rand*2*pi;

x(s,p)=x(s-1,p)+ss*cos(theta);

if x(s,p)bound(1,2)

x(s,p)=x(s-1,p)-2*ss*cos(theta);

end

y(s,p)=y(s-1,p)+ss*sin(theta);

if y(s,p)bound(2,2)

y(s,p)=y(s-1,p)-2*ss*sin(theta);

end

end

end

T1=[1 10 50 100 500 1000 5000 9500 10000];

T2=[10001 10010 10050 10100 10500 11000 15000 19500 20000]; I=length(T1);

for i=1:I

subplot(3,3,i)

plot(x(T1(i),:),y(T1(i),:),'o')

axis([-100 100 -100 100])

xlabel('x');ylabel('y');

title(['t=',num2str(T1(i))]);

end

figure

for i=1:I

subplot(3,3,i)

plot(x(T2(i),:),y(T2(i),:),'o')

axis([-100 300 -100 100])

xlabel('x');ylabel('y');

title(['t=',num2str(T2(i))]);

end

%entropy

E=[0 0];

for s=1:ns

if s<1e+4

bound=[-100 100;-100 100];

else

bound=[-100 300;-100 100];

end

hx=histc(x(s,:),bound(1,1):20:bound(1,2));

hy=histc(y(s,:),bound(2,1):20:bound(2,2));

for i=1:(bound(1,2)-bound(1,1))/20

for j=1:(bound(2,2)-bound(2,1))/20

P(i,j)=hx(i)*hy(j)/np^2;

if P(i,j)==0

E(s,1)=s;

E(s,2)=E(s,2);

else

E(s,1)=s;

E(s,2)=E(s,2)-1.443*(P(i,j)*log(P(i,j)));

end

end

end

end

plot(E(:,1),E(:,2))

五、将1，2，...,10排列成镜像对称的三角型(1 3 5 7 9 10 8 6 4 2)，画出排列图；（1）如果将三角型的任意1对数对调，一共有多少种不同的排列，请列举；

A=[1 3 5 7 9 10 8 6 4 2];

C=A;

plot(A)

n=length(A);

for i=1:n/2

for j=n/2+1:n

B=A;

temp=B(j);

B(j)=B(i);

B(i)=temp;

C=[C;B];

end

end

C

（2）如果将三角型的任意2对数对调，一共能形成多少种新的不同的排列，请列举。 六、自激振荡模型：121 1.9910.0010.001n n n n n n n

u v u v v u σξ++=+++=--????? ，其中ξ是高斯白噪声。取初始值000(,,)(1, 1.995,u v ξ=--小随机数）

。当σ分别取为：0.0015，0.015，0.09，并在n=[1000,9000] 范围内加噪声(其他时间不加噪声），分别画出n=1,12000时的u n -n ， v n -n ； u n - v n 图。

clear,clc

sigma=[0.0015 0.015 0.09];

for i=1:3

A=[];B=[];

u=-1;v=-1.995;

for n=1:999

A=[A;u];

B=[B,v];

u=1.99/(1+u.^2)+v;

v=v-0.001*u-0.001;

end

for n=1000:9000

e=randn;

A=[A;u];

B=[B,v];

u=1.99/(1+u.^2)+v+sigma(i)*e;

v=v-0.001*u-0.001;

end

for n=9001:12000

A=[A;u];

B=[B,v];

u=1.99/(1+u.^2)+v;

v=v-0.001*u-0.001;

end

subplot(3,2,2*i-1)

plot(1:12000,A)

hold on

plot(1:12000,B)

xlabel('n');ylabel('un,vn');title(['\sigma =',num2str(sigma(i))]) axis([-1500 12000 -2.5 0.5])

subplot(3,2,2*i)

plot(A,B)

xlabel('un');ylabel('vn');title(['\sigma =',num2str(sigma(i))]) axis([-2 0.5 -2.16 -1.96])

end