计算物理学刘金远课后习题答案 第6章：偏微分方程数值解法

第6章：偏微分方程数值解法

6.1对流方程

[6.1.1]考虑边值问题

ut=aux,0<x<1,Z>0

z/(x,O)—x

如果取：Av=2/7,=(y-0.5)Ax,y=1,2,3»7=8/49,tk=k\t

求出

【解】采用Crank-Nicolson方法

k~\1/kkkjt+lGJt+l.Ar+l\

2V7-'.2〃i/+"j+l+Uj-1-2«/+%+»

4/\.A

k+\,qjt+lt+1kk,*

uuu

~j-i+3%-wy+|~)-\~j+〃小

由边界条件：/(0/)=0,取幺——=0,U\=u^k=0,1,---

Ax

W(1,Z)=1,=1

*+l

-

--1100•••"w00000…

-s(l+2s)-s0…%s(1-2s)s0…小

0-s(l+2s)-s0…“20s(l-2s)s0…〃2

••••.・—••••・•

-s(l+2s)-ss(l-2s)s

01011

2I

由初始条件：w；=x.=y(j--),j=1,2,3,

1-2(Ar)2

1-「-1°

--11000'“oo0000■Mo

-13-100%1-1100

0-13-10=01-110

u2u2

-I3-1%1-11〃3

01人011

oo,t5

2〃：-=y>-w；+3〃；-Wj=-wj+w?=y»-wl+3〃；-uA一〃，+I=—

7

145!51169

III—Uj—9

91~9191

6.2抛物形方程

[6.2.1)分别用下面方法求定解问题

dud1u

a7=a7

z/(x,O)=4x(1-x)0<x<1,/>0

z/(O,r)=w(l,/)=O

(1)取心=0.2,人=1/6用显式格式计算〃；；

(2)取以=0.2,Ar=0.01用隐式格式计算两个时间步。

【解】

k+\kk、kk

(1)%士=%-弋，令"令/","=1/Ax,则

AZAx-

u：*'=Xw*.,+(1-2九)〃:+九“3,A/=XAx2

=〃：=0,w,°=4.r,(l-x,),x,=Ax/,i=1,2,3,4

^=4^(1-^)=0.8/(1-0.2/)

计算程序difu.m

(2)对于隐式格式

一入〃:二+(1+2X)w*+l-入“3=U-,A/=XAx2

计算程序idifu.m和tri.m

【6.2.2］数值求解扩散方程：

弛二蹙,0<x<i

dtdx~

边界条件：型’?二:A：初始条件：(|)(x,0)=100

