東南大學(xué)數(shù)值分析偏微分方程數(shù)值解法

1、東南大學(xué)《數(shù)值分析》上機(jī)練習(xí)——算法與程序設(shè)計(jì)實(shí)驗(yàn)報(bào)告第七章第七章偏微分方程數(shù)值解法偏微分方程數(shù)值解法——CrankNicolson格式格式（學(xué)號(hào)）（姓名）上機(jī)題目要求見教材P34610題。一、算法原理一、算法原理本文研究下列定解問題(拋物型方程)(1)22()(00)(0)()(0)(0)()(1)()(0)uuafxtxltTtxuxxxluttutttT???????????????????????????的有限差分法，其中為正常

2、數(shù)，為已知函數(shù)，且滿足邊界條件和初af???始條件。關(guān)于式(1)的求解，采用離散化方法，剖分網(wǎng)格，構(gòu)造差分格式。其中，網(wǎng)格剖分是將區(qū)域用兩簇平行直線??00DxltT?????(0)(0)ikxxihiMttkkN????????????分割成矩形網(wǎng)格，其中分別為空間步長和時(shí)間步長。將式(1)中的lThMN???偏導(dǎo)數(shù)使用不同的差商代替，將得到不同的差分格式，如古典顯格式、古典隱格式、CrankNicolson格式等。其中，CrankN

3、icolson格式具有更高的收斂階數(shù)，應(yīng)用更廣泛，故本文采用CrankNicolson格式求解拋物型方程。CrankNicolson格式推導(dǎo)：格式推導(dǎo)：在節(jié)點(diǎn)處考慮式(1)，有()2ikxt??(2)22()()()222ikikikuuxtaxtfxttx????????????對(duì)偏導(dǎo)數(shù)用中心差分展開()2ikuxtt????(3)??2311131()()()()()224kkikikikiikikuuxtuxtuxtxtttt??

4、???????????????將在節(jié)點(diǎn)和表示為22()2ikuxtx????()ikxt1()ikxt?東南大學(xué)《數(shù)值分析》上機(jī)練習(xí)——算法與程序設(shè)計(jì)實(shí)驗(yàn)報(bào)告二、計(jì)算代碼二、計(jì)算代碼Crank_Nicolson格式完整代碼functionU=Crank_Nicolson(fax0xndxt0tmdtgs0sn)%采用Crank_Nicolson格式求解拋物線型偏微分方程%dudtad2udx2=f(xt)%Inputf拋物方程右端函數(shù)%

5、a為二階導(dǎo)系數(shù)%x0xn空間域端點(diǎn)%t0tm時(shí)間域端點(diǎn)%dx為空間步長，dt為時(shí)間步長%g為初始條件函數(shù)%s0sn為邊界條件函數(shù)%OutputU輸出橫坐標(biāo)為空間，縱坐標(biāo)為時(shí)間M=(tmt0)dtN=(xnx0)dx%網(wǎng)格數(shù)x=x0dx:dx:xndxt=t0:dt:tmu=feval(gx)u=ur=adtdxdx%步長比%CrankNicolson格式：Au_(k1)=Bu_kCA=r2[zeros(1N1)eye(N2N2)zero

6、s(N21)]r2...[zeros(N21)eye(N2N2)zeros(1N1)](1r)eye(N1N1)A=inv(A)B=r2[zeros(1N1)eye(N2N2)zeros(N21)]r2...[zeros(N21)eye(N2N2)zeros(1N1)](1r)eye(N1N1)U=[]fk=1:MC=dtfeval(fxt(k)0.5dt)C=CC(1)=C(1)r2(feval(s0t(k))feval(s0t(k1