计算物理讲稿13章.ppt

《计算物理讲稿13章.ppt》由会员分享，可在线阅读，更多相关《计算物理讲稿13章.ppt（53页珍藏版）》请在三一办公上搜索。

1、计算物理第一章引言计算物理（Computational Physics)是一门新兴的边缘学科，是物理学、数学、计算机科学三者相结合的产物。计算物理也是物理学的一个分支，它与理论物理、实验物理有密切联系，但又保持着自己的相对的独立性。可以这样给计算物理下一个定义：计算物理是以计算机为工具，以相关算法为手段，解决复杂物理问题的一门应用科学。计算物理已经给复杂体系,的物理规律、物理性质的研究提供了重要手段，对物理学的发展起着极大的推动作用。90年代初期，在我国许多大学的应用物理系纷纷开设了计算物理课程，受到师生们的欢迎。它对训练学生开阔的思维、培养综合分析的能力和获得广博实用的知识大有益处，同时对理

2、论教学提供了一个实践的过程，使得以往抽象的理论教学形象化。本课程面向本科生教学，学时为42课时，其中需用20,个左右的课时上机实践。本课程编程是基于Matlab程序的，需要学生有良好的Matlab编程能力。计算物理通常涉及各类线性和非线性方程（组）的求根、数值积分、常微分方程和偏微分方程的求解、另外还包括付里叶变换、信号处理等内容。本课程教学将涉及微分方程、偏微分方程的求解、付里叶变换和信号处理（主要介绍滤波）。,第二章物理学中常微分方程及其计算方法 科学计算中常常需要求解中常微分方程的定解问题，这类问题最简单的形式是一阶方程的初值问题：虽然求解常微分方程有各种各样的解析方程，但解析方法只能用

3、来求解一些特殊类型的方程，求解从实际问题中归结出来的微分方程主要靠数值解法。,1、欧拉（Euler)方法 数值方法的第一步就是将微分方程中的导数项y进行离散化。设在区间xn,xn+1的左端点xn,则：y(xn)=f(xn,y(xn)并用差商 替代导数项y(xn),则有或写成,这就是著名的Euler格式，若初值y0已知，在取定步长h后，就可以逐步叠代算出数值解y1，y2.。实际应用中Euler格式存在较大的误差，为此人们又提出了各种改进的Euler格式。其中有一种改进的Euler格式是：,2、龙格库塔（Runge-Kutta)方法2.1 Runge-Kutta方法的设计思想差商考察，根据微分中值

4、定理，存在一点，因此设，称作区间xn，xn+1上的平均斜率，这样，只要对平均斜率提供一种算法，就可以导出相应的一种计算,格式。实际是欧拉格式简单地取点xn处的斜率值f(xn,yn)作为平均斜率K，精度自然很低。改进的欧拉格式是用xn与xn+1两个点的斜率值K1和K2算术平均作为平均斜率K，而xn+1处的斜率值K2则利用已知信息yn通过Euler格式来预报。如果我们设法在xn,xn+1内多报几个点的斜率值，然后将它们加权平均作为平均斜率K，则可能构造出更高精度的计算格式，这就是Runge-Kutta方法的设计思想。,2.2 龙格库塔（RungeKutta）格式如果y(x)在xn,xn+1上有若干

5、钭率值，即所谓的预报钭率值k1,k2kr以及权数a1,a2,.ar则：现在最常用的是四阶RK格式：,这一格式用4个点,的斜率值,加权,平均生成平均斜率。通过计算机语言编程就可以求解这样的常微分方程。在一个实际工作中，选择何种计算方法要根据实际情况而定，基本原则是：(1)满足需要的精度，(2)节省机时。,%微分方程：y=y-2*x/y,0=x=1%初始条件：y(0)=1dyfun=inline(y-2*x/y);x,y=rk4(dyfun,0 1,1,0.1);x,yplot(x,y),龙格库塔法解微分方程程序:,function x,y=rk4(dyfun,xspan,y0,h)x=xspan

6、(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n);%计算dyfun值 k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x;y=y;,第三章常微分方程（组）的Matlab解法及其在物理学中的应用Matlab有着强大的解常微分方程（组）功能，从方法上来讲可分为符号法和数值法

