南京大学计算物理课程课件5

1、2011-3-29E=mcE=mc2 2热传导方程热传导方程 -偏微分方程的数值解法偏微分方程的数值解法( (抛物型抛物型) ) dVtucdtdsnutzyxKdtVS),(对内部无热源物体：对内部无热源物体： dVtucdtdVuKdtVV)(uKtuc)(222222zuyuxutu热传导方程热传导方程Qu(x,y,z)dsdV0)(),()(), 0(0)()0 ,(0, 0),(),(2122ttgtlutgtulxxxuTtxtxuttxu一维热传导方程一维热传导方程l)(1tg)(2tg思路：用差分代替微分思路：用差分代替微分kikikikikikikiuutuhuuuxu,1,

2、2, 1, 1,22|2|空间步长时间步长:h2, 1, 1,1,2huuuuukikikikikiMkkgukguNiihuuhuhuhukNkikikikiki, 1 , 0)(),(1, 2 , 1)()21 (2,1, 00, 12,2, 121,TMhlN,计算步骤：计算步骤：1,2,1, 00. 4)(),()(. 3,. 2,. 1kikNkiukgukguihuMNThl用差分格式计算计算边界值：计算初值：计算给定，xtk+1ki-1 i i+1212h稳定条件稳定条件边界条件差分格式：边界条件差分格式：第一类边界条件第一类边界条件第二类边界条件第二类边界条件第三类边界条件第三

3、类边界条件)(),()(), 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给定计算过程：计算过程：1 . 0, 1, 1hl600161)21(61

4、222hhh1u=0u=0MN,. 2 计算3610MhlN)(),()(. 32,1, 00kgukguihukNki计算边值：计算初值：，10, 2 , 1)1 (4)1 (4)0 ,(0,iihihuxxxui36, 2 , 1 , 0, 00), 1 (), 0(, 0kuututukNk1,. 4kiu用差分格式计算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网格法网格法 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 计算结果：计算结果：二维扩散方程的差分解法二维扩散方程的差分解法砷砷ldlJdtJdJdt线积分变面积分：线积分变面积分：Jt若粒子流仅仅由扩散导致，则若粒子流仅仅由扩散导致，则DJ)(2222yxDt二维扩散方程二维扩散方程建立差分格式：建立差分格式：MijhyNiihxkkt2 , 1 , 02 , 1 , 02 , 1 , 0空间步长时间步长:hM1M2(i, j)MNO对点对点(i, j),在在k时刻有：时刻有：hyhxtkjikjikjikjikjikjikjikjikjikj

7、ikji, 1, 1,2,2, 1, 12,2,1,22)22(, 1, 1, 1, 1,1,hhDkjikjikjikjikjikjikjikji代入扩散方程：代入扩散方程：kjikjikjikjikjikjihDhD, 1, 1, 1, 12,21,)41 (整理得递推公式：整理得递推公式：412hD稳定条件：i-1i+1ikk+1j-1jj+1txNiMjMMMjMMMMjkMikikjNkjNkjkj, 1 , 001, 2 , 1111,2 , 11, 11, 0, 1,21,1, 02211，边界条件：边界条件：流程图流程图开始开始输入输入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计算结果：计算结果： 波动方程波动方程 -偏微分方程的数值解法偏微分方程的数值解法( (双曲双曲型型) )弦线的横振动方程弦线的横振动方程),()(2222txPxyTtyx线密度线密度张力张力外力外力对均匀弦线，无外力的自由振动情况：对均匀弦线，无外力的自由振动情况：22222xyvty为波速其中Tv 一维波动方程一维波动方程)(),()(), 0()()0 ,()()0 ,(0;02122222tgtlytgtyxtxyxxyTtlxxyvty初始条件初始条件边界条件边界条件问题转化为求如下定解问题问题转化为求如下定解问题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初始条件中心差分格式：初始条件中心差分格式：波动方程差分格式波动方程差分格式(k=0)(k=0)一维波动方程定解问题的两种差分格式一维波动方程定解问题的两种差分格式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,计算步骤：计算步骤：1,y. 4. 3/,/. 2,. 1kiTMhlNThlv用差分格式计算计算初值和边值计算给定举例：举例：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流程图流程图开始开始输入输入N,M,V,H初始边界初始边界条件赋值条件赋值K=1,2,MI=1,2,N重新赋边值重新赋边值输出输出jiy,1, 1, 11,kikikikiyyyy计算结果：计算结果： 求解本征值问题求解本征值问题物理问题物理问题 1 ：两端固定的均匀弦自由

15、振动：两端固定的均匀弦自由振动设设代入上述波动方程和边界条件得代入上述波动方程和边界条件得 2000000(0)( )( )ttxxxx ltttua uuux luxux x=0ux=lx)()(),(tTxtxu0)()(0)()0(02tTltTTaT用用 遍除各项即得遍除各项即得 上式成立要求：上式成立要求：即有：即有：Ta22TaT22kTaT. 0) 1()0(222xxkdxdTakdtTd2222和和本征值问题本征值问题. 0) 1()0()0(222xxlxkdxd特点：特点：1 1 只对一些特定的只对一些特定的k k值才能找到满足值才能找到满足 边界条件的解边界条件的解 2

16、 2 方程是线性的和齐次的方程是线性的和齐次的)(:xknn本征函数本征值：xnxnnknnsin)(2 , 1,解析解：数值求解：打靶法数值求解：打靶法1.1.先猜一个试验本征值先猜一个试验本征值2.2.对微分方程作为初值问题求解对微分方程作为初值问题求解3.3.检验所得解是否满足边界条件检验所得解是否满足边界条件4.4.若满足，则该试验本征值为真实本征值，若满足，则该试验本征值为真实本征值，对应的解为本征函数，否则重复对应的解为本征函数，否则重复1,2,31,2,3步步. 0) 1()0()0(222xxlxkdxd对弦振动问题对弦振动问题1. 1. 选取选取k值值2. 2. 从从x=0

17、x=0开始求微分方程，初始条件为开始求微分方程，初始条件为)0( , 0)0(xx3. 3. 求出求出x=lx=l时的时的 值，并判断是否为值，并判断是否为0 04. 4. 重新调整重新调整k k，并再度求解微分方程，直到，并再度求解微分方程，直到x=l x=l 时时 的值为的值为0 0，这时就找到了本征值以及对，这时就找到了本征值以及对 应的本征函数应的本征函数 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计算结果：计算结果：k1=3.141592439138776物理问题物理问题2：一维薛定谔方程的定态解：一维薛定谔方程的定态解ExVdxdm)(2222)(2022222xVEmkkdxdxminxmaxxmxx波函数