计算物理讲稿13章

1、计算物理计算物理 第一章引言第一章引言 计算物理（计算物理（computational physics) 是一门新兴的边缘是一门新兴的边缘 学科，是物理学、数学、计算机科学三者相结合的产物。计算学科，是物理学、数学、计算机科学三者相结合的产物。计算 物理也是物理学的一个分支，它与理论物理、实验物理有密切物理也是物理学的一个分支，它与理论物理、实验物理有密切 联系，但又保持着自己的相对的独立性。可以这样给计算物理联系，但又保持着自己的相对的独立性。可以这样给计算物理 下一个定义：计算物理是以计算机为工具，以相关算法为手段，下一个定义：计算物理是以计算机为工具，以相关算法为手段， 解决复杂物理问题

2、的一门应用科学。计算物理已经给复杂体系解决复杂物理问题的一门应用科学。计算物理已经给复杂体系 的物理规律、物理性质的研究提供了重要手段，对物理学的物理规律、物理性质的研究提供了重要手段，对物理学 的发展起着极大的推动作用。的发展起着极大的推动作用。 90年代初期，在我国许多大学的应用物理系纷纷开设年代初期，在我国许多大学的应用物理系纷纷开设 了计算物理课程，受到师生们的欢迎。它对训练学生开阔了计算物理课程，受到师生们的欢迎。它对训练学生开阔 的思维、培养综合分析的能力和获得广博实用的知识大有的思维、培养综合分析的能力和获得广博实用的知识大有 益处，同时对理论教学提供了一个实践的过程，使得以往益

3、处，同时对理论教学提供了一个实践的过程，使得以往 抽象的理论教学形象化。抽象的理论教学形象化。 本课程面向本科生教学，学时为本课程面向本科生教学，学时为42课时，其中需用课时，其中需用20 个左右的课时上机实践。本课程编程是基于个左右的课时上机实践。本课程编程是基于matlab程序的，程序的， 需要学生有良好的需要学生有良好的matlab编程能力。编程能力。 计算物理通常涉及各类线性和非线性方程（组）的求根、计算物理通常涉及各类线性和非线性方程（组）的求根、 数值积分、常微分方程和偏微分方程的求解、另外还包括付数值积分、常微分方程和偏微分方程的求解、另外还包括付 里叶变换、信号处理等内容。本课

4、程教学将涉及微分方程、里叶变换、信号处理等内容。本课程教学将涉及微分方程、 偏微分方程的求解、付里叶变换和信号处理（主要介绍滤偏微分方程的求解、付里叶变换和信号处理（主要介绍滤 波）。波）。 第二章物理学中常微分方程及其计算方法第二章物理学中常微分方程及其计算方法 科学计算中常常需要求解中常微分方程的定解问题，这科学计算中常常需要求解中常微分方程的定解问题，这 类问题最简单的形式是一阶方程的初值问题：类问题最简单的形式是一阶方程的初值问题： 虽然求解常微分方程有各种各样的解析方程，但解析方法只虽然求解常微分方程有各种各样的解析方程，但解析方法只 能用来求解一些特殊类型的方程，求解从实际问题中归