7、，特别是数值法应用范围更加广泛。3.1 Matlab微分方程符号解法（回顾复习）r=dsolve(eq1,eq2,cond1,cond2,.,v)eq1,eq2,.表示微分方程，v为独立变量（默认的独立变量为 t)，cond1,cond2,.表示初始条件（或者边界条件）.D表示微分算子（d/dt),D2,D3等表示二次、三次等微分，由dsolve()求出的只能是解析解，如没有解析解，需用数值法解,例:解微分方程,，初始条件：,y=dsolve(Dy=1+y2,y(0)=1),y=tan(t+1/4*pi),例:解微分方程,y=dsolve(D2y=cos(2*x)-y,y(0)=1,Dy(0)

8、=0,x)y=4/3*cos(x)-1/3*cos(2*x),例：解微分方程组,初始条件：,x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y,x(0)=0,y(0)=1)x=exp(3*t)*sin(4*t)y=exp(3*t)*cos(4*t),下面讨论受阻力作用时振动系统的运动特征。比较下面三种情况下振子的轨迹：1、欠阻尼状态；2、过阻尼状态；3、临界阻尼状态。弹簧振子系统中质点受弹性力以及液体对其的阻力作用,根据牛顿第二定律,X,是振动系统的固有频率，为阻尼因数，与振动系统以及媒介的性质有关,参数设置：、欠阻尼状态：,、过阻尼状态：,、临界阻尼状态：,x=dsolve(

9、D2x+2*b*Dx+w02*x=0,Dx(0)=5,x(0)=1);w0=5;t=0:0.1:10;b=0.7;x1=eval(x);plot(t,x1)hold onb=14.5;x1=eval(x);plot(t,x1)hold onb=5;x1=eval(x);plot(t,x1),作业:弹簧振子系统由质量m=2kg的滑块，k=400N/m弹簧组成，已知t0时，m在A处，x0=A=0.1m,由静止开始释放，研究：xt,vt,at关系，作图表示。,涉及的微分方程和初始条件：,x=dsolve(D2x+w2*x=0,Dx(0)=0,x(0)=0.1,t)v=diff(x,t);a=diff

