计算物理计算物理第一章 引言第一章 引言 计算物理（ 计算物理（Computational Physics) 是一门新兴的边缘是一门新兴的边缘学科，是物理学、数学、计算机科学三者相结合的产物计算学科，是物理学、数学、计算机科学三者相结合的产物计算物理也是物理学的一个分支，它与理论物理、实验物理有密切物理也是物理学的一个分支，它与理论物理、实验物理有密切联系，但又保持着自己的相对的独立性可以这样给计算物理联系，但又保持着自己的相对的独立性可以这样给计算物理下一个定义：计算物理是以计算机为工具，以相关算法为手段，下一个定义：计算物理是以计算机为工具，以相关算法为手段，解决复杂物理问题的一门应用科学计算物理已经给复杂体系解决复杂物理问题的一门应用科学计算物理已经给复杂体系 的物理规律、物理性质的研究提供了重要手段，对物理学的物理规律、物理性质的研究提供了重要手段，对物理学的发展起着极大的推动作用的发展起着极大的推动作用 90年代初期，在我国许多大学的应用物理系纷纷开设年代初期，在我国许多大学的应用物理系纷纷开设了计算物理课程，受到师生们的欢迎它对训练学生开阔了计算物理课程，受到师生们的欢迎。
它对训练学生开阔的思维、培养综合分析的能力和获得广博实用的知识大有的思维、培养综合分析的能力和获得广博实用的知识大有益处，同时对理论教学提供了一个实践的过程，使得以往益处，同时对理论教学提供了一个实践的过程，使得以往抽象的理论教学形象化抽象的理论教学形象化 本课程面向本科生教学，学时为本课程面向本科生教学，学时为42课时，其中需用课时，其中需用20 个左右的课时上机实践本课程编程是基于个左右的课时上机实践本课程编程是基于Matlab程序的，程序的，需要学生有良好的需要学生有良好的Matlab编程能力编程能力 计算物理通常涉及各类线性和非线性方程（组）的求根、 计算物理通常涉及各类线性和非线性方程（组）的求根、数值积分、常微分方程和偏微分方程的求解、另外还包括付数值积分、常微分方程和偏微分方程的求解、另外还包括付里叶变换、信号处理等内容本课程教学将涉及微分方程、里叶变换、信号处理等内容本课程教学将涉及微分方程、偏微分方程的求解、付里叶变换和信号处理（主要介绍滤波）偏微分方程的求解、付里叶变换和信号处理（主要介绍滤波） 第二章 物理学中常微分方程及其计算方法第二章 物理学中常微分方程及其计算方法 科学计算中常常需要求解中常微分方程的定解问题，这科学计算中常常需要求解中常微分方程的定解问题，这类问题最简单的形式是一阶方程的初值问题：类问题最简单的形式是一阶方程的初值问题：虽然求解常微分方程有各种各样的解析方程，但解析方法只虽然求解常微分方程有各种各样的解析方程，但解析方法只能用来求解一些特殊类型的方程，求解从实际问题中归结出能用来求解一些特殊类型的方程，求解从实际问题中归结出来的微分方程主要靠数值解法。
来的微分方程主要靠数值解法 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方法的设计思想方法的设计思想差商考察 ，根据微分中值定理，存在一点 差商考察 ，根据微分中值定理，存在一点 ，因此 ，因此 设设 ，，称称作作区区间间[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 龙格－库塔（龙格－库塔（Runge－－Kutta）格式）格式如果如果y(x)在在[xn,xn+1]上有若干钭率值，即所谓的预报钭率值上有若干钭率值，即所谓的预报钭率值k1,k2……kr以及权数以及权数a1,a2,…..ar则：则：现在最常用的是四阶现在最常用的是四阶R－－K格式：格式： 这一格式用这一格式用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,y]plot(x,y)龙格－库塔法解微分方程程序龙格－库塔法解微分方程程序: function [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;endx=x';y=y'; 第三章 常微分方程（组）的第三章 常微分方程（组）的Matlab解法解法 及其在物理学中的应用 及其在物理学中的应用 Matlab有着强大的解常微分方程（组）功能，从方法上来有着强大的解常微分方程（组）功能，从方法上来讲可分为符号法和数值法，特别是数值法应用范围更加广泛。
讲可分为符号法和数值法，特别是数值法应用范围更加广泛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+y^2’,’y(0)=1’)y = tan(t+1/4*pi)例例: 解微分方程解微分方程y=dsolve('D2y=cos(2*x)-y','y(0)=1,Dy(0)=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('D2x+2*b*Dx+w0^2*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弹弹簧组成，已知簧组成，已知t＝＝0时，时，m在在A处，处，x0=A=0.1m,由静止开始释由静止开始释放，研究：放，研究：x~t, v~t, a~t关系，作图表示。
关系，作图表示涉及的微分方程和初始条件：涉及的微分方程和初始条件： x=dsolve('D2x+w^2*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、Mathematica等程序设计软件出现以前，常微分方程数值解主要是通过R－K算法由Fortran或C语言编程来完成的，工作量大、效率低，特别是对一些较为复杂的问题求解更是困难在Matlab 中提供了求解微分方程的函数odemn (如ode23,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构造函数odefun[x, y]=ode45(odefun, [0, 1], 1) (2)复杂问题（主要指二阶以上微分方程以及方程组 )此类方程（组）需要首先建立一个关于一阶导数的函数，函数程序由function引导，所以对于二阶、三阶等高阶方程，必须运用数学变量替换将高阶（大于2阶）的方程（组）写成一阶微分方程组，这是求解的关键过程解ode的基本过程如下：u根据物理的规律，用微分方程与初始条件进行描述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’形成由一阶导数组成的向量： u用一阶导数矩阵建立函数文件,如odefun.m。
u建立一个M文件，将函数文件与初始条件传递给求解器（如ODE45），运行后就可得到在指定区间上的解列向量y.首先以前面的阻尼振动为例，讨论二阶常微分方程的数值解的求解器ode45的使用方法微分方程为 需将二阶转化为一阶①构建初如条件向量 ②构建一阶导数向量： y(1)=yy(2)=dy(1)/dtdy(2)/dt=-2*b*y(2)-w^2*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];tit{1}='欠阻尼状态欠阻尼状态';tit{2}='过阻尼状态过阻尼状态';tit{3}='临界阻尼状态临界阻尼状态';y0=[1 5];%这里也可以写成列向量这里也可以写成列向量y0＝＝[1;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(tit{i},'color','b')end 例例1：：研究有空气阻力时抛运动的特征。
比较下面三种情况下的抛体的轨道：①、没有空气阻力；②、空气阻力与速度一次方成正比；③、空气阻力与速度二次方成正比1、以地面为参考系，以抛点为原点O建立直角坐标系OXY,OX沿水平 方 向 ， OY竖 直 向 上 初 始 条 件 ： t=0时 ，x=0,dx/dt=5,y=0,dy/dt=8 Yx 质点受重力和空气阻力作用，而空气阻力包括三种情况质点的运动微分方程可表示为： 投影式：参数值：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).^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;[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点，写出方程: xyomgvOx: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*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)-w0^2*y(1)+h*cos(k*w0*t)]; Resonance.m% x''+2bx'+w0^2x=hcos(wt)%w=k*w0w0=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))];endplot(k,a);xlabel('\omega /\omega_0')ylabel('振幅') 小课题1:散射，重核在原点处 初始条件: %alzssfunmainy0=[-10,1,10,0]plot(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(@alzssfun,[0:.01:32],y0(i,:),[],3);comet(y(:,1),y(:,3))endfunction 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],[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 on xlabel('x'); ylabel('y'); zlabel('z');endfunction 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 Matlab边界值问题的常微分方程的数值解法两点边值微分方程问题，也是物理学和工程中常会遇到的如：如解边值问题：步骤：1、建立一阶微分方程（组）odefun函数2、建立边值函数bcfun 3、给出节点预估值网络(调用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]; xy 作业：输出时的y值，并作y(t)图 。