计算物理课件1–3章.ppt

⑵、主程序： m=1;b=0;p=0; [t,y]=ode45(＠kqzlptfun,[0:0.001:2],[0,5,0,8],[],b,p,m); subplot(3,1,1) plot(y(:,1),y(:,3)); title(无空气阻力抛体运动,color,b) %%加标题 hold on %%保持图形窗口继续画图 subplot(3,1,2) m=1;b=0.1;p=0; [t,y]=ode45(＠kqzlptfun,[0:0.001:2],[0,5,0,8],[],b,p,m); plot(y(:,1),y(:,3)); title(空气阻力与速度一次方成正比,color,b) subplot(3,1,3) m=1;b=0.1;p=1; [t,y]=ode45(＠kqzlptfun,[0:0.001:2],[0,5,0,8],[],b,p,m); plot(y(:,1),y(:,3)); title(空气阻力与速度二次方成正比,color,b) 例2：弹簧小球的非线性运动：质量为m=0.1kg的小球，挂在原长为l=1.0m，劲度系数为k=4.8N/m的轻弹簧的一端，弹簧的另一端固定。 1、建立如图坐标。设小球某时刻t位于p点，写出方程: x y o mg v Ox: Oy: 2、令y(1) =x,　 y(2) = dx/dt,　 y(3) = y,　 y(4) = dy/dt; 　　　m=0.1,k=4.8,l=1.0,g=9.8 3、程序 ⑴、函数程序： function ydot=tanhuangfun(t,y, m,k,l,g); ydot=[y(2); -k*y(1)/m+k*l*y(1)/(m*sqrt(y(1).^2+y(3).^2)); y(4); g-k*y(3)/m+k*l*y(3)/(m*sqrt(y(1).^2+y(3).^2))]; （2）主程序 m=0.1;k=4.8;l=1.0;g=9.8; [t,y]=ode45(＠tanhuangfun,[0:0.001:30],[1,0,0,0],[],m,k,l,g); plot(y(:,1),y(:,3)) title(弹簧小球的非线性运动轨迹,color,b) xlabel(X);ylabel(Y) 　　在具体的实际问题中，常遇到给定初值的常微分方程（组），MATLAB对解决这类问题有着独到之处。对于简单的问题（通常有解析式），我们可以用符号法解，即调用dsolve函数。而对于复杂的问题，有时我们就要花许多时间才能解出最终关系式，或者我们根本就解不出，此时需要用MATLAB的ODE45函数方法去解，我们只要花少量的时间编一下程序，不管多难的问题都可以解出来，并且可以用图像直观地表示出关系来 。 课堂练习：受迫阻尼振动： 初始条件：t=0时 计算区间　0：0.01：100 参数： 程序：znzdfun1main.m,znzdfun1.m 共振曲线：振幅－频率 曲线：振幅－ 程序：resonance.m,znzdfun1.m %znzdfun1main.m % x+2bx+w0^2x=hcos(wt) %w=k*w0 w0=0.3*pi; b=0.02; k=0.5; h=0.4; t=[0:0.01:100]; y0=[0.2 0]; [t,y]=ode45(@znzdfun,t,y0,[],w0,b,k,h); comet(t,y(:,1)); xlabel(t) ylabel(位移) %znzdfun1.m function ydot=znzdfun1(t,y,w0,b,k,h) ydot=[y(2);-2*b*y(2)-w0^2*y(1)+h*cos(k*w0*t)]; Resonance.m % x+2bx+w0^2x=hcos(wt) %w=k*w0 w0=0.3*pi; b=0.02; h=0.4; a=[]; k=[0:0.1:2.5]; t=[0:0.1:100]; y0=[0.2 0]; for i=1:length(k) [t,y]=ode45(@znzdfun,t,y0,[],w0,b,k(i),h); a=[a,max(y(:,1))]; end plot(k,a); xlabel(\omega /\omega_0) ylabel(振幅) 小课题1: 散射， 重核在原点处 初始条件: %alzssfunmain y0=[-10,1,10,0] plot(0,0,r.,MarkerSize,50); text(2,0,靶粒子,fontsize,14 ); xlabel(x);ylabel(y);