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

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

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

3、解析方法只能用來求解一些特殊類型的方程，求解從實際問題中歸結出來的微分方程主要靠數值解法。,1、歐拉（Euler)方法 數值方法的第一步就是將微分方程中的導數項y進行離散化。設在區(qū)間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、想 差商考察，根據微分中值定理，存在一點 ，因此 設，稱作區(qū)間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）格式 如

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

7、tlab有著強大的解常微分方程（組）功能，從方法上來講可分為符號法和數值法，特別是數值法應用范圍更加廣泛。 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

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

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

10、, 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等程序設計軟件

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

12、數。 （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的基本過程如下： 根據物理的規(guī)律，用微分方程與初始條件進行描述 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 形成由一階導數組成的向量：,用一階導數矩陣建立函數文件,如odefun.m。 建立一個M文件，將函數文件與初始條件傳遞給求解器（如ODE45），運行后就可得到在指定區(qū)間上的解列向量y. 首先以前面的阻尼振動為例，討論二階常微分方程的數值解的求解器ode45的使用方法。 微分方程為需將二階轉化為一階 構建初如條件向量構建一階導數向量：,y(1)=y y(2)=dy(1)/dt dy(2)/dt=-2*

14、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=欠阻尼狀態(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：研究有空氣阻力時拋運動的特征。比較下面三種情況下的拋體的軌道：、沒有空氣阻力；、空氣阻力與速度一次方成正比；、空氣阻力與速度二次方成正比。 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表示空氣阻力與速度二次方成

16、正比)。,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

17、) %加標題 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：彈簧小球的非線性運動：質量為m=0.1k

18、g的小球，掛在原長為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(

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

20、B的ODE45函數方法去解，我們只要花少量的時間編一下程序，不管多難的問題都可以解出來，并且可以用圖像直觀地表示出關系來 。,課堂練習：受迫阻尼振動：,初始條件：t=0時,計算區(qū)間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,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:,散射，,重核在原點處,初始條件:,%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,帶電粒子在電磁場中的運動,以場中一點為原點,以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);