北邮计算物理参考代码.docx

1、北邮计算物理参考代码测试第一题：二阶RK2函数：function out2 = 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函数：function out4 = RK4(x,time,tau,driverRK3) F1 = feval(driverRK3,x,time); t_half = time + tau/2;

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函数：function drive = gravrk(y,t)global F m gama1 gama2 k K w;

3、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) accel1 vz(2) accel2;end二阶RK算法主函数：clear allglobal F m gama1 gama2 k K w;F = 1;m

4、 = 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 = 0 0 0 0;for i = 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

5、) = state(4); timeS(i) = time;endplot(timeS,Sx,-b,timeS,Sy,-r,timeS,Sv1,:m,timeS,Sv2,-.k)四阶RK算法主函数：clear allglobal F m gama1 gama2 k K w;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 = 0 0 0 0;for i = 1:nStep state = RK4(state,

6、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;endplot(timeS,Sx,-b,timeS,Sy,-r,timeS,Sv1,:m,timeS,Sv2,-.k)利用ode45主函数：clear allglobal F m gama1 gama2 k K w;F = 1;m = 1.55;gam

7、a1 = 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);幅度计算功能函数：function fupin = fupinq

8、(w)global F m gama1 gama2 k K;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作幅频曲线主函数：clear allfor i = 0:40

9、0 theta = pi*i/10; fu = fupinq(theta); plot(theta,fu,MarkerSize,3); hold onend动图主函数：clear allfor theta1 = 1:800 fu(theta1) = fupinq(theta1/1000); theta(theta1) = theta1;end comet(theta,fu);、测试第二题function out4 = RK4(x,time,tau,driverRK3,beta) F1 = feval(driverRK3,x,time,beta); t_half = time + tau/2;

10、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);endfunction drive = gravrk(y,t,beta) global a b aerfa; a

11、= -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 = accel1 accel2 accel3;endfunction main_3_RK4(beta) global a b aerfa; a = -8/7;b = -5/7;aerfa = 9.0; tau = 0.01; nStep = 20*(1/tau); ti

12、me = 0; state = 0 0.01 0.01; for i = 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);endfunction main_3_1_z a = 24 18 16 15.6 15.2 14; for i =1:6 subplot(2,3,i) main_3_RK4(a(i) endendfunction y

13、dot = main_3_ode45(beta) global a b aerfa; 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); hold on xlabel(beta);

14、ylabel(y);endfunction main_3_1 for a = 14 :0.1:25; main_3_ode45(a) endend测试第三题：clear alle = 0.8;x = zeros(5000,51);x(1,:) = rand(1,51);for n =1:5000x(n,1) = 0.2;endfor m = 1:5000for i = 2:51x(m+1,i) = (1-e)*4*x(m,i)*(1-x(m,i)+e*4*x(m,i-1)*(1-x(m,i-1);endendsubplot(2,4,1)plot(1:50,x(1,2:51),.);subplo

15、t(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);作业题第一题function djdb2 theta = 1/180*pi 5/180*pi 10/180*pi; T = ;

16、 ops=odeset(Events,events); for k = 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-4endfunction ydot=f(t,y) ydot=y(2);-sin(y(1);endfunction value,isterminal,direction=events(t,y) value=(y(2); is

17、terminal=1; direction=1;end作业题第二题function xiangtu 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() subp

18、lot(2,1,2) plot(w1(:,1),w1(:,2),b); hold on plot(w2(:,1),w2(:,2),g:) hold off title() endfunction ydot = fun_2(t,y,beta) ydot = y(2); -2*beta*y(2)-sin(y(1);end作业题第三题function dbyd_3 u=2/3;a=0.5;zq=3*pi;f=1.0726 1.08829 1.095; for i=1:3 ydot=(t,y)y(2); -sin(y(1)-a*y(2)+f(i)*cos(u*t); T,Y=ode45(ydot,0:

19、zq/100:20*zq,-0.8;2); subplot(3,1,i) plot(T,Y(:,1),r) endend作业题第五题function main_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) hold on plot(y2(:,1),y2(:,2),k) title()end作业题第六题function main_6 global u w v u=0.3;w=0.44;T=2*pi/w;

20、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/n2; dt=t(2)-t(1); freq=1/dt*(1:n/2)./n; plot(freq,power) axis(0 0.4 0 0.05) hold onendfunction main_6_2 global u w v 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

21、),y(:,2)endfunction ydot=fun_6(t,y) global u w v x0=1;w0=1; ydot=y(2); u*(x02-y(1)2)*y(2)-y(1)*w02-v*cos(w*t); end作业题第七题function main_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

22、(:,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);endfunction ydot=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)

23、;end作业题第八题function main_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)endfunction ydot = 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作业题第九题function main_9_1 N=1000;X=zeros(N,1);Y=X; n=0.3;u=1.4; X(1)=rand(1);

24、Y(1)=rand(1); for i=1:N X(i+1)=Y(i)+1-u*X(i)2; Y(i+1)=n*X(i); end plot(X,Y,r)endfunction main_9_2 N=500;X=zeros(N,1);Y=X; n=0.3;u=linspace(1.05,1.082,50); xx=;U=; for ii=1:50 uu=u(ii); U=U,uu; X(1)=rand(1);Y(1)=rand(1); for i=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(:,10

25、0:500),.r) title(x-u?)end作业题第十题function main_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,0 0 0); 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)

26、,y)endfunction main_10_2 a=0.3 0.35 0.375 0.386 0.3909 0.398 0.4 0.411;b=2;c=4; X=;Y=;Z=; for ii=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,0 0 0); X=X;y(:,1);Y=Y;y(:,2);Z=Z;y(:,3); end for j=1:4 figure(1) subplot(2,2,j) plot3(X(j,5000:end),Y(j,5000:end

27、),Z(j,5000:end) end for j=5:8 figure(2) subplot(2,2,j-4) plot3(X(j,5000:end),Y(j,5000:end),Z(j,5000:end) endend作业题第十一题function out4 = 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 = fe

28、val(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);endfunction drive = gravrk(y,t) global a b c; 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 = accel1 accel2 accel3;endfunction main_3_RK4 global a b c; a = 10;b = 8/3;c =1 ; tau = 0.01; nStep = 20*(1/tau); time = 0; state = 0 0.08 1; for i = 1:nStep state = RK4(sta