北邮计算物理参考代码.docx
《北邮计算物理参考代码.docx》由会员分享，可在线阅读，更多相关《北邮计算物理参考代码.docx（43页珍藏版）》请在冰豆网上搜索。
北邮计算物理参考代码
测试第一题：
二阶RK2函数：
functionout2=RK2（x,time,tau,driverRK3）
F1=feval（driverRK3,x,time）;
t_full=time+tau/2;
xtemp=x+tau*F1/2;
F2=feval（driverRK3,xtemp,t_full）;
out2=x+tau*F2;
end
四阶RK4函数：
functionout4=RK4（x,time,tau,driverRK3）
F1=feval（driverRK3,x,time）;
t_half=time+tau/2;
xtemp=x+tau*F1/2;
F2=feval（driverRK3,xtemp,t_half）;
xtemp=x+tau*F2/2;
F3=feval（driverRK3,xtemp,t_half）;
xtemp=x+tau*F3;
t_full=time+tau;
F4=feval（driverRK3,xtemp,t_full）;
out4=x+tau/6.*（F1+2*F2+2*F3+F4）;
end
自定义gravrk函数：
functiondrive=gravrk（y,t）
globalFmgama1gama2kKw;
F=1;
m=1.55;
gama1=0.04;
gama2=0.00001;
k=0.62;
K=0.1395;
w=2*pi;
xz=[y
（1）y（3）];
vz=[y
（2）y（4）];
accel1=（F*exp（i*w*t）-k*y
（1）+K*（y（3）-y
（1））-gama1*y
（2）/m）/m;
accel2=（K*（y
（1）-y（3））-k*y（3）-gama2*y（4）/m）/m;
drive=[vz
（1）accel1vz
（2）accel2];
end
二阶RK算法主函数：
clearall
globalFmgama1gama2kKw;
F=1;
m=1.55;
gama1=0.04;
gama2=0.00001;
k=0.62;
K=0.1395;
w=2*pi;
tau=0.01;
nStep=100*（1/tau）;
time=0;
state=[0000];
fori=1:
nStep
state=RK2（state,time,tau,'gravrk2'）;
x=[state
（1）state（3）];
v=[state
（2）state（4）];
time=time+tau;
Sx（i）=state
（1）;
Sy（i）=state（3）;
Sv1（i）=state
（2）;
Sv2（i）=state（4）;
timeS（i）=time;
end
plot（timeS,Sx,'-b',timeS,Sy,'--r',timeS,Sv1,':
m',timeS,Sv2,'-.k'）
四阶RK算法主函数：
clearall
globalFmgama1gama2kKw;
F=1;
m=1.55;
gama1=0.04;
gama2=0.00001;
k=0.62;
K=0.1395;
w=2*pi;
tau=0.01;
nStep=100*（1/tau）;
time=0;
state=[0000];
fori=1:
nStep
state=RK4（state,time,tau,'gravrk2'）;
x=[state
（1）state（3）];
v=[state
（2）state（4）];
time=time+tau;
Sx（i）=state
（1）;
Sy（i）=state（3）;
Sv1（i）=state
（2）;
Sv2（i）=state（4）;
timeS（i）=time;
end
plot（timeS,Sx,'-b',timeS,Sy,'--r',timeS,Sv1,':
m',timeS,Sv2,'-.k'）
利用ode45主函数：
clearall
globalFmgama1gama2kKw;
F=1;
m=1.55;
gama1=0.04;
gama2=0.00001;
k=0.62;
K=0.1395;
w=2*pi;
thzz_fun=@（t,y）[y
（2）/m;（F*exp（i*w*t）-k*y
（1）+K*（y（3）-y
（1））-gama1*y
（2）/m）/m;y（4）/m;（K*（y
（1）-y（3））-k*y（3）-gama2*y（4）/m）/m;];
[t,u2]=ode45（thzz_fun,[0,200],[0,0,0,0]）;
plot（t,u2（:
1）,'-r',t,u2（:
2）,'--b',t,u2（:
3）,'-.g',t,u2（:
4）,':
p'）;
幅度计算功能函数：
functionfupin=fupinq（w）
globalFmgama1gama2kK;
F=1;
m=1.55;
gama1=0.04;
gama2=0.00001;
k=0.62;
K=0.1395;
thzz_fun=@（t,y）[y
（2）/m;（F*exp（i*w*t）-k*y
（1）+K*（y（3）-y
（1））-gama1*y
（2）/m）/m;y（4）/m;（K*（y
（1）-y（3））-k*y（3）-gama2*y（4）/m）/m;];
[t,u2]=ode45（thzz_fun,[0,50],[0,0,0,0]）;
fupin=max（2*u2（3））;
end
作幅频曲线主函数：
clearall
fori=0:
400
theta=pi*i/10;
fu=fupinq（theta）;
plot（theta,fu,'MarkerSize',3）;
holdon
end
动图主函数：
clearall
fortheta1=1:
800
fu（theta1）=fupinq（theta1/1000）;
theta（theta1）=theta1;
end
comet（theta,fu）;
、
测试第二题
functionout4=RK4（x,time,tau,driverRK3,beta）
F1=feval（driverRK3,x,time,beta）;
t_half=time+tau/2;
xtemp=x+tau*F1/2;
F2=feval（driverRK3,xtemp,t_half,beta）;
xtemp=x+tau*F2/2;
F3=feval（driverRK3,xtemp,t_half,beta）;
xtemp=x+tau*F3;
t_full=time+tau;
F4=feval（driverRK3,xtemp,t_full,beta）;
out4=x+tau/6.*（F1+2*F2+2*F3+F4）;
end
functiondrive=gravrk（y,t,beta）
globalabaerfa;
a=-8/7;b=-5/7;aerfa=9.0;
fun31=@（x）（b*x+0.5*（a-b）*（abs（x+1）-abs（x-1）））;
accel1=aerfa*（y
（2）-y
（1）-fun31（y
（1）））;
accel2=y
（1）-y
（2）+y（3）;
accel3=-beta*y
（2）;
drive=[accel1accel2accel3];
end
functionmain_3_RK4（beta）
globalabaerfa;
a=-8/7;b=-5/7;aerfa=9.0;
tau=0.01;
nStep=20*（1/tau）;
time=0;
state=[00.010.01];
fori=1:
nStep
state=RK4（state,time,tau,'gravrk',beta）;
time=time+tau;
Sx（i）=state
（1）;
Sy（i）=state
（2）;
Sz（i）=state（3）;
timeS（i）=time;
end
plot（Sx,Sy）;
end
functionmain_3_1_z
a=[24181615.615.214];
fori=1:
6
subplot（2,3,i）
main_3_RK4（a（i））
end
end
functionydot=main_3_ode45（beta）
globalabaerfa;
a=-8/7;b=-5/7;aerfa=9.0;
fun31=@（x）（b*x+0.5*（a-b）*（abs（x+1）-abs（x-1）））;
fun_3_1=@（t,y）[aerfa*（y
（2）-y
（1）-fun31（y
（1）））;y
（1）-y
（2）+y（3）;-beta*y
（2）];
[t,u]=ode45（fun_3_1,[0:
0.01:
10],[0,0.01,0.01]）;
plot（beta,u（500:
10:
1000,2）,'r*','MarkerSize',1）;
holdon
xlabel（'\beta'）;ylabel（'y'）;
end
functionmain_3_1
fora=14:
0.1:
25;
main_3_ode45（a）
end
end
测试第三题：
clearall
e=0.8;
x=zeros（5000,51）;
x（1,:
）=rand（1,51）;
forn=1:
5000
x（n,1）=0.2;
end
form=1:
5000
fori=2:
51
x（m+1,i）=（1-e）*4*x（m,i）*（1-x（m,i））+e*4*x（m,i-1）*（1-x（m,i-1））;
end
end
subplot（2,4,1）
plot（1:
50,x（1,2:
51）,'.'）;
subplot（2,4,2）
plot（1:
50,x（50,2:
51）,'.'）;
subplot（2,4,3）
plot（1:
50,x（100,2:
51）,'.'）;
subplot（2,4,4）
plot（1:
50,x（500,2:
51）,'.'）;
subplot（2,4,5）
plot（1:
50,x（1000,2:
51）,'.'）;
subplot（2,4,6）
plot（1:
50,x（5000,2:
51）,'.'）;
subplot（2,4,7）
pcolor（0:
49,0:
100,x（10:
110,2:
51））;
作业题第一题
functiondjdb2
theta=[1/180*pi5/180*pi10/180*pi];
T=[];
ops=odeset（'Events',@events）;
fork=1:
3
[t,u]=ode45（@f,[0,40],[theta（k）,0],ops）;%ÊýÖµ½â
T=[T,2*t（end）];
end
k2=sin（theta./2）.^2;%½âÎö½â
[K,E]=ellipke（k2）;
error=4*K-T;
plot（theta,error,'.-'）
title（'Îó²î'）%ͼÖпɼûÎó²îÔÚ10^-4Á¿¼¶
end
functionydot=f（t,y）
ydot=[y
（2）;-sin（y
（1））];
end
function[value,isterminal,direction]=events（t,y）
value=（y
（2））;
isterminal=1;
direction=1;
end
作业题第二题
functionxiangtu
beta1=0.2;
beta2=1.2;
f1=@（t,y）[y
（2）;-2*beta1*y
（2）-sin（y
（1））];
f2=@（t,y）[y
（2）;-2*beta2*y
（2）-sin（y
（1））];
[t1,w1]=ode45（f1,[0,10*pi],[pi/10;0]）;
[t2,w2]=ode45（f2,[0,10*pi],[pi/10;0]）;
subplot（2,1,1）
plot（t1,w1（:
1）,'-.k',t2,w2（:
1）,'--r'）;
xlabel（'t'）,ylabel（'\theta'）;
title（'½ÇλÒÆ'）
subplot（2,1,2）
plot（w1（:
1）,w1（:
2）,'b'）;
holdon
plot（w2（:
1）,w2（:
2）,'g:
'）
holdoff
title（'Ïàͼ'）
end
functionydot=fun_2（t,y,beta）
ydot=[y
（2）;
-2*beta*y
（2）-sin（y
（1））];
end
作业题第三题
functiondbyd_3
u=2/3;a=0.5;zq=3*pi;f=[1.07261.088291.095];
fori=1:
3
ydot=@（t,y）[y
（2）;
-sin（y
（1））-a*y
（2）+f（i）*cos（u*t）];
[T,Y]=ode45（ydot,[0:
zq/100:
20*zq],[-0.8;2]）;
subplot（3,1,i）
plot（T,Y（:
1）,'r'）
end
end
作业题第五题
functionmain_5
k=0.1;
f=@（t,y）[y
（2）;
-k*y
（2）+y
（1）-y
（1）^3];
[t,y1]=ode45（f,[0,200],[0,0.8]）;
[t,y2]=ode45（f,[0,200],[0,-0.8]）;
plot（y1（:
1）,y1（:
2）,'r'）
holdon
plot（y2（:
1）,y2（:
2）,'k'）
title（'ÓÐ×èÄáÎÞÇý¶¯Ïàͼ'）
end
作业题第六题
functionmain_6
globaluwv
u=0.3;w=0.44;T=2*pi/w;v=0;
[t,y]=ode45（@fun_6,[0:
T/1000:
500*T],[5,4]）;
Y=fft（cos（w*t））;
Y
（1）=[];n=length（Y）;m=fix（n/2）;
power=abs（Y（1:
m））.^2/n^2;
dt=t
（2）-t
（1）;
freq=1/dt*（1:
n/2）./n;
plot（freq,power）
axis（[00.400.05]）
holdon
end
functionmain_6_2
globaluwv
u=0.3;
w=0.44;T=2*pi/w;v=0;
[t,y]=ode45（@fun_6,[0:
T/1000:
500*T],[5,4]）;
comet（y（:
1）,y（:
2））
end
functionydot=fun_6（t,y）
globaluwv
x0=1;w0=1;
ydot=[y
（2）;
u*（x0^2-y
（1）^2）*y
（2）-y
（1）*w0^2-v*cos（w*t）];
end
作业题第七题
functionmain_7
[t,y]=ode45（@fun_7,[0,31],[0.20,0.1,eps]）;
subplot（2,2,1）
plot3（y（:
1）,y（:
2）,y（:
3）,'r',y（1,1）,y（1,2）,y（1,3）,'k*'）
xlabel（'x'）,ylabel（'y'）,zlabel（'z'）
view（3）
subplot（2,2,2）
plot（y（:
1）,y（:
2）,'y',y（1,1）,y（1,2）,'k*'）
xlabel（'x'）,ylabel（'y'）
subplot（2,2,3）
plot（y（:
1）,y（:
3）,'b',y（1,1）,y（1,3）,'k*'）
xlabel（'x'）,ylabel（'z'）
subplot（2,2,4）
plot（y（:
2）,y（:
3）,'g',y（1,2）,y（1,3）,'k*'）
xlabel（'y'）,ylabel（'z'）;
end
functionydot=fun_7（t,y）
r=225;
ydot=[-10*y
（1）+10*y
（2）;
r*y
（1）-y
（1）*y（3）-y
（2）;
y
（1）*y
（2）-8/3*y（3）];
end
作业题第八题
functionmain_8
T=2*pi;
[t,y]=ode45（@fun_8,[0,31],[1,0,0]）;
plot（y（:
1）,y（:
3）,'r',y（1,1）,y（1,3）,'k*'）
xlabel（'x'）,ylabel（'z'）
end
functionydot=fun_8（t,y）
r=25;
ydot=[-10*y
（1）+10*y
（2）;
r*y
（1）-y
（1）*y（3）-y
（2）;
y
（1）*y
（2）-8/3*y（3）];
end
作业题第九题
functionmain_9_1
N=1000;X=zeros（N,1）;Y=X;
n=0.3;u=1.4;
X
（1）=rand
（1）;Y
（1）=rand
（1）;
fori=1:
N
X（i+1）=Y（i）+1-u*X（i）^2;
Y（i+1）=n*X（i）;
end
plot（X,Y,'r'）
end
functionmain_9_2
N=500;X=zeros（N,1）;Y=X;
n=0.3;u=linspace（1.05,1.082,50）;
xx=[];U=[];
forii=1:
50
uu=u（ii）;
U=[U,uu];
X
（1）=rand
（1）;Y
（1）=rand
（1）;
fori=1:
N
X（i+1）=Y（i）+1-uu*X（i）^2;
Y（i+1）=n*X（i）;
end
xx=[xx;X'];
end
plot（U,xx（:
100:
500）,'.r'）
title（'x-u¡¤?
?
¨ª?
?
'）
end
作业题第十题
functionmain_10
a=0.32;b=3;c=6;
f=@（t,y）[-y
（2）-y（3）;
y
（1）+a*y
（2）;
b+y
（1）*y（3）-c*y（3）];
[t,y]=ode45（f,[0:
0.01:
100],[000]）;
subplot（2,2,1）
plot3（y（:
1）,y（:
2）,y（:
3）,'r'）;
subplot（2,2,2）
plot（y（:
1）,y（:
2）,'b'）
subplot（2,2,3）
plot（y（:
1）,y（:
3）,'g'）
subplot（2,2,4）
plot（y（:
2）,y（:
3）,'y'）
end
functionmain_10_2
a=[0.30.350.3750.3860.39090.3980.40.411];b=2;c=4;
X=[];Y=[];Z=[];
forii=1:
8
aa=a（ii）;
f=@（t,y）[-y
（2）-y（3）;
y
（1）+aa*y
（2）;
b+y
（1）*y（3）-c*y（3）];
[t,y]=ode45（f,[0:
0.01:
100],[000]）;
X=[X;y（:
1）'];Y=[Y;y（:
2）'];Z=[Z;y（:
3）'];
end
forj=1:
4
figure
（1）
subplot（2,2,j）
plot3（X（j,5000:
end）,Y（j,5000:
end）,Z（j,5000:
end））
end
forj=5:
8
figure
（2）
subplot（2,2,j-4）
plot3（X（j,5000:
end）,Y（j,5000:
end）,Z（j,5000:
end））
end
end
作业题第十一题
functionout4=RK4（x,time,tau,driverRK3）
F1=feval（driverRK3,x,time）;
t_half=time+tau/2;
xtemp=x+tau*F1/2;
F2=feval（driverRK3,xtemp,t_half）;
xtemp=x+tau*F2/2;
F3=feval（driverRK3,xtemp,t_half）;
xtemp=x+tau*F3;
t_full=time+tau;
F4=feval（driverRK3,xtemp,t_full）;
out4=x+tau/6.*（F1+2*F2+2*F3+F4）;
end
functiondrive=gravrk（y,t）
globalabc;
a=10;b=8/3;c=1;
accel1=-a*y
（1）+a*y
（2）;
accel2=b*y
（1）-y
（2）-y
（1）*y（3）;
accel3=y
（1）*y
（2）-c*y（3）;
drive=[accel1accel2accel3];
end
functionmain_3_RK4
globalabc;
a=10;b=8/3;c=1;
tau=0.01;
nStep=20*（1/tau）;
time=0;
state=[00.081];
fori=1:
nStep
state=RK4（sta
升级会员 