北邮计算物理参考代码2.docx
《北邮计算物理参考代码2.docx》由会员分享，可在线阅读，更多相关《北邮计算物理参考代码2.docx（23页珍藏版）》请在冰点文库上搜索。
北邮计算物理参考代码2
第五章解常微分方程
1.
clearall
e=0.2;
a=0.9;
n=[1,100,500,1000,2000,5000];
x=zeros（10000,101）;
x（1,:
）=rand（1,101）;
form=1:
10000
fori=2:
101
ifi~=101
x（m+1,i）=（1-e）*（1-a*x（m,i）^2）+e*（2-a*x（m,i-1）^2+a*x（m,i+1）^2）/2;
else
x（m+1,i）=（1-e）*（1-a*x（m,i）^2）+e*（2-a*x（m,i-1）^2+a*x（m,1）^2）/2;
end
end
end
plot（1:
100,x（n
（1）,2:
101））;
plot（1:
100,x（n
（2）,2:
101））;
plot（1:
100,x（n（3）,2:
101））;
plot（1:
100,x（n（4）,2:
101））;
plot（1:
100,x（n（5）,2:
101））;
plot（1:
100,x（n（6）,2:
101））;
2.
clear
rm=10;
r=2:
0.1:
rm;
theta=（0:
360）*pi/180;
[TH,R]=meshgrid（theta,r）;
[X,Y]=pol2cart（TH,R）;
U=2*X./R.^3;
mesh（X,Y,U）,holdon;
plot3（[-2;2],[0;0],[0;0],'-o','LineWidth',2）,boxon;
3.
k=pi;
v=1;
A1=-1;
A2=1;
x=0:
0.01:
4;
t=0:
0.01:
10;
for i=1:
length（t）-1
y1=A1*sin（k*（x-v*t（i）））;
y2=A2*sin（k*（x+v*t（i）））;
y=y1+y2;
plot（x,y,'r'）,axis（[-1,5,-2.2,2.2]）;
grid on;
xlabel（'驻波'）;
pause（0.01）;
end
4.
clearall
e=0.3;
a=1.402;
x=zeros（20000,601）;
x（1,:
）=rand（1,601）;
x（:
1）=（sqrt（4*a+1）-1）/（2*a）;
form=1:
20000
fori=2:
601
x（m+1,i）=（1-e）*（1-4*x（m,i）*（1-x（m,i）））+e*4*x（m,i-1）*（1-x（m,i-1））;
end
end
x1=x（10000:
10031,:
）;
surf（x1）
P266-T10.
（1）
Y=dsolve（'Dy=-1000*（y-sin（t））+cos（t）','y（0）=1','t'）
t=0:
0.01:
1;
yy=subs（Y,t）;
plot（t,yy）;
Y=
exp（-1000*t）+sin（t）
（2）（3）（4）
F=@（t,y）[-1000*（y-sin（t））+cos（t）];
[x1,y1]=ode23（F,[0,1],1）;
disp（'ode23执行步数：
'）
length（y1）
[x2,y2]=ode23s（F,[0,1],1）;
disp（'ode23s执行步数：
'）
length（y2）
plot（x1,y1,'x',x2,y2,'o'）
ode23执行步数：
ans=
417
ode23s执行步数：
ans=
58
（5）
变化缓慢处ode23s步长更大。
（6）
变化剧烈处ode23步长更大
P266T11
（1）
clearall
r0=300;f0=150;a0=0.01;a=zeros（1,100）;n=1;
Timep_fun=@（t,y）[2*y
（1）-a0*y
（1）*y
（2）;-y
（2）+a0*y
（1）*y
（2）];
[ty]=ode45（Timep_fun,[0,15],[r0,f0]）;
plot（t,y（:
1）,'b',t,y（:
2）,'g'）;
holdon
z1=2*y（:
1）-a0*y（:
1）.*y（:
2）;
z2=-y（:
2）+a0*y（:
1）.*y（:
2）;
plot（t,z1,'r',t,z2,'y'）;
fori=2:
length（z1）-1
if（z1（i）>z1（i-1）&z1（i）>z1（i+1））
a（n）=t（i）;
n=n+1;
end
end
a
（2）-a
（1）
ans=
4.9738
clearall
r0=15;f0=22;a0=0.01;a=zeros（1,100）;n=1;
Timep_fun=@（t,y）[2*y
（1）-a0*y
（1）*y
（2）;-y
（2）+a0*y
（1）*y
（2）];
[ty]=ode45（Timep_fun,[0,15],[r0,f0]）;
plot（t,y（:
1）,'b',t,y（:
2）,'g'）;
holdon
z1=2*y（:
1）-a0*y（:
1）.*y（:
2）;
z2=-y（:
2）+a0*y（:
1）.*y（:
2）;
plot（t,z1,'r',t,z2,'y'）;
fori=2:
length（z1）-1
if（z1（i）>z1（i-1）&z1（i）>z1（i+1））
a（n）=t（i）;
n=n+1;
end
end
a
（2）-a
（1）
ans=
6.5967
clearall
r0=102;f0=198;a0=0.01;a=zeros（1,100）;n=1;
Timep_fun=@（t,y）[2*y
（1）-a0*y
（1）*y
（2）;-y
（2）+a0*y
（1）*y
（2）];
[ty]=ode45（Timep_fun,[0,15],[r0,f0]）;
plot（t,y（:
1）,'b',t,y（:
2）,'g'）;
holdon
z1=2*y（:
1）-a0*y（:
1）.*y（:
2）;
z2=-y（:
2）+a0*y（:
1）.*y（:
2）;
plot（t,z1,'r',t,z2,'y'）;
fori=2:
length（z1）-1
if（z1（i）>z1（i-1）&z1（i）>z1（i+1））
a（n）=t（i）;
n=n+1;
end
end
a
（2）-a
（1）
ans=
4.4671
clearall
r0=100;f0=50;a0=0.01;a=zeros（1,100）;n=1;
Timep_fun=@（t,y）[-y
（2）;2*y
（1）];
[ty]=ode45（Timep_fun,[0,15],[r0,f0]）;
plot（t,y（:
1）,'b',t,y（:
2）,'g'）;
holdon
z1=2*y（:
1）-a0*y（:
1）.*y（:
2）;
z2=-y（:
2）+a0*y（:
1）.*y（:
2）;
plot（t,z1,'r',t,z2,'y'）;
fori=2:
length（z1）-1
if（z1（i）>z1（i-1）&z1（i）>z1（i+1））
a（n）=t（i）;
n=n+1;
end
end
disp（'周期为：
'）
T=a
（2）-a
（1）
周期为：
T=
4.3781
pendulum
（1）
disp（'所求解的方程：
'）
disp（'Domega（t）/Dt=（-g/length）*theta（t）'）
disp（'Dtheta（t）/Dt=omrga（t）'）
所求解的方程：
Domega（t）/Dt=（-g/length）*theta（t）
Dtheta（t）/Dt=omrga（t）
clearall
[x,y]=dsolve（'Domega=（-9.8）*theta','Dtheta=omega','omega（0）=0,theta（0）=0.2','t'）
x=
-（7*5^（1/2）*sin（（7*5^（1/2）*t）/5））/25
y=
cos（（7*5^（1/2）*t）/5）/5
clearall
fun=@（t,y）[-9.8*y
（2）;y
（1）];
[ty]=ode45（fun,[0,10],[0,0.2]）;
y1=y（:
1）.^2/2;
y2=（1-cos（y（:
2）））*9.8;
plot（t,y1,'r',t,y2,'b',t,y1+y2,'g'）
行星轨迹
程序代码：
main函数：
clearall;
h=1;
r0=1;v0=h*（2*pi*r0）;
r=[r0,0];v=[0,v0];
state=[r
（1）,r
（2）,v
（1）,v
（2）];
GM=4*pi^2;
mass=1.;
tau=0.001;
nStep=6*（1/tau）;
time=0;
state=rk4（state,time,tau,'gravrk',GM）;
r=[state
（1）state
（2）];
v=[state（3）state（4）];
time=time+tau;
fori=1:
nStep
state=rk4（state,time,tau,'gravrk',GM）;
r=[state
（1）state
（2）];
v=[state（3）state（4）];
time=time+tau;
Sx（i）=r
（1）;
Sy（i）=r
（2）;
timeS（i）=time;
kinetic（i）=.5*mass*（v
（1）^2+v
（2）^2）;
portential（i）=-GM*mass/（r
（1）^2+r
（2）^2）^.5;
Lz（i）=mass*（r
（1）*v
（2）-r
（2）*v
（1））;
end
功能函数：
functionxout=rk4（x,t,tau,derivsRK,param）
half_tau=0.5*tau;
F1=feval（derivsRK,x,t,param）;
t_half=t+half_tau;
xtemp=x+half_tau*F1;
F2=feval（derivsRK,xtemp,t_half,param）;
xtemp=x+half_tau*F2;
F3=feval（derivsRK,xtemp,t_half,param）;
t_full=t+tau;
xtemp=x+tau*F3;
F4=feval（derivsRK,xtemp,t_full,param）;
xout=x+tau/6.*（F1+F4+2.*（F2+F3））;
functionderiv=gravrk（s,t,GM）
r=[s
（1）s
（2）];
v=[s（3）s（4）];
accel=-GM*r/norm（r）^3;
deriv=[v
（1）v
（2）accel
（1）accel
（2）];
clearall;
GM=4*pi^2;
mass=1;
L=2*pi*（1^2）;
r=zeros（2,600）;
r（1,:
）=（0:
599）/120;
r（2,:
）=（0:
599）/120;
R=（r（1,:
）.^2+r（2,:
）.^2）.^0.5;
R2=R.^2;
portential=-GM*mass./R;
Uli=L^2./（2*mass*R2）;
Ueff=Uli+portential;
plot（[-R,R],[portential,portential],'r',[-R,R],[Uli,Uli],'g',[-R,R],[Ueff,Ueff],'b'）;holdon;
plot（0,-20,'.r',0,-10,'.r',0,0,'.r',0,30,'.r'）;
[x,y]=meshgrid（-4:
0.005:
4）;
z=4*pi*（1./（2*（（x.^2+y.^2）.^.5）.^2）-1./（（x.^2+y.^2）.^.5））;
z（z>3）=NaN;
mesh（x,y,z）
holdon;
t=linspace（0,2*pi,1000）;
x1=1.9*sin（t）+1;
y1=1.9*cos（t）+1;
z1=4*pi*（1./（2*（（x1.^2+y1.^2）.^.5）.^2）-1./（（x1.^2+y1.^2）.^.5））;
plot3（x1,y1,z1,'w','LineWidth',3）;
h1=line（'Color',[0,0,0],'Marker','.','MarkerSize',20,'EraseMode','xor'）;
i=1;j=1;
while1;
set（h1,'xdata',x1（i）,'ydata',y1（i）,'zdata',z1（i））;
pause（0.01）;
i=i+1;
ifi>length（z1）
i=1;j=j+1;
end
end
单摆EJP
rm=0.0083;g=9.8;l=1.36;theta0=14*pi/180;fy0=0.545;
ydot=@（t,y）[y
（2）;
-rm*y
（2）+sin（y
（1））*cos（y
（1））*y（4）^2-g/l*sin（y
（1））;
y（4）;
-rm*y（4）-2*cot（y
（1））*y
（2）*y（4）];
[t,y]=ode45（ydot,[0,300],[theta0;0;0;fy0]）;
x=l.*sin（y（:
1））.*cos（y（:
3））;
y=l.*sin（y（:
1））.*sin（y（:
3））;
plot（x,y,'y'）
axis（[-0.40.4-0.40.4]）;
set（gca,'color',[000]）
rm=0.0083;g=9.8;l=1.36;theta0=20*pi/180;fy0=0.028;
ydot=@（t,y）[y
（2）;
-rm*y
（2）+sin（y
（1））*cos（y
（1））*y（4）^2-g/l*sin（y
（1））;
y（4）;
-rm*y（4）-2*cot（y
（1））*y
（2）*y（4）];
[t,y]=ode45（ydot,[0,300],[theta0;0;0;fy0]）;
x=l.*sin（y（:
1））.*cos（y（:
3））;
y=l.*sin（y（:
1））.*sin（y（:
3））;
plot（x,y,'y'）
axis（[-0.40.4-0.40.4]）;
set（gca,'color',[000]）
rm=0.0083;g=9.8;l=1.36;theta0=25*pi/180;fy0=0.733;
ydot=@（t,y）[y
（2）;
-rm*y
（2）+sin（y
（1））*cos（y
（1））*y（4）^2-g/l*sin（y
（1））;
y（4）;
-rm*y（4）-2*cot（y
（1））*y
（2）*y（4）];
[t,y]=ode45（ydot,[0,300],[theta0;0;0;fy0]）;
x=l.*sin（y（:
1））.*cos（y（:
3））;
y=l.*sin（y（:
1））.*sin（y（:
3））;
plot（x,y,'y'）
axis（[-0.40.4-0.40.4]）;
set（gca,'color',[000]）
升级会员 