第五講——顯式差分和隱式差分(5)

1、回顧1. 有限差分法基礎(chǔ)2. 差分格式3. 差分方程4. 邊界條件的處理5. 相容性、穩(wěn)定性和收斂性回顧1. 有限差分法的相容性、穩(wěn)定性和收斂性相容性：針對(duì)差分格式而言，在時(shí)間步長(zhǎng)和空間步長(zhǎng)趨近于零的情況下，如果差分格式的截?cái)嗾`差（差分格式與原有偏微分方程之差）的模趨近于零，則該差分格式與原偏微分方程是相容的，或稱該差分方程與原偏微分方程具有相容性。穩(wěn)定性（stability）：如果偏微分方程的嚴(yán)格解析解有界，差分格式給出的解也有界，稱該差分格式是穩(wěn)定的；如果差分格式給出的解是無界的，則稱該差分格式是不穩(wěn)定的。穩(wěn)定性反映了差分格式在計(jì)算中控制誤差傳遞的能力收斂性（convergence）：如果

2、當(dāng)時(shí)間和空間步長(zhǎng)趨于零時(shí)，F(xiàn)DE解趨于PDE解，稱該差分格式是收斂的。如果則稱該差分格式是收斂的。, , ,0,0lim ()0mmi j ki j khtUu 收斂性描述的是當(dāng)差分網(wǎng)格無限細(xì)化時(shí)，差分方程的解是否具有無限逼近偏微分方程的解的能力Lax等價(jià)定理（Lax equivalence theorem）：如果逼近一個(gè)給定問題的差分格式是相容的，那么該差分格式的收斂性與穩(wěn)定性互為充分必要條件。相容性是比較容易滿足的。在此基礎(chǔ)上，如果滿足了穩(wěn)定性條件，差分格式的收斂性就自動(dòng)滿足。U=0U=0U=100U=020U2143658710912111413152.5 有限差分法實(shí)例(i,j)(i+

3、1,j)(i-1,j)01234(i,j)(i+1,j-1)(i-1,j-1)(i,j+1)(i+1,j+1)(i-1,j+1)i-1ii+1j-1jj+1h1h3h2h4(1, )(1, )( ,1)( ,1)4 ( , )0U ijU ijU i jU i jU i j(i,j)(i+1,j)(i-1,j)01234(i,j)(i+1,j-1)(i-1,j-1)(i,j+1)(i+1,j+1)(i-1,j+1)i-1ii+1j-1jj+1h1h3h2h4for j=2:n-1for i=2:m-1; a(j-1)*m+i,(j-1)*m+i+1)=1; a(j-1)*m+i,(j-1)*m

4、+i-1)=1; a(j-1)*m+i,j*m+i)=1; a(j-1)*m+i,(j-2)*m+i)=1; a(j-1)*m+i,(j-1)*m+i)=-4; endend(1, )(1, )( ,1)( ,1)4 ( , )0U ijU ijU i jU i jU i j內(nèi)部節(jié)點(diǎn)：邊界節(jié)點(diǎn)：A矩陣非零系數(shù)減少，同時(shí)引入第一類邊界，方程右端項(xiàng)B向量出現(xiàn)非零元素。局部節(jié)點(diǎn)編號(hào)總體節(jié)點(diǎn)編號(hào)AXB(135,135)AA(135,1)XX(135,1)BB組建A和B矩陣，求解線性方程組得到X%Matlab 2Dclear;clc;figure(color,w); a=zeros(135,135);f

5、or i=1:135 a(i,i)=1;end;for i=1:7 a(15*i+1,15*i+2)=-0.25; a(15*i+1,15*i+16)=-0.25; a(15*i+1,15*i-14)=-0.25;endfor i=1:7 a(15*i+15,15*i+14)=-0.25; a(15*i+15,15*i+30)=-0.25; a(15*i+15,15*i)=-0.25; Enda(1,2)=-0.25;a(1,16)=-0.25;a(121,122)=-0.25;a(121,106)=-0.25;a(135,134)=-0.25;a(135,120)=-0.25;a(15,14

6、)=-0.25;a(15,30)=-0.25;for i=2:14 a(i,i-1)=-0.25; a(i,i+1)=-0.25; a(i,i+15)=-0.25; endfor i=122:134 a(i,i-1)=-0.25; a(i,i+1)=-0.25; a(i,i-15)=-0.25; endfor i=1:7for j=2:14;a(15*i+j,15*i+j-1)=-0.25;a(15*i+j,15*i+j+1)=-0.25;a(15*i+j,15*i+j+15)=-0.25;a(15*i+j,15*i+j-15)=-0.25;endendb=a(-1);c=zeros(135,

7、1);for i=121:135 c(i,1)=25;endd=b*c;s=zeros(11,17);for i=2:16 s(11,i)=100; endfor i=1:9for j=1:15;s(i+1,j+1)=d(15*(i-1)+j,1);endend subplot(1,2,1),mesh(s)axis(0,17,0,11,0,100)subplot(1,2,2),contour(s,32)05101505100204060801005101512345678910112.5 應(yīng)用實(shí)例應(yīng)用實(shí)例南加州一次未來大地震的強(qiáng)地面運(yùn)動(dòng)的數(shù)值模擬盆地效應(yīng)Cui, 2013Cui, 2013Cu

8、i, 2013Cui, 2013總結(jié)：總結(jié)：1、有限差分方法給出的數(shù)值解的精度取決于所用的差分形式（向前、向后、中心）。2、偏微分方程的顯式有限差分格式通常是有條件穩(wěn)定的，為了保證得到精確的數(shù)值解，最關(guān)鍵的是需要根據(jù)穩(wěn)定性條件選取正確的空間和時(shí)間步長(zhǎng)。顯式與隱式差分格式主講人：胡才博中國(guó)科學(xué)院大學(xué)地球科學(xué)學(xué)院中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室 顯式差分格式（explicit difference scheme） 差分方法中可逐層逐點(diǎn)分別求解的格式。 特點(diǎn) 1. 不聯(lián)立解方程； 2.時(shí)間步長(zhǎng)和空間步長(zhǎng)的選擇受限制。通常要求時(shí)間步長(zhǎng)足夠小。隱式差分格式（implicit difference sc

9、heme）特點(diǎn)1. 時(shí)間步長(zhǎng)和空間步長(zhǎng)的選擇不受限制；2. 需要聯(lián)立解方程組顯式和隱式：求解問題與時(shí)間相關(guān)()( ( )Y ttF Y t( ( ), ()0G Y t Y tt例子：1. 顯式差分格式：左端：n+1時(shí)刻的值；右端：n時(shí)刻的值。特點(diǎn)：結(jié)構(gòu)簡(jiǎn)潔，直接求解，求解速度快。但是，時(shí)間步長(zhǎng)需滿足：顯式差分格式才能得到穩(wěn)定的數(shù)值解，否則，數(shù)值解將會(huì)不穩(wěn)定而振蕩。顯示差分格式示意圖2. 隱式差分格式：時(shí)間一階精度空間二階精度隱式有限差分格式111111112(2)(2)2()nnnnnnnniiiiiiiiTTTTTTTTtx111111112(2)(2)2()nnnnnnnniiiiiii

10、iTTTTTTTTtx1111111(22 )(22 )nnnnnniiiiiisTs TsTsTs TsTCrank-Nicolson 隱式差分格式Crank-Nicolson 隱式差分格式1111111(22 )(22 )nnnnnniiiiiisTs TsTsTs TsTForward-Time Central-Space methodBackward -Time Central -Space method1/2Crank-Nicolson 隱式差分格式一般差分格式312546求解區(qū)域：邊界條件：初始條件：一種隱式差分格式的程序?qū)崿F(xiàn)A = sparse(nx,nx);for i=2:nx

11、-1A(i,i-1) = -s;A(i,i ) = (1+2*s);A(i,i+1) = -s;endA(1 ,1 ) = 1;A(nx,nx) = 1;rhs = zeros(nx,1);rhs(2:nx-1) = Told(2:nx-1);rhs(1) = Tleft; rhs(nx) = Tright;內(nèi)部節(jié)點(diǎn)：邊界節(jié)點(diǎn)：載荷項(xiàng)：內(nèi)部邊界312546邊界條件：初始條件：Crank-Nicolson 隱式差分格式的程序?qū)崿F(xiàn)1111111(22 )(22 )nnnnnniiiiiisTs TsTsTs TsT(1)( )nnAcBc100000220000220000220000220000

12、01ssssssAssssss11121(1)3141516nnnnnnnTTTcTTT 10000022000022000022000022000001ssssssBssssss12( )3456nnnnnnnTTTcTTT A = sparse(nx,nx);for i=2:nx-1A(i,i-1) = -s;A(i,i ) = (2+2*s);A(i,i+1) = -s;endA(1 ,1 ) = 1;A(nx,nx) = 1;內(nèi)部節(jié)點(diǎn)：邊界節(jié)點(diǎn)：1111111(22 )(22 )nnnnnniiiiiisTs TsTsTs TsT(1)( )nnAcBc1000002200002200

13、0022000022000001ssssssAssssss10000022000022000022000022000001ssssssBssssss12( )3456nnnnnnnTTTcTTT B= sparse(nx,nx);for i=2:nx-1B(i,i-1) = s;B(i,i ) = (2-2*s);B(i,i+1) = s;endB(1 ,1 ) = 1;B(nx,nx) = 1;內(nèi)部節(jié)點(diǎn)：邊界節(jié)點(diǎn)：例子：牛頓冷卻定律：溫度高于周圍環(huán)境的物體向周圍媒質(zhì)傳遞熱量逐漸冷卻時(shí)所遵循的規(guī)律。當(dāng)物體表面與周圍存在溫度差時(shí)，單位時(shí)間從單位面積散失的熱量與溫度差成正比。Tair一階常微分方程

14、的數(shù)值解( , )dTf T tdt首先對(duì)時(shí)間和溫度進(jìn)行離散：0, ( )jjjttj t TT t 利用向前差分形式：1()jjjt tTTdTOtdtt得到以下的顯式差分格式：1(,)jjjjTTtf T t Tcap利用向前差分格式：1()jjjt tTTdTOtdtt1(,)jjjjTTtf T t 1-=-jjjjttTTTT（1）0-jjtTT （1）現(xiàn)在改用向后差分形式進(jìn)行近似，得到隱式差分格式：1()jjjt tTTdTOtdtt1(,)jjjjTTtf T t 11(1)jjtTT0(1)jjtTT可以驗(yàn)證，當(dāng)時(shí)間步長(zhǎng)趨近于零時(shí)，以上近似解趨于解析解。因此，該格式收斂。穩(wěn)定性

15、條件：T單調(diào)減小的條件1011t( , )Tf T t 顯式差分格式隱式差分格式當(dāng)dt=1.25，tau=0.7時(shí)，顯式差分格式不穩(wěn)定，結(jié)果振蕩；隱式差分格式穩(wěn)定，結(jié)果不精確。11(1)jjtTT隱式差分格式1011t無條件穩(wěn)定重點(diǎn)考察差分格式的收斂性當(dāng)dt=1，tau=0.7時(shí)，顯式差分格式不穩(wěn)定，結(jié)果振蕩；隱式差分格式穩(wěn)定，結(jié)果不精確。當(dāng)dt=0.5，tau=0.7時(shí)，顯式差分格式穩(wěn)定，隱式差分格式穩(wěn)定，結(jié)果不精確，兩者都不精確。當(dāng)dt=0.1，tau=0.7時(shí)，顯式差分格式穩(wěn)定；隱式差分格式穩(wěn)定；結(jié)果都比較精確。當(dāng)dt=0.01，tau=0.7時(shí)，顯式差分格式穩(wěn)定；隱式差分格式穩(wěn)定；結(jié)果

16、都相當(dāng)精確。2231132322322()226 ()24 12jjjjjjjjjjjjTTffdTdtd Tdtd TO dtdtdtdtdtdtdfdtd ffO dtdtdtdTdtdfdt3332() ()jjjTO dtdtdTfO dtdt當(dāng)dt和tau都大于零時(shí)，該式無條件滿足，因此混合差分格式無條件穩(wěn)定。 xt(nt+1)=nt*dt; plot(xt,T_e,b.-, xt,T_i,g.-, xt,T_m,m.-, xt,T_a,r.-,);hold on % set(gca,DataAspectRatio,(max(xt)-min(xt)/(max(T_e)-min(T_e

17、)/3 1 1); xlabel(Time (s),Fontname,times new roman,FontSize,14); ylabel(Temperature,Fontname,times new roman,FontSize,14); title(dt=0.01 tau=0.7);%Malab-1Dclear;clc;figure(color,w); t0=1; % initial temperaturetau=0.7; % time constantdt=0.01; % time intervalt_total=10;nt=round(t_total/dt); % total ti

18、me steps T_e(1)=t0; T_i(1)=t0; T_m(1)=t0; for i=1:nt; xt(i)=(i-1)*dt;T_e(i+1)=T_e(i)*(1-dt/tau); %explicitT_i(i+1)=T_i(i)/(1+dt/tau); % implicitT_m(i+1)=T_m(i)*(1-dt/2/tau)/(1+dt/2/tau); %mix T_a(i)=t0*exp(-xt(i)/tau); %analytical resultsend不同差分格式的matlab程序混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比混合差分格式精度最高！不同差分格式計(jì)算結(jié)果

19、對(duì)比混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比顯式差分格式1.對(duì)步長(zhǎng)有要求；2.無需解方程二階精度%Matlab 2Dclear;clc;figure(color,w); lx=17;ly=11; %v1=zeros(ly,lx); %for j=2:lx-1 v1(ly,j)=100;end %v2=v1;maxt=1;t=0;k=0;while(maxt1e-6) %k=k+1maxt=0;for i=2:ly-1for j=2:lx-1;v2(i,j)=(v1(i,j+1)+v1(i+1,j)

20、+v1(i-1,j)+v1(i,j-1)/4; t=abs(v2(i,j)-v1(i,j);if(tmaxt) maxt=t;endendendv1=v2;end %subplot(1,2,1),mesh(v1)axis(0,17,0,11,0,100)subplot(1,2,2),contour(v1,32)迭代解法:K:迭代步數(shù)K=4190510150510020406080100510151234567891011 總結(jié)總結(jié)：顯式格式算法簡(jiǎn)單、易于編程，可以從給定的初始條件開始，在時(shí)間上逐層前進(jìn)求解。一些與時(shí)間有關(guān)的偏微分方程的求解，需要用到隱式差分格式，在時(shí)間上計(jì)算數(shù)值解的傳播時(shí)，需要

21、求解線性方程組。通常在計(jì)算的每一個(gè)時(shí)間步，需要求矩陣的逆矩陣。因此，隱式格式算法相對(duì)于顯式格式更復(fù)雜，編程更困難。顯式格式通常比隱式格式的穩(wěn)定性差，如果時(shí)間步長(zhǎng)取得過大，可能會(huì)給出物理上不正確的結(jié)果。某些隱式格式的優(yōu)點(diǎn)是其無條件穩(wěn)定性，因此時(shí)間步長(zhǎng)可以取得大一些，但是并不能保證精度很高?？梢岳蔑@式-隱式混合格式（如： Crank-Nicholson scheme ），它們無條件穩(wěn)定，而且精度高于相應(yīng)的顯式和隱式格式。 總結(jié)有限差分法的主要內(nèi)容包括：1、根據(jù)問題的特點(diǎn)將定解區(qū)域作網(wǎng)格剖分；在所有網(wǎng)格節(jié)點(diǎn)上用有限差分格式對(duì)導(dǎo)數(shù)求近似，對(duì)函數(shù)、初始和邊界條件求近似；把原方程離散化為代數(shù)方程組，即有

22、限差分方程組.2、有限差分模型形態(tài)的理論分析，以保證計(jì)算過程可行及計(jì)算結(jié)果正確：解的相容性：對(duì)于一個(gè)微分方程建立的各種差分格式，為了有實(shí)用意義，一個(gè)基本要求是它們能夠任意逼近微分方程，這就是相容性要求。解的穩(wěn)定性：因?yàn)椴罘指袷降挠?jì)算過程是逐層推進(jìn)的，在計(jì)算第m1層的近似值時(shí)要用到第m層的近似值，直到與初始值有關(guān)。前面各層若有舍入誤差，必然影響到后面各層的值，如果誤差的影響越來越大，以致差分格式的精確解的面貌完全被掩蓋，這種格式是不穩(wěn)定的，相反如果誤差的傳播是可以控制的，就認(rèn)為格式是穩(wěn)定的。只有在這種情形下，差分格式在實(shí)際計(jì)算中的近似解才可能任意逼近差分方程的精確解。解的收斂性：一個(gè)差分格式是否有用，最終要看