10、(x,t,2);k=400;m=2;w=sqrt(k/m);t=0:0.01:0.9;x1=eval(x);v1=eval(v);a1=eval(a);subplot(3,1,1)plot(t,x1)subplot(3,1,2)plot(t,v1)subplot(3,1,3)plot(t,a1),3.2 Matlab初值问题的常微分方程的数值解法在Matlab、Mathematica等程序设计软件出现以前，常微分方程数值解主要是通过RK算法由Fortran或C语言编程来完成的，工作量大、效率低，特别是对一些较为复杂的问题求解更是困难。在Matlab 中提供了求解微分方程的函数odemn(如od

11、e23,ode45等），通过函数调用，可以很方便的求解各种类型的线性和非线性常微分方程（或者方程组）。Matlab 可以求解一阶常微分方程的初值问题和边界问题，通过改变方程的,形式，可以求解高阶问题。基本格式x,y=ode45(odefun,xspace,y0,p1,p2,.)odefun是与微分方程有关的函数,xspace是变量取值范围，y0是初如条件，p1,p2,.是传递参数。（1）简单显式：如 y=f(t,y)，可以通过inline函数ode函数例如：y=y-2x/y,y(0)=1odefun=inline(y-2*x/y);用inline构造函数odefunx,y=ode45(odef

12、un,0,1,1),(2)复杂问题（主要指二阶以上微分方程以及方程组)此类方程（组）需要首先建立一个关于一阶导数的函数，函数程序由function引导，所以对于二阶、三阶等高阶方程，必须运用数学变量替换将高阶（大于2阶）的方程（组）写成一阶微分方程组，这是求解的关键过程。解ode的基本过程如下：根据物理的规律，用微分方程与初始条件进行描述F(y,y,y,y(n)=0,y(0)=y0,y(0)=y1,y(n-1)(0)=yn-1,y0=y0,y1,yn-1;%初始条件向量变量替换：y1=y,y2=y1,y3=y2,.yn=yn-1形成由一阶导数组成的向量：,用一阶导数矩阵建立函数文件,如odef

13、un.m。建立一个M文件，将函数文件与初始条件传递给求解器（如ODE45），运行后就可得到在指定区间上的解列向量y.首先以前面的阻尼振动为例，讨论二阶常微分方程的数值解的求解器ode45的使用方法。微分方程为需将二阶转化为一阶构建初如条件向量构建一阶导数向量：,y(1)=yy(2)=dy(1)/dtdy(2)/dt=-2*b*y(2)-w2*y(1)初值向量：y0=1,5一阶导数向量：ydot=y(2);-2*b*y(2)-w*y(1),函数文件function ydot=znzdfun(t,y,w,b)ydot=y(2);-2*b*y(2)-w.2*y(1);主程序w=5;b=0.7,24.

14、2,5;tit1=欠阻尼状态;tit2=过阻尼状态;tit3=临界阻尼状态;y0=1 5;%这里也可以写成列向量y01;5for i=1:3t,y=ode45(znzdfun,0:0.1:10,y0,w,b(i);subplot(3,1,i)plot(t,y(:,1);title(titi,color,b)end,例1：研究有空气阻力时抛运动的特征。比较下面三种情况下的抛体的轨道：、没有空气阻力；、空气阻力与速度一次方成正比；、空气阻力与速度二次方成正比。1、以地面为参考系，以抛点为原点O建立直角坐标系OXY,OX沿水平方向，OY竖直向上。初始条件：t=0时，x=0,dx/dt=5,y=0,d

15、y/dt=8,Y,x,质点受重力和空气阻力作用，而空气阻力包括三种情况。质点的运动微分方程可表示为：,投影式：,参数值：b=0,0.1,0.1,p=0,0,1，(b=0,p=0表示无空气阻力；b=0.1,p=0表示空气阻力与速度一次方成正比；b=0.1,p=1表示空气阻力与速度二次方成正比)。,2、用ODE解常微分方程令，将二阶微分方程组方程组写成一阶微分方程组,3、程序、函数文件：function ydot=kqzlptfun(t,y,b,p,m)ydot=y(2);-b/m*y(2)*(y(2).2+y(4).2).(p/2);y(4);-9.8-b/m*y(4)*(y(2).2+y(4)

16、.2).(p/2);,、主程序：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;

17、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点，写出方程:,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

18、,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对解决这类问题有着独到

19、之处。对于简单的问题（通常有解析式），我们可以用符号法解，即调用dsolve函数。而对于复杂的问题，有时我们就要花许多时间才能解出最终关系式，或者我们根本就解不出，此时需要用MATLAB的ODE45函数方法去解，我们只要花少量的时间编一下程序，不管多难的问题都可以解出来，并且可以用图像直观地表示出关系来。,课堂练习：受迫阻尼振动：,初始条件：t=0时,计算区间0：0.01：100,参数：,程序：znzdfun1main.m,znzdfun1.m,共振曲线：振幅频率,曲线：振幅,程序：resonance.m,znzdfun1.m,%znzdfun1main.m%x+2bx+w02x=hcos(w

20、t)%w=k*w0w0=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.mfunction ydot=znzdfun1(t,y,w0,b,k,h)ydot=y(2);-2*b*y(2)-w02*y(1)+h*cos(k*w0*t);,Resonance.m%x+2bx+w02x=hcos(wt)%w=k*w0w0=0.3*pi;b=0.02;h=0.4;a=;k=0:0.1:2.5;t=0:0

21、.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);endplot(k,a);xlabel(omega/omega_0)ylabel(振幅),小课题1:,散射，,重核在原点处,初始条件:,%alzssfunmainy0=-10,1,10,0plot(0,0,r.,MarkerSize,50);text(2,0,靶粒子,fontsize,14);xlabel(x);ylabel(y);hold onaxis(-10 20-20 20)for i=1:1 t,y=ode23(alzs

22、sfun,0:.01:32,y0(i,:),3);comet(y(:,1),y(:,3)end,function ydot=alzssfun(t,y,p)ydot=y(2);p*y(1)/sqrt(y(1).*y(1)+y(3).*y(3).3;y(4);p*y(3)/sqrt(y(1).*y(1)+y(3).*y(3).3;,小课题2,带电粒子在电磁场中的运动,以场中一点为原点,以E为oy方向,B为oz方向建立oxyz系,初始条件,%dccfunmain.mq=1.6e-2;m=0.02;B=2;1;0;E=1;0;1;for i=1:3 t,y=ode45(dccfun,0:0.1:20,

23、0,0.01,0,6,0,0.01,q,m,B(i),E(i);subplot(1,3,i)plot3(y(:,1),y(:,3),y(:,5),linewidth,2);grid onxlabel(x);ylabel(y);zlabel(z);end,function ydot=dccfun(t,y,q,m,b,e)ydot=y(2);q*b*y(4)/m;y(4);q*e/m-q*b*y(2)/m;y(6);0;,作业：求解描述振荡器的经典的Ver der Pol微分方程d2y/dt2-u(1-y)dy/dt+y=0,y(0)=1,y(0)=0.试分析t与y,t与y关系(取u=7),3.3

24、 Matlab边界值问题的常微分方程的数值解法两点边值微分方程问题，也是物理学和工程中常会遇到的如：,如解边值问题：,步骤：,1、建立一阶微分方程（组）odefun函数2、建立边值函数bcfun3、给出节点预估值网络(调用bvpinit函数）4、调用bvp4c函数完成求解,边界条件是由边值函数描述的，基本格式：res=bcfun(ya,yb),ya和yb分别表示y(a)和y(b)计算中要先给出预估初值，基本格式：solinit=bvpinit(tinit,yinit),y(1)=z,y(2)=dz/dt,y(1)=y(2)y(2)=-|y(1)|,%twobcode.msolinit=bvpinit(linspace(0,4,20),1;0);sol=bvp4c(twoode,twobc,solinit);x=linspace(0,4,10);y=deval(sol,x);plot(x,y(1,:),sol.x,sol.y(1,:),o);_function dydx=twoode(x,y)dydx=y(2)-abs(y(1);_ function res=twobc(ya,yb)res=ya(1);yb(1)+2;,x,y,作业：,输出,时的y值，并作y(t)图,