南京大學(xué)計(jì)算物理課程課件5

1、2011-3-29E=mcE=mc2 2熱傳導(dǎo)方程熱傳導(dǎo)方程 -偏微分方程的數(shù)值解法偏微分方程的數(shù)值解法( (拋物型拋物型) ) dVtucdtdsnutzyxKdtVS),(對(duì)內(nèi)部無(wú)熱源物體：對(duì)內(nèi)部無(wú)熱源物體： dVtucdtdVuKdtVV)(uKtuc)(222222zuyuxutu熱傳導(dǎo)方程熱傳導(dǎo)方程Qu(x,y,z)dsdV0)(),()(), 0(0)()0 ,(0, 0),(),(2122ttgtlutgtulxxxuTtxtxuttxu一維熱傳導(dǎo)方程一維熱傳導(dǎo)方程l)(1tg)(2tg思路：用差分代替微分思路：用差分代替微分kikikikikikikiuutuhuuuxu,1,

2、2, 1, 1,22|2|空間步長(zhǎng)時(shí)間步長(zhǎng):h2, 1, 1,1,2huuuuukikikikikiMkkgukguNiihuuhuhuhukNkikikikiki, 1 , 0)(),(1, 2 , 1)()21 (2,1, 00, 12,2, 121,TMhlN,計(jì)算步驟：計(jì)算步驟：1,2,1, 00. 4)(),()(. 3,. 2,. 1kikNkiukgukguihuMNThl用差分格式計(jì)算計(jì)算邊界值：計(jì)算初值：計(jì)算給定，xtk+1ki-1 i i+1212h穩(wěn)定條件穩(wěn)定條件邊界條件差分格式：邊界條件差分格式：第一類(lèi)邊界條件第一類(lèi)邊界條件第二類(lèi)邊界條件第二類(lèi)邊界條件第三類(lèi)邊界條件第三

3、類(lèi)邊界條件)(),()(), 0(21tgtlutgtu)(),()(), 0(21tgxtlutgxtu)(),()(),()(), 0()(), 0(2211tgtlutxtlutgtutxtu)()(2,1, 0kgukgukNk)()(2, 1,1, 0, 1kghuukghuukNkNkk)()()()(2,2, 1,1, 01, 0, 1kgukhuukgukhuukNkNkNkkk舉例：舉例：ttutuxxxxutxxutu00), 1 (, 0), 0(10)1 (4)0 ,(0, 1022Thl,. 1給定計(jì)算過(guò)程：計(jì)算過(guò)程：1 . 0, 1, 1hl600161)21(61

4、222hhh1u=0u=0MN,. 2 計(jì)算3610MhlN)(),()(. 32,1, 00kgukguihukNki計(jì)算邊值：計(jì)算初值：，10, 2 , 1)1 (4)1 (4)0 ,(0,iihihuxxxui36, 2 , 1 , 0, 00), 1 (), 0(, 0kuututukNk1,. 4kiu用差分格式計(jì)算kikikikiuuuu, 1, 11,613261xtki-1 i i+10, 20, 30, 41 , 30, 10, 20, 31 , 20, 00, 10, 21 , 1613261613261613261uuuuuuuuuuuu1 , 21 , 31 , 42,

5、 31 , 11 , 21 , 32, 21 , 01 , 11 , 22, 1613261613261613261uuuuuuuuuuuu網(wǎng)格法網(wǎng)格法 DIMENSION U(0:200,0:200) OPEN(10,FILE=CONDUCT1.DAT) A=1/6. H=0.005 DO K=0,200 U(0,K)=0. U(200,K)=0. ENDDO DO I=1,199 U(I,0)=4*I*H*(1-I*H) ENDDO DO K=0,200 DO I=1,199 U(I,K+1)=A*U(I+1,K)+(1-2*A)*U(I,K)+A*U(I-1,K) ENDDO ENDDO

6、 DO I=0,200 WRITE(10,*)I,U(I,200) ENDDO CLOSE(10) END 計(jì)算結(jié)果：計(jì)算結(jié)果：二維擴(kuò)散方程的差分解法二維擴(kuò)散方程的差分解法砷砷ldlJdtJdJdt線(xiàn)積分變面積分：線(xiàn)積分變面積分：Jt若粒子流僅僅由擴(kuò)散導(dǎo)致，則若粒子流僅僅由擴(kuò)散導(dǎo)致，則DJ)(2222yxDt二維擴(kuò)散方程二維擴(kuò)散方程建立差分格式：建立差分格式：MijhyNiihxkkt2 , 1 , 02 , 1 , 02 , 1 , 0空間步長(zhǎng)時(shí)間步長(zhǎng):hM1M2(i, j)MNO對(duì)點(diǎn)對(duì)點(diǎn)(i, j),在在k時(shí)刻有：時(shí)刻有：hyhxtkjikjikjikjikjikjikjikjikjikj

7、ikji, 1, 1,2,2, 1, 12,2,1,22)22(, 1, 1, 1, 1,1,hhDkjikjikjikjikjikjikjikji代入擴(kuò)散方程：代入擴(kuò)散方程：kjikjikjikjikjikjihDhD, 1, 1, 1, 12,21,)41 (整理得遞推公式：整理得遞推公式：412hD穩(wěn)定條件：i-1i+1ikk+1j-1jj+1txNiMjMMMjMMMMjkMikikjNkjNkjkj, 1 , 001, 2 , 1111,2 , 11, 11, 0, 1,21,1, 02211，邊界條件：邊界條件：流程圖流程圖開(kāi)始開(kāi)始輸入輸入N,M,M1,M2,D,H,MT初始邊界初

8、始邊界條件賦值條件賦值S= TD / H2K=1,2,MT重新賦邊值重新賦邊值輸出輸出ji,kjikjikjikjikjikjihDhD, 1, 1, 1,1, 12,21,)41 (I=2,3,M-1J=2,3,M-1 DIMENSION C1(51,71),C2(51,71) OPEN(10,FILE=DIFFUSION.DAT) READ *, N,M,M1,M2,D,H,T,MT N1=N+1 MX=M+1 S=T*D/(H*H) IF(S.GT.0.25)STOP IF(N1.GT.51.OR.MX.GT.71)STOP DO I=1,N1 DO J=1,MXC1(I,J)=0.0

9、ENDDO ENDDO DO I=M1,M2 C1(1,I)=1.0 ENDDODO K=1,MT DO I=2,N DO J=2,M C2(I,J)=(1-4.0*S)*C1(I,J)+S*(C1(I+1),J)+C1(I-1,J)+C1(I,J+1)+C1(I,J-1) ENDDO ENDDO DO J=2,M DO I=2,N C1(I,J)=C2(I,J) ENDDO C1(N1,J)=C2(N,J) ENDDO DO I=2,M1 C1(1,I)=C2(2,I) ENDDO DO I=M2,M C1(1,I)=C2(2,I) ENDDO WRITE(10,*)(I,J,C1(I,J)

10、,I=1,N1),J=1,MX) ENDDOEND計(jì)算結(jié)果：計(jì)算結(jié)果： 波動(dòng)方程波動(dòng)方程 -偏微分方程的數(shù)值解法偏微分方程的數(shù)值解法( (雙曲雙曲型型) )弦線(xiàn)的橫振動(dòng)方程弦線(xiàn)的橫振動(dòng)方程),()(2222txPxyTtyx線(xiàn)密度線(xiàn)密度張力張力外力外力對(duì)均勻弦線(xiàn)，無(wú)外力的自由振動(dòng)情況：對(duì)均勻弦線(xiàn)，無(wú)外力的自由振動(dòng)情況：22222xyvty為波速其中Tv 一維波動(dòng)方程一維波動(dòng)方程)(),()(), 0()()0 ,()()0 ,(0;02122222tgtlytgtyxtxyxxyTtlxxyvty初始條件初始條件邊界條件邊界條件問(wèn)題轉(zhuǎn)化為求如下定解問(wèn)題問(wèn)題轉(zhuǎn)化為求如下定解問(wèn)題xtki-1 i

11、i+122222xyvty思路：用差分代替微分思路：用差分代替微分建立差分格式：建立差分格式：21,1,222, 1, 12222kikikikikikiyyytyhyyyxy差分格式：差分格式：1, 1, 12,21,)()()(1 (2kikikikikiyyyhvyhvy1hv收斂條件：k-1k+1初始條件差分格式：初始條件差分格式：)()0 ,()()0 ,(xtxyxxyNiyytyiii, 1 , 00,1 ,0,)(0 ,1 ,ihyyii向前差分：向前差分：初始條件向前差分格式初始條件向前差分格式Niyytyiii, 1 , 021,1 ,0,2)(1,1 ,ihyyii1,0

12、, 10, 120,21 ,)()()(1 (2iiiiiyyyvhyvhy中心差分：中心差分：)()()(21)(1 (0, 10, 120,21 ,ihyyvhyvhyiiii初始條件中心差分格式：初始條件中心差分格式：波動(dòng)方程差分格式波動(dòng)方程差分格式(k=0)(k=0)一維波動(dòng)方程定解問(wèn)題的兩種差分格式一維波動(dòng)方程定解問(wèn)題的兩種差分格式MkkgyMkkgyNiihihyNiihyMkNiyyyhvyhvykNkiikikikikiki, 1 , 0)(, 1 , 0)(, 1 , 0)()(, 1 , 0)(1, 2 , 1; 1, 2 , 1)()()(1 (22,1, 01 ,0,1

13、, 1, 12,21,MkkgyMkkgyNiihhihihvihhvyNiihyMkNiyyyhvyhvykNkiikikikikiki, 1 , 0)(, 1 , 0)(1, 1 , 0)() 1() 1()(21)()(1 (, 1 , 0)(1, 2 , 1; 1, 2 , 1)()()(1 (22,1, 0221 ,0 ,1, 1, 12,21,計(jì)算步驟：計(jì)算步驟：1,y. 4. 3/,/. 2,. 1kiTMhlNThlv用差分格式計(jì)算計(jì)算初值和邊值計(jì)算給定舉例：舉例：ttytyxxxtxyxxytxxyty00), 1 (), 0(10)1 ()0 ,(sin)0 ,(01022

14、22011hv取2 . 0h2 . 0/vh5/1100hNM取建立差分格式：建立差分格式：, 2 , 104 , 3 , 2 , 1)1 (sin5 , 3 , 2 , 1sin, 2 , 14 , 3 , 2 , 1, 1, 01 ,0,1, 1, 11,kyyiihihihyiihykiyyyykkiikikikiki流程圖流程圖開(kāi)始開(kāi)始輸入輸入N,M,V,H初始邊界初始邊界條件賦值條件賦值K=1,2,MI=1,2,N重新賦邊值重新賦邊值輸出輸出jiy,1, 1, 11,kikikikiyyyy計(jì)算結(jié)果：計(jì)算結(jié)果： 求解本征值問(wèn)題求解本征值問(wèn)題物理問(wèn)題物理問(wèn)題 1 ：兩端固定的均勻弦自由

15、振動(dòng)：兩端固定的均勻弦自由振動(dòng)設(shè)設(shè)代入上述波動(dòng)方程和邊界條件得代入上述波動(dòng)方程和邊界條件得 2000000(0)( )( )ttxxxx ltttua uuux luxux x=0ux=lx)()(),(tTxtxu0)()(0)()0(02tTltTTaT用用 遍除各項(xiàng)即得遍除各項(xiàng)即得 上式成立要求：上式成立要求：即有：即有：Ta22TaT22kTaT. 0) 1()0(222xxkdxdTakdtTd2222和和本征值問(wèn)題本征值問(wèn)題. 0) 1()0()0(222xxlxkdxd特點(diǎn)：特點(diǎn)：1 1 只對(duì)一些特定的只對(duì)一些特定的k k值才能找到滿(mǎn)足值才能找到滿(mǎn)足 邊界條件的解邊界條件的解 2

16、 2 方程是線(xiàn)性的和齊次的方程是線(xiàn)性的和齊次的)(:xknn本征函數(shù)本征值：xnxnnknnsin)(2 , 1,解析解：數(shù)值求解：打靶法數(shù)值求解：打靶法1.1.先猜一個(gè)試驗(yàn)本征值先猜一個(gè)試驗(yàn)本征值2.2.對(duì)微分方程作為初值問(wèn)題求解對(duì)微分方程作為初值問(wèn)題求解3.3.檢驗(yàn)所得解是否滿(mǎn)足邊界條件檢驗(yàn)所得解是否滿(mǎn)足邊界條件4.4.若滿(mǎn)足，則該試驗(yàn)本征值為真實(shí)本征值，若滿(mǎn)足，則該試驗(yàn)本征值為真實(shí)本征值，對(duì)應(yīng)的解為本征函數(shù)，否則重復(fù)對(duì)應(yīng)的解為本征函數(shù)，否則重復(fù)1,2,31,2,3步步. 0) 1()0()0(222xxlxkdxd對(duì)弦振動(dòng)問(wèn)題對(duì)弦振動(dòng)問(wèn)題1. 1. 選取選取k值值2. 2. 從從x=0

17、x=0開(kāi)始求微分方程，初始條件為開(kāi)始求微分方程，初始條件為)0( , 0)0(xx3. 3. 求出求出x=lx=l時(shí)的時(shí)的 值，并判斷是否為值，并判斷是否為0 04. 4. 重新調(diào)整重新調(diào)整k k，并再度求解微分方程，直到，并再度求解微分方程，直到x=l x=l 時(shí)時(shí) 的值為的值為0 0，這時(shí)就找到了本征值以及對(duì)，這時(shí)就找到了本征值以及對(duì) 應(yīng)的本征函數(shù)應(yīng)的本征函數(shù) IMPLICIT NONE REAL*8 L,H,K,TOLK,DK,PHIOLD,PHI INTEGER N COMMON /PHI_PSI/PHI,H,N,K PRINT *, input string length L REA

18、D *, L N=100 H=L/N TOLK=1E-06 K=0 DK=0.5 CALL R_K PHIOLD=PHI DO WHILE(ABS(DK).GT.TOLK) K=K+DK CALL R_K IF(PHI*PHIOLD.GT.0)GOTO 200 K=K-DK DK=DK/2.200 CONTINUE ENDDO PRINT *, EIGEN VALUE =,K END SUBROUTINE R_K IMPLICIT NONE REAL*8 PHI0,PSI0,PHI,PSI,K1, & K2,K3,K4,L1,L2,L3,L4,H,K INTEGER I,N COMMO

19、N /PHI_PSI/PHI,H,N,K PHI0=0 PSI0=0.01 DO I=1,N-1 K1=-K*K*PHI0 L1=PSI0 K2=-K*K*(PHI0+0.5*H*L1) L2=PSI0+0.5*H*K1 K3=-K*K*(PHI0+0.5*H*L2) L3=PSI0+0.5*H*K2 K4=-K*K*(PHI0+H*L3) L4=PSI0+H*K3 PHI=PHI0+H*(K1+2*K2+2*K3+K4)/6. PSI=PSI0+H*(L1+2*L2+2*L3+L4)/6. PHI0=PHI PSI0=PSI ENDDO RETURN END計(jì)算結(jié)果：計(jì)算結(jié)果：k1=3.141592439138776物理問(wèn)題物理問(wèn)題2：一維薛定諤方程的定態(tài)解：一維薛定諤方程的定態(tài)解ExVdxdm)(2222)(2022222xVEmkkdxdxminxmaxxmxx波函數(shù)