計算物理課件1-3章.ppt

1、計算物理 第一章引言 計算物理（Computational Physics) 是一門新興的邊緣學(xué)科，是物理學(xué)、數(shù)學(xué)、計算機(jī)科學(xué)三者相結(jié)合的產(chǎn)物。計算物理也是物理學(xué)的一個分支，它與理論物理、實驗物理有密切聯(lián)系，但又保持著自己的相對的獨(dú)立性?？梢赃@樣給計算物理下一個定義：計算物理是以計算機(jī)為工具，以相關(guān)算法為手段，解決復(fù)雜物理問題的一門應(yīng)用科學(xué)。計算物理已經(jīng)給復(fù)雜體系,的物理規(guī)律、物理性質(zhì)的研究提供了重要手段，對物理學(xué)的發(fā)展起著極大的推動作用。 90年代初期，在我國許多大學(xué)的應(yīng)用物理系紛紛開設(shè)了計算物理課程，受到師生們的歡迎。它對訓(xùn)練學(xué)生開闊的思維、培養(yǎng)綜合分析的能力和獲得廣博實用的知識大有益處，

2、同時對理論教學(xué)提供了一個實踐的過程，使得以往抽象的理論教學(xué)形象化。 本課程面向本科生教學(xué)，學(xué)時為42課時，其中需用20,個左右的課時上機(jī)實踐。本課程編程是基于Matlab程序的，需要學(xué)生有良好的Matlab編程能力。 計算物理通常涉及各類線性和非線性方程（組）的求根、數(shù)值積分、常微分方程和偏微分方程的求解、另外還包括付里葉變換、信號處理等內(nèi)容。本課程教學(xué)將涉及微分方程、偏微分方程的求解、付里葉變換和信號處理（主要介紹濾波）。,第二章物理學(xué)中常微分方程及其計算方法 科學(xué)計算中常常需要求解中常微分方程的定解問題，這類問題最簡單的形式是一階方程的初值問題： 雖然求解常微分方程有各種各樣的解析方程，但

3、解析方法只能用來求解一些特殊類型的方程，求解從實際問題中歸結(jié)出來的微分方程主要靠數(shù)值解法。,1、歐拉（Euler)方法 數(shù)值方法的第一步就是將微分方程中的導(dǎo)數(shù)項y進(jìn)行離散化。設(shè)在區(qū)間xn,xn+1的左端點(diǎn)xn,則： y(xn)=f(xn,y(xn) 并用差商 替代導(dǎo)數(shù)項y(xn),則有 或?qū)懗?這就是著名的Euler格式，若初值y0已知，在取定步長h后，就可以逐步疊代算出數(shù)值解y1，y2.。 實際應(yīng)用中Euler格式存在較大的誤差，為此人們又提出了各種改進(jìn)的Euler格式。其中有一種改進(jìn)的Euler格式是：,2、龍格庫塔（Runge-Kutta)方法 2.1 Runge-Kutta方法的設(shè)計思

4、想 差商考察，根據(jù)微分中值定理，存在一點(diǎn) ，因此 設(shè)，稱作區(qū)間xn，xn+1上的平均斜率，這樣，只要對平均斜率提供一種算法，就可以導(dǎo)出相應(yīng)的一種計算,格式。實際是歐拉格式簡單地取點(diǎn)xn處的斜率值f(xn,yn)作為平均斜率K，精度自然很低。改進(jìn)的歐拉格式是用xn與xn+1兩個點(diǎn)的斜率值K1和K2算術(shù)平均作為平均斜率K，而xn+1處的斜率值K2則利用已知信息yn通過Euler格式來預(yù)報。 如果我們設(shè)法在xn,xn+1內(nèi)多報幾個點(diǎn)的斜率值，然后將它們加權(quán)平均作為平均斜率K，則可能構(gòu)造出更高精度的計算格式，這就是Runge-Kutta方法的設(shè)計思想。,2.2 龍格庫塔（RungeKutta）格式 如

