第五講顯式差分和隱式差分

第五講顯式差分和隱式差分第1頁(yè)，共52頁(yè)，2023年，2月20日，星期一回顧1.有限差分法的相容性、穩(wěn)定性和收斂性相容性：針對(duì)差分格式而言，在時(shí)間步長(zhǎng)和空間步長(zhǎng)趨近于零的情況下，如果差分格式的截?cái)嗾`差（差分格式與原有偏微分方程之差）的模趨近于零，則該差分格式與原偏微分方程是相容的，或稱該差分方程與原偏微分方程具有相容性。第2頁(yè)，共52頁(yè)，2023年，2月20日，星期一穩(wěn)定性（stability）：如果偏微分方程的嚴(yán)格解析解有界，差分格式給出的解也有界，稱該差分格式是穩(wěn)定的；如果差分格式給出的解是無(wú)界的，則稱該差分格式是不穩(wěn)定的。穩(wěn)定性反映了差分格式在計(jì)算中控制誤差傳遞的能力第3頁(yè)，共52頁(yè)，2023年，2月20日，星期一收斂性（convergence）：如果當(dāng)時(shí)間和空間步長(zhǎng)趨于零時(shí)，F(xiàn)DE解趨于PDE解，稱該差分格式是收斂的。如果則稱該差分格式是收斂的。收斂性描述的是當(dāng)差分網(wǎng)格無(wú)限細(xì)化時(shí)，差分方程的解是否具有無(wú)限逼近偏微分方程的解的能力Lax等價(jià)定理（Laxequivalencetheorem）：如果逼近一個(gè)給定問(wèn)題的差分格式是相容的，那么該差分格式的收斂性與穩(wěn)定性互為充分必要條件。相容性是比較容易滿足的。在此基礎(chǔ)上，如果滿足了穩(wěn)定性條件，差分格式的收斂性就自動(dòng)滿足。第4頁(yè)，共52頁(yè)，2023年，2月20日，星期一U=0U=0U=100U=02143658710912111413152.5有限差分法實(shí)例(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+1h1h3h2h4第5頁(yè)，共52頁(yè)，2023年，2月20日，星期一(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+1h1h3h2h4forj=2:n-1fori=2:m-1;a((j-1)*m+i,(j-1)*m+i+1)=1;a((j-1)*m+i,(j-1)*m+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內(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)第6頁(yè)，共52頁(yè)，2023年，2月20日，星期一組建A和B矩陣，求解線性方程組得到X第7頁(yè)，共52頁(yè)，2023年，2月20日，星期一%Matlab2Dclear;clc;figure('color','w');

a=zeros(135,135);fori=1:135a(i,i)=1;end;fori=1:7a(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;endfori=1:7a(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)=-0.25;a(15,30)=-0.25;fori=2:14a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i+15)=-0.25;endfori=122:134a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i-15)=-0.25;endfori=1:7forj=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,1);fori=121:135c(i,1)=25;endd=b*c;s=zeros(11,17);fori=2:16s(11,i)=100;endfori=1:9forj=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)第8頁(yè)，共52頁(yè)，2023年，2月20日，星期一第9頁(yè)，共52頁(yè)，2023年，2月20日，星期一2.5應(yīng)用實(shí)例南加州一次未來(lái)大地震的強(qiáng)地面運(yùn)動(dòng)的數(shù)值模擬第10頁(yè)，共52頁(yè)，2023年，2月20日，星期一盆地效應(yīng)第11頁(yè)，共52頁(yè)，2023年，2月20日，星期一Cui,2013第12頁(yè)，共52頁(yè)，2023年，2月20日，星期一Cui,2013第13頁(yè)，共52頁(yè)，2023年，2月20日，星期一Cui,2013第14頁(yè)，共52頁(yè)，2023年，2月20日，星期一Cui,2013第15頁(yè)，共52頁(yè)，2023年，2月20日，星期一總結(jié)：1、有限差分方法給出的數(shù)值解的精度取決于所用的差分形式（向前、向后、中心）。2、偏微分方程的顯式有限差分格式通常是有條件穩(wěn)定的，為了保證得到精確的數(shù)值解，最關(guān)鍵的是需要根據(jù)穩(wěn)定性條件選取正確的空間和時(shí)間步長(zhǎng)。第16頁(yè)，共52頁(yè)，2023年，2月20日，星期一顯式與隱式差分格式主講人：胡才博中國(guó)科學(xué)院大學(xué)地球科學(xué)學(xué)院中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室第17頁(yè)，共52頁(yè)，2023年，2月20日，星期一顯式差分格式（explicitdifferencescheme）

差分方法中可逐層逐點(diǎn)分別求解的格式。特點(diǎn)1.不聯(lián)立解方程；2.時(shí)間步長(zhǎng)和空間步長(zhǎng)的選擇受限制。通常要求時(shí)間步長(zhǎng)足夠小。隱式差分格式（implicitdifferencescheme）特點(diǎn)時(shí)間步長(zhǎng)和空間步長(zhǎng)的選擇不受限制；需要聯(lián)立解方程組顯式和隱式：求解問(wèn)題與時(shí)間相關(guān)第18頁(yè)，共52頁(yè)，2023年，2月20日，星期一例子：1.顯式差分格式：左端：n+1時(shí)刻的值；右端：n時(shí)刻的值。特點(diǎn)：結(jié)構(gòu)簡(jiǎn)潔，直接求解，求解速度快。但是，時(shí)間步長(zhǎng)需滿足：顯式差分格式才能得到穩(wěn)定的數(shù)值解，否則，數(shù)值解將會(huì)不穩(wěn)定而振蕩。第19頁(yè)，共52頁(yè)，2023年，2月20日，星期一顯示差分格式示意圖第20頁(yè)，共52頁(yè)，2023年，2月20日，星期一2.隱式差分格式：時(shí)間一階精度空間二階精度第21頁(yè)，共52頁(yè)，2023年，2月20日，星期一隱式有限差分格式第22頁(yè)，共52頁(yè)，2023年，2月20日，星期一Crank-Nicolson隱式差分格式第23頁(yè)，共52頁(yè)，2023年，2月20日，星期一Crank-Nicolson隱式差分格式第24頁(yè)，共52頁(yè)，2023年，2月20日，星期一Forward-TimeCentral-SpacemethodBackward-Time

Central

-SpacemethodCrank-Nicolson隱式差分格式一般差分格式第25頁(yè)，共52頁(yè)，2023年，2月20日，星期一312546求解區(qū)域：邊界條件：初始條件：一種隱式差分格式的程序?qū)崿F(xiàn)第26頁(yè)，共52頁(yè)，2023年，2月20日，星期一A=sparse(nx,nx);fori=2:nx-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)部邊界第27頁(yè)，共52頁(yè)，2023年，2月20日，星期一312546邊界條件：初始條件：Crank-Nicolson隱式差分格式的程序?qū)崿F(xiàn)第28頁(yè)，共52頁(yè)，2023年，2月20日，星期一A=sparse(nx,nx);fori=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)：B=sparse(nx,nx);fori=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)：第29頁(yè)，共52頁(yè)，2023年，2月20日，星期一例子：牛頓冷卻定律：溫度高于周?chē)h(huán)境的物體向周?chē)劫|(zhì)傳遞熱量逐漸冷卻時(shí)所遵循的規(guī)律。當(dāng)物體表面與周?chē)嬖跍囟炔顣r(shí)，單位時(shí)間從單位面積散失的熱量與溫度差成正比。Tair一階常微分方程的數(shù)值解首先對(duì)時(shí)間和溫度進(jìn)行離散：利用向前差分形式：得到以下的顯式差分格式：Tcap第30頁(yè)，共52頁(yè)，2023年，2月20日，星期一利用向前差分格式：現(xiàn)在改用向后差分形式進(jìn)行近似，得到隱式差分格式：可以驗(yàn)證，當(dāng)時(shí)間步長(zhǎng)趨近于零時(shí)，以上近似解趨于解析解。因此，該格式收斂。穩(wěn)定性條件：T單調(diào)減小的條件顯式差分格式隱式差分格式第31頁(yè)，共52頁(yè)，2023年，2月20日，星期一當(dāng)dt=1.25，tau=0.7時(shí)，顯式差分格式不穩(wěn)定，結(jié)果振蕩；隱式差分格式穩(wěn)定，結(jié)果不精確。隱式差分格式無(wú)條件穩(wěn)定重點(diǎn)考察差分格式的收斂性第32頁(yè)，共52頁(yè)，2023年，2月20日，星期一當(dāng)dt=1，tau=0.7時(shí)，顯式差分格式不穩(wěn)定，結(jié)果振蕩；隱式差分格式穩(wěn)定，結(jié)果不精確。第33頁(yè)，共52頁(yè)，2023年，2月20日，星期一當(dāng)dt=0.5，tau=0.7時(shí)，顯式差分格式穩(wěn)定，隱式差分格式穩(wěn)定，結(jié)果不精確，兩者都不精確。第34頁(yè)，共52頁(yè)，2023年，2月20日，星期一當(dāng)dt=0.1，tau=0.7時(shí)，顯式差分格式穩(wěn)定；隱式差分格式穩(wěn)定；結(jié)果都比較精確。第35頁(yè)，共52頁(yè)，2023年，2月20日，星期一當(dāng)dt=0.01，tau=0.7時(shí)，顯式差分格式穩(wěn)定；隱式差分格式穩(wěn)定；結(jié)果都相當(dāng)精確。第36頁(yè)，共52頁(yè)，2023年，2月20日，星期一第37頁(yè)，共52頁(yè)，2023年，2月20日，星期一當(dāng)dt和tau都大于零時(shí)，該式無(wú)條件滿足，因此混合差分格式無(wú)條件穩(wěn)定。第38頁(yè)，共52頁(yè)，2023年，2月20日，星期一

xt(nt+1)=nt*dt;

plot(xt,T_e,'b.-',xt,T_i,'g.-',xt,T_m,'m.-',xt,T_a,'r.-',);holdon

%

set(gca,'DataAspectRatio',[(max(xt)-min(xt))/(max(T_e)-min(T_e))/311]);xlabel('Time(s)','Fontname','timesnewroman','FontSize',14);

ylabel('Temperature','Fontname','timesnewroman','FontSize',14);title('dt=0.01tau=0.7');%Malab-1Dclear;clc;figure('color','w');

t0=1;%initialtemperaturetau=0.7;%timeconstantdt=0.01;%timeintervalt_total=10;nt=round(t_total/dt);%totaltimesteps

T_e(1)=t0;T_i(1)=t0;T_m(1)=t0;

fori=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);%analyticalresultsend不同差分格式的matlab程序第39頁(yè)，共52頁(yè)，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比第40頁(yè)，共52頁(yè)，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比第41頁(yè)，共52頁(yè)，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比第42頁(yè)，共52頁(yè)，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比第43頁(yè)，共52頁(yè)，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計(jì)算結(jié)果對(duì)比第44頁(yè)，共52頁(yè)，2023年，2月20日，星期一顯式差分格式1.對(duì)步長(zhǎng)有要求；2.無(wú)需解方程第45頁(yè)，共52頁(yè)，2023年，2月20日，星期一二階精度第46頁(yè)，共52頁(yè)，2023年，2月20日，星期一第47頁(yè)，共52頁(yè)，2023年，2月20日，星期一%Matlab2Dclear;clc;figure('color','w');

lx=17;ly=11;%v1=zeros(ly,lx);%forj=2:lx-1v1(ly,j)=100;end%v2=v1;maxt=1;t=0;k=0;while(maxt>1e-6)%k=k+1maxt=0;fori=2:ly-1forj=2:lx-1;v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4;t=abs(v2(i,j)-v1(i,j));if(t>maxt)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=419第48頁(yè)，共52頁(yè)，2023年，2月20日，星期一第49頁(yè)，共52頁(yè)，2023年，2月20日，星期一總結(jié)：顯式格式算法簡(jiǎn)單、易于編程，可以從給定的初始條件開(kāi)始，在時(shí)間上逐層前進(jìn)求解。一些與時(shí)間有關(guān)的偏微分方程的求解，需要用到隱式差分格式，在時(shí)間上計(jì)算數(shù)值解的傳播時(shí)，需要求解線性方程組。通常在計(jì)算的每一個(gè)時(shí)間步，需要求矩陣的逆矩陣。因此，隱式格式算法相對(duì)于顯式格式更復(fù)雜，編程更困難。顯式格式通常比隱式格式的穩(wěn)定性差，如果時(shí)間步長(zhǎng)取得過(guò)大，可能會(huì)給出物理上不正確的結(jié)果。某些隱式格式的優(yōu)點(diǎn)是其無(wú)條件穩(wěn)定性，因此時(shí)間步長(zhǎng)可以取得大一些，但是并不能保證精度很高。可以利用顯式-隱式混合格式（如：Crank-Nicholsonscheme