第四章 抛物型微分方程有限差分法1设已知初边值问题22, 01, 0<(,0)sin , 01(0,)(1,)0, 0 u ux t t x u x x x u t u t t T π⎧∂∂=<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩T ≤, 试用最简显格式求上述问题的数值解。

取h=0.1,r=0.1.0 1/10 2/10 … 1 T 2τ τt解: 1.矩形网格剖分区域. 取空间步长1, 时间2510h =0.00τ=以及0.01τ=的矩形网格剖分区域, 用节点)表示坐标点(,j k (,)(,)j k x t jh k τ=, 0,1,...1/; 0,1,...,/j h k T τ==, 如图所示.显然, 我们需要求解这(1/1)(/1)h T τ+×+个点对应的函数值. 事实上由已知初边界条件蓝标附近的点可直接得到, 所以只要确定微分方程的解在其它点上的取值即可. 沿用记号[]k(,)j j k u x t =。

u 2. 建立差分格式, 对于11,...1; 0,1,...,1Tj k hτ=−=−, 用向前差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式:1122k k k k k1jj j j u u u u u h ++−+=. 变形j τ−−有:1112(12) (k k k kj j j j u ru r u ru r h τ+−+=+−+=(4.1)用向后差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式最简隐格式:111122k k k k k j jj j j u u u u u h τ++++−−+=11+−1kj +,变形有:1111(12) k k k j j j ru r u ru u ++−−−++−= (4.2)(4.1)*0.5+(4.2)*0.5得CN 格式为:111112222k k k k k k k k j jj j j j j j u u u u u u u u h τ+++−+−−++−+=111++−1kj +x x变形有:111111(22)(22) k k k k k j j j j j ru r u ru ru r u ru ++−−+−−++−=+−+ (4.3)3 初边界点差分格式处理.对于初始条件u x (,0)sin , 01=π≤≤h 离散为(4.4)0sin 0,1,...1/j u jh j π==对于边界条件离散为(0,)(1,)0, 0 u t u t t T ==≤≤00 0,1,.../k k N u u k T τ===(4.5)总结: 联立方程(4.1)(4.4)(4.5)得到已知问题的最简显格式差分方程组:11100(12)1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N u ru r u ru T j k h u jh j h u u k T τπτ+−+⎧=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.2)( 4.4)( 4.5)得到已知问题的最简隐格式差分方程组:1111100(12) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N ru r u ru u T j k h u jh j h u u k T τπτ++−−+⎧−++−=⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.3)( 4.4)( 4.5)得到已知问题的CN 格式差分方程组:11111100(22)(22) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k k j j j j j jk k N ru r u ru ru r u ru T j k h u jh j h u u k T τπτ++−−+−⎧−++−=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩1k j + 4 求解并显示结果利用软件计算(Matlab)如上最简显格式差分方程组.h=1/10;tau=0.0025;T=0.5; r=tau/h^2;M=1/h+1;N=T/tau+1; u=zeros(M,N);for m=1:Mu(m,1)=sin((m-1)*h*pi); endu(1,1:N)=0;u(M,1:N)=0;for n=1:N-1for m=2:M-1u(m,n+1)=r*(u(m+1,n)+u(m-1,n))+(1-2*r)*u(m,n); end end u=u’ 这样我们就计算出不同时刻不同位置k t j x 对应的函数值(,)j k u x t 取tau=0.0025, 即r=0.25绘图, 取tau=0.01, r=1再绘图,如图()图4.2 习题1数值解图示(左r=0.25, 右r=1)2.试构造初边值问题 ()()()()(), 0.51, 0,,0, 0.51,0.5,0, 1,0.51,, 0u u x x x T t x x u x x x u ⎪∂u t t u t t T x ϕ⎧∂∂∂⎛⎞=<<<≤⎜⎟⎪∂∂∂⎝⎠⎪⎪=≤≤⎨⎪==−≤≤⎪∂⎩的显格式，并给出其按最大范数稳定的充分条件。

