经典偏微分方程课后习题答案

第四章 抛物型微分方程有限差分法

1设已知初边值问题22, 01, 0<(,0)sin , 01(0,)(1,)0, 0 u u

x 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

, 时间2510

h =

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. 建立差分格式, 对于1

1,...1; 0,1,...,

1T

j k h

τ

=−=−, 用向前差商代

替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式:

112

2k k k k k

1

j

j j j u u u u u h ++−+=

. 变形j τ

−−有:

1112

(12) (k k k k

j j j j u ru r u ru r h τ

+−+=+−+=

(4.1)

用向后差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式最简隐格式:

1

1112

2k k k k k j j

j j j u u u u u h τ

++++−−+=

1

1

+−1k

j +,变形有:

1111(12) k k k j j j ru r u ru u ++−−−++−= (4.2)

(4.1)*0.5+(4.2)*0.5得CN 格式为:

111112

222k k k k k k k k j j

j j j j j j u u u u u u u u h τ+++−+−−++−+=

11

1

++−1k

j +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 j

k 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 j

k 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 j

k 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:M

u(m,1)=sin((m-1)*h*pi); end

u(1,1:N)=0;u(M,1:N)=0;

for n=1:N-1

for m=2:M-1

u(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再绘图,如图()