5、果y(x)在xn,xn+1上有若干鈄率值，即所謂的預(yù)報鈄率值k1,k2kr以及權(quán)數(shù)a1,a2,.ar則： 現(xiàn)在最常用的是四階RK格式：,這一格式用4個點(diǎn),的斜率值,加權(quán),平均生成平均斜率。通過計算機(jī)語言編程就可以求解這樣的常微分方程。在一個實際工作中，選擇何種計算方法要根據(jù)實際情況而定，基本原則是：(1)滿足需要的精度，(2)節(jié)省機(jī)時。,%微分方程： 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=rk

6、4(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解法 及其在物理學(xué)中的應(yīng)用 Ma

7、tlab有著強(qiáng)大的解常微分方程（組）功能，從方法上來講可分為符號法和數(shù)值法，特別是數(shù)值法應(yīng)用范圍更加廣泛。 3.1 Matlab微分方程符號解法（回顧復(fù)習(xí)） r=dsolve(eq1,eq2,cond1,cond2,.,v) eq1,eq2,.表示微分方程，v為獨(dú)立變量（默認(rèn)的獨(dú)立變量為 t )，cond1, cond2, . 表示初始條件（或者邊界條件）. D表示微分算子（d/dt), D2, D3等表示二次、三次等微分，由 dsolve()求出的只能是解析解，如沒有解析解，需用數(shù)值法解,例: 解微分方程,，初始條件：,y=dsolve(Dy=1+y2,y(0)=1),y = tan(t+1

8、/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),下面討論受阻力作用時振動系統(tǒng)的運(yùn)動特征。比較下面三種情況下振子的軌跡： 1、欠阻尼狀態(tài)； 2、過阻尼狀態(tài)； 3、臨界阻尼狀態(tài)。 彈簧振子系統(tǒng)中質(zhì)點(diǎn)受彈性力以及液體對其的阻力作用 ,根據(jù) 牛頓第二定律,X,是振動系

9、統(tǒng)的固有頻率， 為阻尼因數(shù)，與振動系統(tǒng)以及媒介的性質(zhì)有關(guān),參數(shù)設(shè)置： 、欠阻尼狀態(tài)：,、過阻尼狀態(tài)：,、臨界阻尼狀態(tài)：,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),作業(yè):彈簧振子系統(tǒng)由質(zhì)量m=2kg的滑塊，k=400N/m彈簧組成，已知t0時，m在A處，x0=A=0.1m,由靜止開始釋放，研究：xt, vt

10、, at關(guān)系，作圖表示。,涉及的微分方程和初始條件：,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初值問題的常微分方程的數(shù)值解法 在Matlab、Mathematica等程序設(shè)計軟件