代替偏导数建立显格式, 也可以利用积分插对区域进行矩形网格剖分, 解: 此即可以利用差商值法. 下面后一种方法建立差分格式. 0.5, j k x jh t k τ=+=, 0,1,...1/:, j h N =N 0,1,...,/:k T ττ==.= 假设(j , k )为内点，固定t , 取关于x 的对偶区间22j j h hx x x −≤≤+, 在此区间上对方程两边进行积分得：112211j j xx dx x t x −−=⎜∂∂∂j j xx uu dx x ++∂∂∂⎛⎞⎟⎝⎠ (4.6)对于等式左侧积分利用中矩形公式, 有∫∫ 1213u ∂⎡⎤()j j xxju dx h O h t t +−∂=+⎢∂∂⎣⎦∫, 用向前差商代替关于时间的一阶偏导数, 有1[][]()k k kj jju u u O t ττ+−∂⎡⎤=+⎢⎥∂⎣⎦, 即得121213[][]()j j k k x j j xu u u dx h O h t τ+−+−∂=+∂∫(4.7) 6)右端直接积分, 有对于(4.121211112222j j xxj j x −⎢∂⎣⎦j j u u u x dx x x x x x +−++−⎛⎞∂∂∂∂⎛⎞⎡⎤⎡⎤⎜⎟=−⎜⎟⎢⎥⎜⎟∂∂∂⎝⎠⎣⎦⎝⎠∫对上式右端, 用一阶中心差商代替关于空间的一阶偏导数, 得11221122[][][][](), ()k k k k j j j j j j u u u u u u O h O h x h x +−+−−−h ∂∂⎡⎤⎡⎤=+=+⎢⎥⎢∂∂⎣⎦⎣⎦ (4.8) 联立(4.6)-(4.10), 得()()1[][]k k u u 211112221[][][][]()j jk k k kj j j j j j x u u x u u O h h ττ+−+−=−−−++⎢⎥⎣⎦ 舍去误差项+−⎡⎤[]kk, 用j j u u 代替, 那么可以定义差分格式()()111112221, 1,,1; 0,,1k k j jk k k k j j j j j j u u x u u x u u j N k N h ττ++−+−−⎡⎤=−−−=−=−⎢⎥⎣⎦"". (4.9) x x ϕ=对于初始条件() 0.51,u x ≤≤和左边界条件(),0,()0.5u ,0, 0t t T =≤≤可直接定义差分格式;N N ()0,0,1,,j j u x j ϕ=="00, 0,1,,ku k τ=="此时应有成立. 而对于右边界条件()0()0N x x τϕϕ==()()1,0.51,, 0ut u t t x∂T =−≤∂≤, 利用增设虚点法,一方面用一阶中心差商代替一阶偏导数, 即211[][]()k kN N Nu u u O h x h +−−∂⎡⎤=+⎢⎥∂⎣⎦, 由此定义差分方程110.5, 0,1,,k kk N N N u u u k N hτ+−−=−=" (4.10) 另一方面把看作内点由(4.9)得()1,t ()()111112,N N u u −−221k kk k k kNN N N N N u u x u u x h τ+++−⎡⎤−=−−⎢⎥⎣⎦.11)边界条件的离散方程为(4.11)将(4.10)和(4联立, 消去1kN u +, 得右()()1, 0,1,,1k N u k N τ−=−" 111122210.5k kk k k kNN N N N NN N u u x u h u u x u h τ+−+−⎡⎤−=−−−−⎢⎥⎣⎦到此建立的最简显格式为()()()()()11111222011111222, 1,,1; 0,, 1 , 0,1,,;0, 1,,10.5, 0,1,,1k k j k kj j j j j j j j kk kk k k k k N N N N N N N N N u u j N k N h u x j N u k N u u x u h u u x u u k N h ττττϕτ++−+−+−−+−⎤−=−=−⎪⎣⎦⎪==⎪⎨==⎪⎪⎡⎤−=−−−−=−⎢⎥⎣⎦⎩"""""⎪⎪1j k k u x u u x u ⎧−⎡=−−⎪⎢⎥取2r h τ=, 则(4.9)等价于()111112212k k kk j j j j j j u rx u rx u rxu +−+−+=+−+j (4.12)注意到0.51x <<, 所以11220, 0, 0j j j x xx−+>>. >如果120j rx −≥这时(4.12)右端各项系数非负，则有：()11111k j +⎟⎠2212k k kj j j j j j u rx u r x u x u +−−+⎛⎞≤−++⎜⎝()112212kk j cck j j rx u r x u xu −+⎛⎞≤−++⎜⎟⎝⎠c()111222kk j j j c crx u r x h x h u ⎛⎞=−+−++⎜⎟⎝⎠()122kkj j ccrx u rx u =−+k cu = 所以，当时，按最大范数稳定. 此时120j rx −≥11min22jj r x ≤≤是其按最大范数稳定的充分条件. 3. 试构造初边值问题:221,01,0,1t x x x (,0)(),01,(0,)0,(1,)0,0u u ux t T u x x x ut u t t T x ϕ⎧∂∂∂=+<<<≤⎪=≤≤⎨⎪∂⎪==≤≤∂⎪⎩的显格式,并给出其按最大范数稳定的充解：对区域进行矩形网格剖分, ∂+∂∂⎪⎪分条件., x jh t k τ==, 0,1,...1/, j h = 0,1,...,/k T τ=. 在内点(j , k )处，用向前差于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 用一阶中心差商代替关于空间的一阶偏导数, 有数值微分公式商代替关1[][][]O()k kj j k j u u u t ττ+−∂=+∂, 22222[][]O()j x j u h 1k ku x h δ=+∂, ∂11[][]k k j j k u u u []O()2j h x h=+∂, 其中+−−∂211[][]2[][]k k k k x j j j u u u u δ+−=−+. j 将其代入方程, 则有111222[][][][]11[]()12k k k kj jj j k x j u u u u u O h jh hδττ++−−−≡++++h 去掉误差项, 用那么可以定义差分格式[]k kj j u u 代替, 1112211, 1,...1/1, 0,1,...,/112h jh h hk k k kj jj j k x j u u u u u j k T δτ++−−−≡+=−=−,对于初始条件τ+(,0)(),01,u x x x ϕ=≤≤和右边界条件(1,)0,0u t t T =≤≤可直接定义差分格式;h ()0,0,1,,1/j j u x j ϕ=="0, 0,1,,1/kN u k τ=="此时应有()0N x τϕ=成立. 而对于左边界条件(0,)0, 0u∂t t T x=≤≤∂, 利用增设虚点法,一方面用一阶中心差商代替一阶偏导数, 即211⎣⎦[][]()k k u u u O h x h −−∂⎡⎤=+⎢⎥∂, 由此定义差分方程11[][],0,1,,k ku u k N τ−==" (4.10) 另一方面把看作内点由(4.9)得()1,t ()()11122N N N x u u x h τ++⎢⎥⎣1121,k kk k k k NN N N N u u u u +−−⎡⎤−=−−−⎦(4.11) 将(4.10)和(4.11)联立, 消去得右边界条件的离散方程为1kN u +,()()115k kk k k k kNN u u 11112220., 0,1,,N N N N N N N x u h u u x u u k N h ττ−−+−−=⎢⎥⎣⎦"到此建立的最简显格式为+⎡⎤−=−−−()()()()()111112220111112221, 1,,1; 0,, 1 , 0,1,,;0, 0,1,,10.5, 0,1,,k k j j k k k kj j j j j j j j kk kk k k k k N N N N N N N N N u u x u u x u u j N k N h u x j N u k N u u x u h u u x u u k N h ττττϕτ++−+−+−−+−⎧−⎡⎤=−−−=−=⎪⎢⎥⎪⎣⎦⎪==⎪⎨==⎪⎪⎡⎤−=−−−−=⎢⎥⎣⎦⎩"""""⎪⎪−2r h τ=, 那么(4.13)等价于()记()()11112112121k kk k j j j h h u r j u r u jh jh ++−⎞⎛⎞=+−⎜⎟⎜⎟⎜⎟⎜⎟++⎝⎠⎝⎠(4.14)下面讨论格式的稳定性. 如果u r ⎛−++12r ≤，则有 120r −≥, ()1021hr ⎛−>, jh ⎞⎜⎟⎜⎟+⎝⎠故()()()11112112121k kk jjj j j h h ur u r u r x x +k j u +−⎛⎞⎛⎞⎜⎟⎜⎟≤−+++−⎜⎟⎜⎟++⎝⎠⎝⎠()()()12112121kk ccj j hhr u r u r u x x ⎛⎞⎛⎞⎜⎟⎜⎟≤−+++−⎜⎟⎜⎟++⎝⎠⎝⎠k c kcu =()01j N ≤≤−1k k c c u u +≤由此可得，即此时最简显格式按最大范数关于初值稳定。