计算物理讲稿1-3章.ppt

计算物理 第一章 引言 计算物理（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 龙格库塔（rungekutta）格式 如果y(x)在xn,xn+1上有若干钭率值，即所谓的预报钭率值k1,k2kr以及权数a1,a2,ar则： 现在最常用的是四阶rk格式：,这一格式用4个点,的斜率值,加权,平均生成平均斜率。通过计算机语言编程就可以求解这样的常微分方程。在一个实际工作中，选择何种计算方法要根据实际情况而定，基本原则是：(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),龙格库塔法解微分方程程序:,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; end x=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+y2,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+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 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关系，作图表示。,涉及的微分方程和初始条件：,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、mathematica等程序设计软件出现以前，常微分方程数值解主要是通过rk算法由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)=1 odefun=inline(y-2*x/y); 用inline构造函数odefun x, y=ode45(odefun, 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 形成由一阶导数组成的向量：,用一阶导数矩阵建立函数文件,如odefun.m。 建立一个m文件，将函数文件与初始条件传递给求解器（如ode45），运行后就可得到在指定区间上的解列向量y. 首先以前面的阻尼振动为例，讨论二阶常微分方程的数值解的求解器ode45的使用方法。 微分方程为 需将二阶转化为一阶 构建初如条件向量 构建一阶导数向量：,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=临界阻尼状态; 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沿水平方向，oy竖直向上。初始条件：t=0时，x=0,dx/dt=5,y=0,dy/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).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点，写出方程:,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+w02x=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)-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: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); hold on axis(-10 20 -20 20) for i=1:1 t,y=ode23(alzssfun,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.m q=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)