11、出現(xiàn)以前，常微分方程數(shù)值解主要是通過RK算法由Fortran或C語言編程來完成的，工作量大、效率低，特別是對一些較為復(fù)雜的問題求解更是困難。在Matlab 中提供了求解微分方程的函數(shù)odemn (如ode23,ode45等），通過函數(shù)調(diào)用，可以很方便的求解各種類型的線性和非線性常微分方程（或者方程組）。Matlab 可以求解一階常微分方程的初值問題和邊界問題，通過改變方程的,形式，可以求解高階問題。 基本格式x, y=ode45(odefun, xspace, y0, ,p1,p2,.) odefun是與微分方程有關(guān)的函數(shù), xspace是變量取值范圍， y0是初如條件，p1,p2,.是傳遞參

12、數(shù)。 （1）簡單顯式：如 y=f(t, y)，可以通過inline函數(shù)ode函數(shù) 例如：y=y-2x/y,y(0)=1 odefun=inline(y-2*x/y);用inline構(gòu)造函數(shù)odefun x, y=ode45(odefun, 0, 1, 1),(2)復(fù)雜問題（主要指二階以上微分方程以及方程組 ) 此類方程（組）需要首先建立一個關(guān)于一階導(dǎo)數(shù)的函數(shù)，函數(shù)程序由function引導(dǎo)，所以對于二階、三階等高階方程，必須運(yùn)用數(shù)學(xué)變量替換將高階（大于2階）的方程（組）寫成一階微分方程組，這是求解的關(guān)鍵過程。解ode的基本過程如下： 根據(jù)物理的規(guī)律，用微分方程與初始條件進(jìn)行描述 F(y,y,y

13、,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 形成由一階導(dǎo)數(shù)組成的向量：,用一階導(dǎo)數(shù)矩陣建立函數(shù)文件,如odefun.m。 建立一個M文件，將函數(shù)文件與初始條件傳遞給求解器（如ODE45），運(yùn)行后就可得到在指定區(qū)間上的解列向量y. 首先以前面的阻尼振動為例，討論二階常微分方程的數(shù)值解的求解器ode45的使用方法。 微分方程為需將二階轉(zhuǎn)化為一階 構(gòu)建初如條件向量構(gòu)建一階導(dǎo)數(shù)向量：,y(1)=y y(2)=dy(1)/dt dy(2)/dt=-2*

14、b*y(2)-w2*y(1) 初值向量：y0=1,5 一階導(dǎo)數(shù)向量：ydot=y(2); -2*b*y(2)-w*y(1),函數(shù)文件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=欠阻尼狀態(tài); tit2=過阻尼狀態(tài); tit3=臨界阻尼狀態(tài); 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,

15、color,b) end,例1：研究有空氣阻力時拋運(yùn)動的特征。比較下面三種情況下的拋體的軌道：、沒有空氣阻力；、空氣阻力與速度一次方成正比；、空氣阻力與速度二次方成正比。 1、以地面為參考系，以拋點(diǎn)為原點(diǎn)O建立直角坐標(biāo)系OXY,OX沿水平方向，OY豎直向上。初始條件：t=0時，x=0,dx/dt=5,y=0,dy/dt=8,Y,x,質(zhì)點(diǎn)受重力和空氣阻力作用，而空氣阻力包括三種情況。質(zhì)點(diǎn)的運(yùn)動微分方程可表示為：,投影式：,參數(shù)值：b=0,0.1,0.1,p=0,0,1，(b=0,p=0表示無空氣阻力；b=0.1, p=0表示空氣阻力與速度一次方成正比；b=0.1,p=1表示空氣阻力與速度二次方成

16、正比)。,2、用ODE解常微分方程 令 ，將二階微分方程組方程組寫成一階微分方程組,3、程序 、函數(shù)文件： 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(無空氣阻力拋體運(yùn)動,color,b

17、) %加標(biāo)題 hold on %保持圖形窗口繼續(xù)畫圖 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：彈簧小球的非線性運(yùn)動：質(zhì)量為m=0.1k

18、g的小球，掛在原長為l=1.0m，勁度系數(shù)為k=4.8N/m的輕彈簧的一端，彈簧的另一端固定。 1、建立如圖坐標(biāo)。設(shè)小球某時刻t位于p點(diǎn)，寫出方程:,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、程序 、函數(shù)程序： 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(

19、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(彈簧小球的非線性運(yùn)動軌跡,color,b) xlabel(X);ylabel(Y),在具體的實際問題中，常遇到給定初值的常微分方程（組），MATLAB對解決這類問題有著獨(dú)到之處。對于簡單的問題（通常有解析式），我們可以用符號法解，即調(diào)用dsolve函數(shù)。而對于復(fù)雜的問題，有時我們就要花許多時間才能解出最終關(guān)系式，或者我們根本就解不出，此時需要用MATLA

20、B的ODE45函數(shù)方法去解，我們只要花少量的時間編一下程序，不管多難的問題都可以解出來，并且可以用圖像直觀地表示出關(guān)系來 。,課堂練習(xí)：受迫阻尼振動：,初始條件：t=0時,計算區(qū)間0：0.01：100,參數(shù)：,程序：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,y

21、0,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

22、),h); a=a,max(y(:,1); end plot(k,a); xlabel(omega /omega_0) ylabel(振幅),小課題1:,散射，,重核在原點(diǎn)處,初始條件:,%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

23、,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,帶電粒子在電磁場中的運(yùn)動,以場中一點(diǎn)為原點(diǎn),以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);