北邮计算物理参考代码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])