建议显式格式取：Ax=0.1,r=1/2,隐式格式：Zkx=0.1,/*=1<>

解析解为：①（X"）=-sin（〃兀x）exp（-〃2兀7）,n=2^+1

[6.2.3]求下列热传导方程的数值解

U,=67Wxr,0<X<Xf,O<t<tf

〃(0")=d"),〃0)=为”)

〃(x,O)=〃o(x)

【解】

（1）采用显式FTCS方法

整理得：

〃/=口心+(1-2厂)〃：+r〃二，/=-1,r=aAt/Ax2

计算程序EFTCS.m和dcmo_EFTCS.m

(2)采用隐式FTCS方法

整理得:

2

一厂〃3+（1+2r）u'-r心=U：T,r=abtiAx,i=1,2,…，A/—1

若〃o，〃”的值满足Dirichlet型边界条件、则得到矩阵形式方程

1+2r-r

r1+2r-r

-r1+2/

计算程序：IBTCS.m和demo_IBTCS.m

对于给定牛顿边值问题（即一端给定一阶导数值，另一端给定函数值）

“。"）二标），*

得到:

心二〃;-2&叫；a“),

由上面式得：

-ru\4-(1+2rM-，*=%"，

得到：

(1+2r)u"-ru；=+ru\=14T+ru"-2rAxZ>0(/„),

整理得：

(1+2rM-2rwf=〃:”-2rAxb()(t,l),

这时矩阵形式方程为:

对于边界导数值其他定义的矩阵方程会发生相应的变化。

(3)采用两层加权平均即Grank-Nicholson格式

〃+1nI«­!■!cw+1.n+1〃c〃.〃

u1一"ia-2〃,+u,〃h]—Lii+//,)

-------=----------;---1---+-1---i--;-----

A/21_Ax'Ax"

整理得：

一2”+2(1+r)-/*：；=尸〃3+2(1-r)u'-+，心,r=a\t/Ax2,Z=-1

对于第一类边界条件：已知〃0，〃”，有关系

+2(1+/)〃；一1-rw7+,=U〃:+2(1-r)u"+ru\

_哂)2+2(1+r)〃仁一W；丁=2+2(1-尸)〃；"+ruy

得到：

M+l

2+2尸

对于给定边界导数值的矩阵方程要相应发生变化。

【6.2.4】求解偏微分方程定解问题：

cu3、〃

~=~0<x<2,?>0

100x,0<x<l

w(x,0)=•

200-100x,l<x<2

w(0,/)=0;i/(l,r)=0

(1)写出显式和隐式FTCS格式

(2)时间步长取T=0.08,空间步长取〃=0.4时，用显式FTCS格式计算公，

(保留4位小数)。

【解】(1)

显式格式：

女+1k.

U—U1,,,

L

~~-(wy+i-2wy+wy_1)=0

z/；+,='"+]+(1-2人+入〃;t,X=­r

h~

隐式格式：

A+I_k

^^^-2(媪-2叱+明)=0

Tn

一Q，：；：+(1+2X)My+l-入〃言=7T

h

2

(2)时间步长取0.08,空间步长取0・4,力=——，n=4,%=ih,i=0,1,2,3,4,5

〃+1

0.08利用显式计算，〃7］吟+〃3)

=0.5

-0.42

*=(0408080400),〃；=(0406060400),

片=(0305050300)

【6.2.5］给出热传导方程：

ATa2T

-----1/•

dl-dx2

显式和隐式差分(FTCS格式)及稳定性条件

【解工显式FTCS格式：

j二

即：7；"＜邛+入(7；\—2邛+7；、)

MAx

X=KA//(AX)2,收敛和稳定性条件是九41/2

隐式FTCS格式:

邛“一T.由二3al乜X

-K)

\tA.V

即：一乜竽+(1+2X)7；/+,-九=T：

[6.2.6]采用有限差分方法求解

11

Jcu/dt=du/dx

z/(0,r)=z/(l,/)=0,M(X,0)=sinnx

§63椭圆方程

[6.3.1]求解laplace方程

边界是：x轴，y轴和直线x+y=l；

边界条件是：u(x,y)=x+/；

取Ax=Ay=〃=1/4,采用有限差分方法求三个内点％J,〃2/，%,2的值。

【解】采用中心差分方法，得

1z、

+〃川j++WI,y+I)

其中UiJ=u(ihJh),/J=0,1,2,3,4

由边界条件

…"）二章

1I51

=记",。=〃。♦2=『〃/=心=?肛,2=5

“11

4〃u=Wo.i+5+«1.0+小=V+忆+%,2

o

“H

4%=工

,H

4〃2,1="1.1+〃31+l/2.0+“2,2=Wi.l+~3~

o

3355

得：=x必2=6;而

[6.3.2]数值求解波方程

①〃=①n,0<x<1J>0

黑晨壮。;初始条件沏

边界条件为：0<E

①(1J)=O

解析解为：①(x,/)=sin(兀x)cos(兀Z)

【6.3.3］数值求解Helmholtz方程

r)2u(x,y)f)2u(x,y),、,、,,、，，，，

—TT^+—TT^+g(x/)〃(x/)=/(x,刃，x<x<x,y<y<y

dxdy0f0f

u(x,y)=b(y)u(x,y)=b(x)

边界条件为:0x0Qy0

u(x,y)=b(y)'

fxfu(x,yf)=byf(x)

【解】

毛=Xo+iAr，M=%+iAy,其中­二(》/一/)/%,^y=(yf-y0)/My

〃/+lj-2Hp+“1/

“+A/

整理得:

%J=G(〃卬+〃ZJ)+G(/J+I+%j_i)+q，(gi/M,j一/J)

〃0J="0(乂)，〃如,广既(%)Ax2

"=2("+AJ)2(Ar2+A/)

%0=4o(x)〃i.My=byf(X,)

2(Ax2+Ay2)

数值计,算取：0<^<4,0<j^<4,

22

〃(0,y)=/,IZ(4,y)=16cos(y),f(x,y)=x+y

z/(x,0)=x2,w(x,4)=16cos(x)g(x,y)=yfx

计算程序；Heimho1tz.m和demoHeimho1tz.m

如果将g=inline。sqrt(x)'x'y');

换为g二inlineCO,,,x，/y,)；即可求解泊松方程。

§6.4非线性偏微分方程

[6.4.1]用泰勒展开求下面方程的显示差分格式

1、2

【解】w(x,f+Af)=w(x,t)+\tu,(x,/)+—(Ar)-u„(x,/)

2

uH=(u-ux)t=2uu,-uxl

22

=2W(W-WV)-(W-WV)V

=2,-4wwv+%

所以,

22

w(x,/+A/)=M(X,r)+A/(w-wx)+—(A/)-+〃.)

=(〃+△〃/+(Af))/)一(AZ+2(Afy〃)〃,+y(A/)2z/

A+I)2+(加)2(〃：)3-善(1+2AM)(嗔一〃3)+

U：+i-+心)

、72Ax

【6.4.2】二维泊松方程:

暮+.=-")

axay

在x和y方向上分别取空间步长心,Ay,a=(zVv/Ay)2,采用中心差分方法给出

雅可比迭代格式，高斯-塞德尔迭代格式和超松弛迭代格式，超松弛因子3取值

范围

【解】：⑴设忆=岫,匕),£」=/(再，为)，

G+i,/-2a,j+4-1,/a.2+i-2忆+必

+.-)

(Ax)2(缈)

aJ二」+ij+Mj+a(©j+i+")+(—)",