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

第五講顯式差分和隱式差分第1頁，共52頁，2023年，2月20日，星期一回顧1.有限差分法的相容性、穩(wěn)定性和收斂性相容性：針對差分格式而言，在時間步長和空間步長趨近于零的情況下，如果差分格式的截斷誤差（差分格式與原有偏微分方程之差）的模趨近于零，則該差分格式與原偏微分方程是相容的，或稱該差分方程與原偏微分方程具有相容性。第2頁，共52頁，2023年，2月20日，星期一穩(wěn)定性（stability）：如果偏微分方程的嚴格解析解有界，差分格式給出的解也有界，稱該差分格式是穩(wěn)定的；如果差分格式給出的解是無界的，則稱該差分格式是不穩(wěn)定的。穩(wěn)定性反映了差分格式在計算中控制誤差傳遞的能力第3頁，共52頁，2023年，2月20日，星期一收斂性（convergence）：如果當時間和空間步長趨于零時，F(xiàn)DE解趨于PDE解，稱該差分格式是收斂的。如果則稱該差分格式是收斂的。收斂性描述的是當差分網(wǎng)格無限細化時，差分方程的解是否具有無限逼近偏微分方程的解的能力Lax等價定理（Laxequivalencetheorem）：如果逼近一個給定問題的差分格式是相容的，那么該差分格式的收斂性與穩(wěn)定性互為充分必要條件。相容性是比較容易滿足的。在此基礎上，如果滿足了穩(wěn)定性條件，差分格式的收斂性就自動滿足。第4頁，共52頁，2023年，2月20日，星期一U=0U=0U=100U=02143658710912111413152.5有限差分法實例(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頁，共52頁，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é)點：邊界節(jié)點：A矩陣非零系數(shù)減少，同時引入第一類邊界，方程右端項B向量出現(xiàn)非零元素。局部節(jié)點編號總體節(jié)點編號第6頁，共52頁，2023年，2月20日，星期一組建A和B矩陣，求解線性方程組得到X第7頁，共52頁，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頁，共52頁，2023年，2月20日，星期一第9頁，共52頁，2023年，2月20日，星期一2.5應用實例南加州一次未來大地震的強地面運動的數(shù)值模擬第10頁，共52頁，2023年，2月20日，星期一盆地效應第11頁，共52頁，2023年，2月20日，星期一Cui,2013第12頁，共52頁，2023年，2月20日，星期一Cui,2013第13頁，共52頁，2023年，2月20日，星期一Cui,2013第14頁，共52頁，2023年，2月20日，星期一Cui,2013第15頁，共52頁，2023年，2月20日，星期一總結(jié)：1、有限差分方法給出的數(shù)值解的精度取決于所用的差分形式（向前、向后、中心）。2、偏微分方程的顯式有限差分格式通常是有條件穩(wěn)定的，為了保證得到精確的數(shù)值解，最關鍵的是需要根據(jù)穩(wěn)定性條件選取正確的空間和時間步長。第16頁，共52頁，2023年，2月20日，星期一顯式與隱式差分格式主講人：胡才博中國科學院大學地球科學學院中國科學院計算地球動力學重點實驗室第17頁，共52頁，2023年，2月20日，星期一顯式差分格式（explicitdifferencescheme）

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

Central

-SpacemethodCrank-Nicolson隱式差分格式一般差分格式第25頁，共52頁，2023年，2月20日，星期一312546求解區(qū)域：邊界條件：初始條件：一種隱式差分格式的程序?qū)崿F(xiàn)第26頁，共52頁，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é)點：邊界節(jié)點：載荷項：內(nèi)部邊界第27頁，共52頁，2023年，2月20日，星期一312546邊界條件：初始條件：Crank-Nicolson隱式差分格式的程序?qū)崿F(xiàn)第28頁，共52頁，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é)點：邊界節(jié)點：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é)點：邊界節(jié)點：第29頁，共52頁，2023年，2月20日，星期一例子：牛頓冷卻定律：溫度高于周圍環(huán)境的物體向周圍媒質(zhì)傳遞熱量逐漸冷卻時所遵循的規(guī)律。當物體表面與周圍存在溫度差時，單位時間從單位面積散失的熱量與溫度差成正比。Tair一階常微分方程的數(shù)值解首先對時間和溫度進行離散：利用向前差分形式：得到以下的顯式差分格式：Tcap第30頁，共52頁，2023年，2月20日，星期一利用向前差分格式：現(xiàn)在改用向后差分形式進行近似，得到隱式差分格式：可以驗證，當時間步長趨近于零時，以上近似解趨于解析解。因此，該格式收斂。穩(wěn)定性條件：T單調(diào)減小的條件顯式差分格式隱式差分格式第31頁，共52頁，2023年，2月20日，星期一當dt=1.25，tau=0.7時，顯式差分格式不穩(wěn)定，結(jié)果振蕩；隱式差分格式穩(wěn)定，結(jié)果不精確。隱式差分格式無條件穩(wěn)定重點考察差分格式的收斂性第32頁，共52頁，2023年，2月20日，星期一當dt=1，tau=0.7時，顯式差分格式不穩(wěn)定，結(jié)果振蕩；隱式差分格式穩(wěn)定，結(jié)果不精確。第33頁，共52頁，2023年，2月20日，星期一當dt=0.5，tau=0.7時，顯式差分格式穩(wěn)定，隱式差分格式穩(wěn)定，結(jié)果不精確，兩者都不精確。第34頁，共52頁，2023年，2月20日，星期一當dt=0.1，tau=0.7時，顯式差分格式穩(wěn)定；隱式差分格式穩(wěn)定；結(jié)果都比較精確。第35頁，共52頁，2023年，2月20日，星期一當dt=0.01，tau=0.7時，顯式差分格式穩(wěn)定；隱式差分格式穩(wěn)定；結(jié)果都相當精確。第36頁，共52頁，2023年，2月20日，星期一第37頁，共52頁，2023年，2月20日，星期一當dt和tau都大于零時，該式無條件滿足，因此混合差分格式無條件穩(wěn)定。第38頁，共52頁，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頁，共52頁，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計算結(jié)果對比第40頁，共52頁，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計算結(jié)果對比第41頁，共52頁，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計算結(jié)果對比第42頁，共52頁，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計算結(jié)果對比第43頁，共52頁，2023年，2月20日，星期一混合差分格式精度最高！不同差分格式計算結(jié)果對比第44頁，共52頁，2023年，2月20日，星期一顯式差分格式1.對步長有要求；2.無需解方程第45頁，共52頁，2023年，2月20日，星期一二階精度第46頁，共52頁，2023年，2月20日，星期一第47頁，共52頁，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頁，共52頁，2023年，2月20日，星期一第49頁，共52頁，2023年，2月20日，星期一總結(jié)：顯式格式算法簡單、易于編程，可以從給定的初始條件開始，在時間上逐層前進求解。一些與時間有關的偏微分方程的求解，需要用到隱式差分格式，在時間上計算數(shù)值解的傳播時，需要求解線性方程組。通常在計算的每一個時間步，需要求矩陣的逆矩陣。因此，隱式格式算法相對于顯式格式更復雜，編程更困難。顯式格式通常比隱式格式的穩(wěn)定性差，如果時間步長取得過大，可能會給出物理上不正確的結(jié)果。某些隱式格式的優(yōu)點是其無條件穩(wěn)定性，因此時間步長可以取得大一些，但是并不能保證精度很高?？梢岳蔑@式-隱式混合格式（如：Crank-Nicholsonscheme