东南大学数值分析第七章偏微分方程数值解法.docxVIP

第七章偏微分方程数值解法 Crank-Nicolson 格式 **** （学号） ***** （姓名） 上机题LI要求见教材P346, 10题。 一、算法原理 本文研究下列定解问题（抛物型方程） (0x/,0rT)(0x/)(0rT)du d2u (0x/,0rT) (0x/) (0rT) ot ox u(x,Q) = p(x) ?(0j) = a(0, (lJ) = 0(f) 的有限差分法，其中〃为正常数，fg0为已知函数，且满足边界条件和初始 条件。关于式（1）的求解，采用离散化方法，剖分网格，构造差分格式。其中， 网格剖分是将区域￡ = {0%/,0/7}用两簇平行直线 x = a; =ih （0z M） t = tk = kr （0k N） 分割成矩形网格，其中// = —,r = -分别为空间步长和时间步长。将式（1）中的 M N 偏导数使用不同的差商代替，将得到不同的差分格式，如古典显格式、古典隐格 式、Crank-Nicolson格式等。其中，Crank-Nicolson格式具有更高的收敛阶数， 应用更广泛，故本文采用Crank-Nicolson格式求解抛物型方程。 Crank-Nicolson格式推导：在节点（xprA +-）处考虑式（1）,有 2 毁（也+an z 毁（也+ an z “ 对偏导数竺（兀山+￡）用中心差分展开 ot 2 号(兀‘ 4 + 才)=￡[u(xi‘ L+i)- 3，r*+i)]_ — #(兀,Hi) (L v 4+]) (3) 将炉也+￡）在节点（也）和3,心）表示为 T 1 d2u d1u 卫二二|_左5山）+左Ww） r2 d4u / 人、 ， A 、 飞乔討M）（55） 对以上两个偏导数用二阶差分展开 （X.，4） = * W（也,L） 一 2心,4） + , 4）］ -書养（轨4） 3V ＜和） q2“ j 丽a 心）=乔［”（g，心）一 2心，心）+心,心）］ -￡#（轨心）aW+j -时）_和［（心 -2出 + ） + （将式⑷⑸⑹分别代入式（3）, -时）_和［（心 -2出 + ） + （ i-\ ）］ = f 兀山+ = 用矩阵的形式表示为:1 + r --2r （ r2 2「 M 1M1 用矩阵的形式表示为: 1 + r -- 2 r （ r 2 2 「 M 1 M1 ■ ■ ■ 1-r /? 2 r ~2 l-r r 2 ? ? ? ? ? ? 「 k 1 wi ■ ? ? ? 喘2 ? ? ? r l-r r 心2 Ml 2 2 _4-| _ r l-r 2 曰+却巩匚）+巩心）） Tf 心_2，?+了 ”（兀w-14 +牙卜空（0（4）+ 0仏+1）） Crank-Nicolson格式的截断误差为R = O（t2+h2）,具有较高的精度。 二、计算代码 Crank_Nicolson格式完整代码 functio n U=CranNicolson(f/alxOlxn/dxrtO/tmfdtfg,sO/sn) %采用Crank.Nicolson格式求解抛物线型偏微分方程 % du/dt-a*d2u/dx2=f(x;t) %Input - f抛物方程右端函数 % ?a为二阶导系数 % - xO and xn空间域端庶 % - tO and tm时间域端点 % ?dx为空间步长，dt为时间步长 % -g为初始条件函数 % ?sO,sn为边界条件函数 %Output?U输出,横坐标为空间f纵坐标为时间 M=(tm-tO)/dt;N=(xn-xO)/dx;% 网格数 x=xO+dx:dx:xn-dx; t=tO:dt:tm; u=feval(g/x);u=u,; r=a*dt/dx/dx;% 步长比 %Crank-Nicolson格式：AAuJk+l)=B*u.k+C A=-r/2*[zeros(lrN-l);eye(N-ZN-2)fzeros(N-2J)]-r/2*... [zeros(N2Ueye(N-2,N-2);zeros(：LN-l)H(l+rreye(NJNJ); A=inv(A); B=r/2*[zeros(l/N-l);eye(N-2/N-2)fzeros(N-2fl)]+r/2*... [zeros(N-2/l)feye(N-2/N-2);zeros(lfN-l)]+(l-r)*eye(N-l/N-l); U=D; for k=l:M C=dt*feval(ffx,t(k)+0.5*dt); C=C； C(1)=C ⑴+r/2%feval(s(U(k))+feval(sO,t(k+l))); C(N-l)=C(N-l)+r/2*(feval(snft(k))+teval(sn/t(k+l))); u=A*(B*u+C); U=[U;u]; en