5、结出能用来求解一些特殊类型的方程，求解从实际问题中归结出 来的微分方程主要靠数值解法。来的微分方程主要靠数值解法。 00 )( ),( yxy yxfy 1、欧拉（、欧拉（euler)方法方法 数值方法的第一步就是将微分方程中的导数项数值方法的第一步就是将微分方程中的导数项y进行离散进行离散 化。设在区间化。设在区间xn,xn+1的左端点的左端点xn,则：则： y(xn)=f(xn,y(xn) 并用差商并用差商 替代导数项替代导数项y(xn),则有则有 或写成或写成 h xyxy nn )()( 1 )(,()()( 1nnnn xyxhfxyxy , 2 , 1 , 0),( 1 nyxhf

6、yy nnnn 这就是著名的这就是著名的euler格式，若初值格式，若初值y0已知，在取定步长已知，在取定步长h后，就后，就 可以逐步叠代算出数值解可以逐步叠代算出数值解y1，y2.。 实际应用中实际应用中euler格式格式 存在较大的误差，为此人们又提出了各种改进的存在较大的误差，为此人们又提出了各种改进的euler格式。格式。 其中有一种改进的其中有一种改进的euler格式是：格式是： ),( ),( )( 112 1 2121 hkyxfk yxfk kkyy nn nn h nn 2、龙格库塔（、龙格库塔（runge-kutta)方法方法 2.1 runge-kutta方法的设计思想方

7、法的设计思想 差商考察，根据微分中值定理，存在一点差商考察，根据微分中值定理，存在一点 ，因此，因此 设，称作区间设，称作区间xn， ，xn+1上的平均斜率，这样， 上的平均斜率，这样， 只要对平均斜率提供一种算法，就可以导出相应的一种计算只要对平均斜率提供一种算法，就可以导出相应的一种计算 h xyxy nn )()( 1 1 , nn xx )( )()( 1 y h xyxy nn )(,()()( 1 yhfxyxy nn )(,(*yfk 格式。实际是欧拉格式简单地取点格式。实际是欧拉格式简单地取点xn处的斜率值处的斜率值f(xn,yn)作为平作为平 均斜率均斜率k ，精度自然很低。

8、改进的欧拉格式是用 ，精度自然很低。改进的欧拉格式是用xn与与xn+1两个两个 点的斜率值点的斜率值k1和和k2算术平均作为平均斜率算术平均作为平均斜率k ， ，而而xn+1处的斜率处的斜率 值值k2则利用已知信息则利用已知信息yn通过通过euler格式来预报。格式来预报。 如果我们设法在如果我们设法在xn,xn+1内多报几个点的斜率值，然后将它们内多报几个点的斜率值，然后将它们 加权平均作为平均斜率加权平均作为平均斜率k ， ，则可能构造出更高精度的计算格则可能构造出更高精度的计算格 式，这就是式，这就是runge-kutta方法的设计思想。方法的设计思想。 2.2 龙格库塔（龙格库塔（ru

9、ngekutta）格式）格式 如果如果y(x)在在xn,xn+1上有若干钭率值，即所谓的预报钭率值上有若干钭率值，即所谓的预报钭率值 k1,k2kr以及权数以及权数a1,a2,.ar则：则： 现在最常用的是四阶现在最常用的是四阶rk格式：格式： i i inn kahyy 1 ),( ),( ),( ),( )22( 314 222/13 122/12 1 432161 hkyxfk kyxfk kyxfk yxfk kkkkyy nn h nn h nn nn h nn 这一格式用这一格式用4个点个点 1 , 2 1 2 1 n nn n xxxx 的斜率值的斜率值 4321 ,kkkk加权

10、加权 平均生成平均斜率。通过计算机语言编程就可以求解这样的常平均生成平均斜率。通过计算机语言编程就可以求解这样的常 微分方程。在一个实际工作中，选择何种计算方法要根据实际微分方程。在一个实际工作中，选择何种计算方法要根据实际 情况而定，基本原则是：情况而定，基本原则是：(1)满足需要的精度，满足需要的精度，(2)节省机时。节省机时。 %微分方程： y=y-2*x/y, 0=x=1 %初始条件： y(0)=1 dyfun=inline(y-2*x/y); x,y=rk4(dyfun,0 1,1,0.1); x,y plot(x,y) 龙格库塔法解微分方程程序龙格库塔法解微分方程程序: funct

11、ion x,y=rk4(dyfun,xspan,y0,h) x=xspan(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; end x=x;y=y; 第三章常微分方程（组）的第三章常微分方程（组）的

12、matlab解法解法 及其在物理学中的应用及其在物理学中的应用 matlab有着强大的解常微分方程（组）功能，从方法上来有着强大的解常微分方程（组）功能，从方法上来 讲可分为符号法和数值法，特别是数值法应用范围更加广泛。讲可分为符号法和数值法，特别是数值法应用范围更加广泛。 3.1 matlab微分方程符号解法（回顾复习）微分方程符号解法（回顾复习） r=dsolve(eq1,eq2,cond1,cond2,.,v) eq1,eq2,.表示微分方程，表示微分方程，v为独立变量（默认的独立变量为为独立变量（默认的独立变量为 t )，cond1, cond2, . 表示初始条件（或者边界条件）表示

13、初始条件（或者边界条件）. d表示微分算子（表示微分算子（d/dt), d2, d3等表示二次、三次等微分，由等表示二次、三次等微分，由 dsolve()求出的只能是解析解，如没有解析解，需用数值法解求出的只能是解析解，如没有解析解，需用数值法解 2 1 dy y dx 例例: 解微分方程解微分方程 ，初始条件：，初始条件： 0 |1 x y y=dsolve(dy=1+y2,y(0)=1) y = tan(t+1/4*pi) 例例: 解微分方程解微分方程 2 00 2 cos(2 ),|1,|0 xx d ydy xyy dxdx y=dsolve(d2y=cos(2*x)-y,y(0)=1

14、,dy(0)=0,x) y = 4/3*cos(x)-1/3*cos(2*x) 例：解微分方程组例：解微分方程组 43 34 dy xy dt dx xy dt 初始条件：初始条件：(0)0, (0)1xy 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、过阻尼状态；、过阻尼状态；

15、 3、临界阻尼状态。、临界阻尼状态。 弹簧振子系统中质点受弹性力以及液体对其的阻力作用弹簧振子系统中质点受弹性力以及液体对其的阻力作用 ,根据根据 牛顿第二定律牛顿第二定律 x 2 2 0 2 20 d xdx x dtdt 是振动系统的固有频率，是振动系统的固有频率， 为阻尼因数，与振动系统以及为阻尼因数，与振动系统以及 媒介的性质有关媒介的性质有关 0 0 0 0 参数设置：参数设置： 、欠阻尼状态：、欠阻尼状态： 、过阻尼状态：、过阻尼状态： 、临界阻尼状态：、临界阻尼状态： 5 0 5 0 5 0 7 . 0 5 .14 5 x=dsolve(d2x+2*b*dx+w02*x=0,dx

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

17、2 2 0 0 (0)0.1 |0 t d x x dt x dx dt x=dsolve(d2x+w2*x=0,dx(0)=0,x(0)=0.1,t) v=diff(x,t); a=diff(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、mathemat

18、ica等程序设计软件出现以前，常微分方 程数值解主要是通过rk算法由fortran或c语言编程来完成的， 工作量大、效率低，特别是对一些较为复杂的问题求解更是困 难。在matlab 中提供了求解微分方程的函数odemn (如 ode23,ode45等），通过函数调用，可以很方便的求解各种类 型的线性和非线性常微分方程（或者方程组）。matlab 可以求 解一阶常微分方程的初值问题和边界问题，通过改变方程的 形式，可以求解高阶问题。 基本格式x, y=ode45(odefun, xspace, y0, ,p1,p2,.) odefun是与微分方程有关的函数, xspace是变量取值范围， y0

19、是初如条件，p1,p2,.是传递参数。 （1）简单显式：如 y=f(t, y)，可以通过inline函数ode函数 例如：y=y-2x/y,y(0)=1 odefun=inline(y-2*x/y);用inline构造函数odefun x, y=ode45(odefun, 0, 1, 1) (2)复杂问题（主要指二阶以上微分方程以及方程组 ) 此类方程（组）需要首先建立一个关于一阶导数的函数，函 数程序由function引导，所以对于二阶、三阶等高阶方程，必 须运用数学变量替换将高阶（大于2阶）的方程（组）写成一 阶微分方程组，这是求解的关键过程。解ode的基本过程如下： u根据物理的规律，用

20、微分方程与初始条件进行描述 f(y,y,y,y(n)=0, y(0)=y0, y(0)=y1,y(n-1)(0)=yn-1 y0=y0, y1,yn-1 ; %初始条件向量 u变量替换：y1 =y, y2=y1,y3=y2,.yn=yn-1 形成由一阶导数组成的向量： . . . . . . 2 1 )( n n y y y y y y ydot u用一阶导数矩阵建立函数文件,如odefun.m。 u建立一个m文件，将函数文件与初始条件传递给求解器（如 ode45），运行后就可得到在指定区间上的解列向量y. 首先以前面的阻尼振动为例，讨论二阶常微分方程的数值解的 求解器ode45的使用方法。

21、微分方程为需将二阶转化为一阶 构建初如条件向量构建一阶导数向量： 02 2 2 2 0 y dt dy dt yd y(1)=y y(2)=dy(1)/dt dy(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.2,5; tit1=欠阻尼状态欠阻尼状态; tit2=过阻尼状态过阻尼状态; tit3=临界阻尼状态临界阻尼状态;

22、 y0=1 5;%这里也可以写成列向量这里也可以写成列向量y01;5 for i=1:3 t,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沿水 平 方 向 ， o y 竖 直 向 上 。 初 始 条 件 ： t = 0 时 ， x=0,dx/dt=5,

23、y=0,dy/dt=8 y x 质点受重力和空气阻力作用，而空气阻力包括三种情况。质点 的运动微分方程可表示为： 2 222 () 2 y p bv dy gvv xy m dt 2 222 () 2 p bv dx x vv xy m dt 投影式： vvbgm dt rd m p | 2 2 参数值：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解常微分方程 令 ，将二阶微分方程组方程组写成一阶 微分方程组 (1)(2)(3)(4)y =x,y

24、=dx/dt,y =y,y =dy/dt (1) (2) (2)(2) 22 2 (2)(4) (3) (4) (4)(4) 22 2 (2)(4) () () p p dy y dt dyby yy dtm dy y dt dyby gyy dtm 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).2).(p/2); 、主程序： m=1;b=0;p=0; t,y=ode45(kqzlptfun,

25、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,

26、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 2 22 d xklx mkx dt xy 2 22 d ykly mmgky dt xy 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

27、22 22 ) 3() 1 ( ) 3() 3()4( )4( ) 3( ) 3() 1 ( ) 1 () 1 ()2( )2( ) 1 ( yym kly m ky g dt dy y dt dy yym kly m ky dt dy y dt dy 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

28、.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函数方法去 解，我们只要花少量的时间编

29、一下程序，不管多难的问题都可以 解出来，并且可以用图像直观地表示出关系来 。 课堂练习：受迫阻尼振动： )cos(2 2 0 2 2 thx dt dx dt xd 初始条件：t=0时 0/ 2 .0 dtdx x 计算区间0：0.01：100 参数： 005 .0,3 .0 , 4 .0,02.0 pi h 程序：znzdfun1main.m,znzdfun1.m 共振曲线：振幅频率 )5 . 2:1 . 0:0( 曲线：振幅 程序：resonance.m,znzdfun1.m %znzdfun1main.m % x+2bx+w02x=hcos(wt) %w=k*w0 w0=0.3*pi;

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

31、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: 32 2 r r m k dt rd 散射，重核在原点处 2/3222 2 )(yx x m k dt xd 2/3222 2 )(yx y m k dt yd dt dy y yy dt dx y xy )4( )3( )2( ) 1 ( 2/322 2/322 )3() 1 ( )3()4( )4( )3

32、( )3() 1 ( ) 1 ()2( )2( ) 1 ( yy y m k dt dy y dt dy yy y m k dt dy y dt dy 初始条件: 0 , 4 , 1 ,100: ,0 0000 y vyvxy yx 取 3/mk %alzssfunmain y0=-10,1,10,0 plot(0,0,r.,markersize,50); text(2,0,靶粒子,fontsize,14 ); xlabel(x);ylabel(y); hold on axis(-10 20 -20 20) for i=1:1 t,y=ode23(alzssfun,0:.01:32,y0(i,

33、:),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 带电粒子在电磁场中的运动 bvqeq dt rd m 2 2 以场中一点为原点,以e为oy方向,b为oz方向建立oxyz系 0 / / 2 2 2 2 2 2 dt zd dt dx mqbmqe dt yd dt dy mqb dt xd 1 2 2 4 3 4 4 2 5 6 6 / / 0 d y y d t d y q bm y d t d y y d t d y q emq bm y d t d y y d t d y d t dtdzy zy dtdyy yy dtdxy xy / / / 6 5 4 3 2 1 初始条件 20:1.0:0 01