二阶偏微分方程的分类

范文一：二阶偏微分方程的分类

§3 二阶偏微分方程的分类

一、二阶偏微分方程的分类、标准形式与特征方程考虑二阶偏微分方程

(1) 式中aij(x)=aij(x1,x2,…,xn)为x1,x2,…,xn的已知函数.

[特征方程·特征方向·特征曲面·特征平面·特征锥面] 代数方程

称为二阶方程(1)的特征方程；这里a1,a2,…,an是某

些参数，且有特征方程，即

.如果点x?=(x1?,x2?,…,xn?)满足

则过x?的平面的法线方向

l:(a1,a2,…,an)称为二阶方程的特征方向；如果一个(n)维曲面，其每点的法线方向都是特征方向，则称此曲面为特征曲面；过一点的(n)维平面，如其法线方向为特征方向，则称这个平面为特征平面，在一点由特征平面的包络组成的锥面称为特征锥面. [n个自变量方程的分类与标准形式] 在点P(x1?,x2?,…,xn?)，根据二次型

(ai为参量)

的特征根的符号，可将方程分为四类：

(i) 特征根同号，都不为零，称方程在点P为椭圆型.

(ii) 特征根都不为零，有n个具有同一种符号，余下一个符号相反，称方程在点P为双曲型.

(iii) 特征根都不为零，有点P为超双曲型.

个具有同一种

符号(n>m>1)，其余m个具有另一种符号，称方程在

(iv) 特征根至少有一个是零，称方程在点P为抛物型.

若在区域D内每一点方程为椭圆型，双曲型或抛物型，则分别称方程在区域D内是椭圆型、双曲型或抛物型.

在点P作自变量的线性变换可将方程化为标准形式：

椭圆型：

双曲型：

超双曲型：

抛物型：

式中Φ为不包含二阶导数的项.

[两个自变量方程的分类与标准形式] 方程的一般形式为

(2)

a11,a12,a22为x,y的二次连续可微函数，不同时为零. 方程

a11dy

a12dxdy+a22dx=0

称为方程(2)的特征方程.特征方程的积分曲线称为二阶方程(2)的特征曲线.

在某点P(x0,y0)的邻域D内，根据Δ=a122－a11a12的符号将方程分类：

当Δ>0时，方程为双曲型；当Δ=0时，方程为抛物型；当Δ

在点P的邻域D内作变量替换，可将方程化为标准形式：

（i）（i）双曲型：因Δ>0，存在

两族实特征曲线

，

式

或

和

，

，作变换

（ii）（ii）抛物型: 因Δ=0，只存

在一族实的特征曲线数

，使

，取二次连续可微函

，

，作变换

方程化为标准形式

（iii）（iii）椭圆型：因Δ

在实特征曲线，设

为换

，

的积分，

不同时为零，作变量替

范文二：二阶抛物型偏微分方程

表3( 数学学院(系、所) 硕士研究生课程简介

课程名称:二阶抛物型偏微分方程课程代码:011.538 英文名称:Second Order Parabolic Differential Equations 课程类型:?讲授课程 ?实践,实验、实习，课程 ?研讨课程 ?专题讲座 ?其它考核方式: 考试教学方式: 讲授适用专业:基础数学，应用数学，计算数学适用层次: 硕士 ? 博士 ? 开课学期: 秋总学时讲授学时: 4848 学分:3 / / 先修课程要求:广义函数与SOBOLEV空间;偏微分方程概论

课程组教师姓名职称专业年龄学术专长

汤燕斌教授应用数学 42 偏微分方程

杨茵教授应用数学 48 偏微分方程

段志文副教授应用数学 41 偏微分方程

韩淑霞讲师应用数学 35 偏微分方程

课程教学目标:

本课程的目标是掌握二阶抛物型偏微分方程的基本理论;掌握二阶抛物算子

pq的基本解、抛物型偏微分方程的极值原理和-衰减估计;掌握抛物型偏微LL

分方程定解问题的古典解和弱解的适定性和解析半群方法，为线性算子半群和微分算子理论、非线性发展方程、反应扩散方程、偏微分方程数值计算以及调和分析等领域的研究打下坚实的基础。

教学大纲:(章节目录)

第一章二阶线性抛物型方程

1.1 预备知识以及抛物型方程定解问题介绍

1.2 弱解的存在性

1.3 解的正则性

1.4 极大值原理

第二章半线性抛物型方程

2.1 半线性抛物方程的抽象理论

2.2 半线性抛物方程的初值问题

2.3 半线性抛物方程的初边值问题

第三章拟线性抛物型方程

3.1 分数幂算子与分数幂空间

3.2 解的整体存在性

3.3 解的正则性

第四章全空间中的抛物型方程

4.1 极大值原理

4.2 柯西问题与半群方法

4.3 抛物型HOLDER空间以及抛物型插值不等式

4.4 热方程在HOLDER空间中的可解性

4.5 SCHAUDER估计

4.6 解的存在性

4.7 内估计

4.8 解的正则性

教材:

1. Lawrence C. Evans，《Partial Differential Equations》，Providence，RI: American

Mathematical Society，1998。

2. 王明新编著，《算子半群与发展方程》，北京:科学出版社，2006。

主要参考书:

1. Gary M.Lieberman，《Second order parabolic differential equations》， Singapore:

World Scientific Publishing,1996.

2. 陈亚浙编著，《二阶抛物型偏微分方程》，北京:北京大学出版社，2003。 3. 伍卓群、尹景学、王春朋著，《椭圆与抛物型方程引论》，科学出版社，2003。

该课程所属基层教学组织(教研室、系)专家小组意见:

该课程适合硕士、博士研究生培养的需要，不与其他课程重复，有稳定授课教师队伍。

专家组长

专家 2007年12月25日

范文三：偏微分方程数值解

一类二阶常微分方程的差分格式

1. 引言

考虑二阶常微分边值问题： ?最简形：

d 2u

Lu =-2+q (x ) u (x ) =f a

u (a ) =α, u (b ) =β, (1.2)

其中q (x ) ≥0, f (x ) 为[a , b ]上已知的连续函数，α,β为已知常数。 ?变系数：

d 2u

Lu =-2+p (x ) u '(x ) +q (x ) u (x ) =f a

u (a ) =α, u (b ) =β, (2.2)

其中p (x ), q (x ) ≥0, f (x ) 为[a , b ]上已知的连续函数，α,β为已知常数。 ?守恒形：

Lu =-

d du

[a (x ) ]+q (x ) u (x ) =f a

u (a ) =α, u (b ) =β, (2.4)

其中a (x ), q (x ) ≥0, f (x ) 为[a , b ]上已知的连续函数，α,β为已知常数。

2. 差分格式

将区间[a , b ]分成N 等分（N 为正整数），则对步长h =?最简形：

对充分光滑的解u ，由Taylor 展开式在x i 处展开得：

a -b

，将I =[a , b ]进行网格剖分，则网格节点为x i =a +ih ,0≤i ≤N . N

u (x i +1) =u (x i ) +u '(x i )(x i +1-x i ) +

111

u ''(x i )(x i +1-x i ) 2+u '''(x i )(x i +1-x i ) 3+u (4)(ξ1)(x i +1-x i ) 42! 3! 4!

1114) u ''(x i ) h 2+u '''(x i ) h 3+u ((ξ1) h (x 4i <>

111

u (x i -1) =u (x i ) +u '(x i )(x i -1-x i ) +u ''(x i )(x i -1-x i ) 2+u '''(x i )(x i -1-x i ) 3+u (4)(ξ2)(x i -1-x i ) 4

2! 3! 4! 111(4)234

=u (x i ) -u '(x i ) h +u ''(x i ) h -u '''(x i ) h +u (ξ2) h (x i -1<>

2! 3! 4!

=u (x i ) +u '(x i ) h +由 ① + ②：

h 4(4)

u (x i +1) +u (x i -1) =2u (x i ) +u ''(x i ) h +u (ξ)(x i -1≤ξ≤x i +1, ξ=max{ξ1, ξ2})

③

u (x i +1

h 4(4)

) -2u x ((1=u ) ''x i h (+u ξx i (≤) ξ(≤x i +ξi +) u x i --11=

, ξm ξ1a x , }) { ④2

u (x i +1) -2u (x i ) +u (x i -1) h 2(4)

=u ''(x i ) +u (ξ)(x i -1≤ξ≤x i +1, ξ=max{ξ1, ξ2}) ⑤

h 212

用二阶差商代替x i 点处二阶导数得到：

u (x i +1) -2u (x i ) +u (x i -1) h 2(4)

u ''(x i ) =-u (ξ)(x i -1≤ξ≤x i +1, ξ=max{ξ1, ξ2}) ⑥ 2

h 12

带入原方程：

Lu i =-

u (x i +1) -2u (x i ) +u (x i -1)

+q (x i ) u (x i ) =f (x i ) +R i (x u ) ， ⑦ 2

h 2(4)

u (ξi ) . 其中R n (x i ) =12

则当h →0时，R n (x i ) 是h 的二阶无穷小量，舍去R n (x i ) 即得到逼近方程的差分方程：

L h u i =-

u i +1-2u i +u i -1

+q i u i =f i (i =0,1, , N -1) ⑧ 2

R n (x i ) 即为截断误差。

考虑到边值条件u 0=α, u N =β，组成关于u i 的差分格式：

u i +1-2u i +u i -1

+q i u i =f i (i =0,1, , N -1)

⑨ h 2

u 0=α, u N =βL h u i =-

其解u i 是u (x ) 于x =x i 的近似，称⑨为逼近。 ?变系数：

对充分光滑的解u ，由Taylor 展开式在x i 处展开得：

①

u (x i +1) =u (x i ) +u '(x i )(x i +1-x i ) +

②

u ''(ξ1)(x i +1-x i ) 2=u (x i ) +u '(x i ) h +u ''(ξ1) h 2(x i <>

u ''(ξ2)(x i -1-x i ) 2=u (x i ) -u '(x i ) h +u ''(ξ2) h 2(x i -1<>

u (x i -1) =u (x i ) +u '(x i )(x i -1-x i ) +

由 ① - ②：

u (x i +1) -u (x i -1) =2u '(x i ) +u ''(ξ) h 2(x i -1≤ξ≤x i +1, ξ=max{ξ1, ξ2})

③

u (x i +1) -u (x i -1) h

=u '(x i ) +u ''(ξ)(x i -1≤ξ≤x i +1, ξ=max{ξ1, ξ2}) ④

2h 4

用一阶差商代替x i 点处一阶导数得到：

u '(x i ) =

u (x i +1) -u (x i -1) h

-u ''(ξ)(x i -1≤ξ≤x i +1, ξ=max{ξ1, ξ2}) ⑤

2h 4

带入原方程得到⑥：

Lu i =-

u (x i +1) -2u (x i ) +u (x i -1) u (x i +1) -u (x i -1)

+p (x ) +q (x i ) u (x i ) =f (x i ) +R i (x u ) i

h 22h

h 2(4)h

u (ξi ) +u ''(ξi ) . 其中R n (x i ) =124

则当h →0时，R n (x i ) 是h 的二阶无穷小量，舍去R n (x i ) 即得到逼近方程的差分方程：

L h u i =-

u i +1-2u i +u i -1u i +1-u i -1

+p +q i u i =f i (i =0,1, , N -1) ⑦ i

h 22h

R n (x i ) 即为截断误差。

考虑到边值条件u 0=α, u N =β，组成关于u i 的差分格式：

u i +1-2u i +u i -1u i +1-u i -1

+p +q i u i =f i (i =0,1, , N -1) i

⑧ h 22h

u 0=α, u N =βL h u i =-

其解u i 是u (x ) 于x =x i 的近似，称⑨为逼近。 ?守恒形：

对充分光滑的解u ，由Taylor 展开式在x i 处展开得：

u (x i +1) =u (x =u (x

i +

'1) +u (x

i +

1)(x i +1-x i ) +22

u ''(x 1)(x i +1-x i ) 2+o ((x i +1-x 1) 2)

i +i +2! 22

h h

') +u (x ) +u ''(x 1) +o (h 2) 11

i +i +i +28222

①

u (x i ) =u (x =u (x

i +

'1) +u (x

i +

1)(x i -x 2

i +

1) ++2

u ''(x 1)(x i -x 1) 2+o ((x i -x 1) 2)

i +i +i +2! 222

h h

') -u (x ) +u ''(x 1) +o (h 2) 11

i +i +i +28222

②

由 ① - ②：

u (x i +1) -u (x i ) =hu '(x 1) +o (h 2)

i +2

③

u (x i +1) -u (x i )

=u '(x 1) +o (h ) ④

i +h 2

用一阶差商代替x

i +12

点处一阶导数得到：

u '(x

i +

1) =2

u (x i +1) -u (x i )

-o (h ) ⑤

i -12

同理得到x

点处一阶导数⑥：

u '(x

i -

1) =2

u (x i ) -u (x i -1)

-o (h )

用二阶差商代替x i 点处二阶导数得到⑦：

d du [a (x ) ]i =dx dx

[a (x )

u (x i +1) -u (x i ) u (x i ) -u (x i -1) du du

]x i +1/2-[a (x ) ]x i -1/2a (x i +1/2) -a (x i -1/2) =

h h

带入原方程得到：

Lu i =-

a (x i +1/2)

u (x i +1) -u (x i ) u (x i ) -u (x i -1) -a (x i -1/2)

+q (x i ) u (x i ) =f (x i ) +R i (x u ) ，

h 2(4)

u (ξi ) . 其中R n (x i ) =12

则当h →0时，R n (x i ) 是h 的二阶无穷小量，舍去R n (x i ) 即得到逼近方程的差分方程：

L h u i =-

a i +1/2

u i +1-u i u -u

-a i -1/2i i -1

+q i u i =f i (i =0,1, , N -1) ⑧

R n (x i ) 即为截断误差。

考虑到边值条件u 0=α, u N =β，组成关于u i 的差分格式：

u i +1-u i u -u

-a i -1/2i i -1

L h u i =-+q i u i =f i (i =0,1, , N -1) ⑨

u 0=α, u N =β

a i +1/2

其解u i 是u (x ) 于x =x i 的近似，称⑨为逼近。

3. 格式求解

?最简形：

-[u i +1-2u i +u i -1]+h 2q i u i =h 2f i (i =0,1, , N -1) -u i +1+(2+h q i ) u i -u i -1=h f i (i =0,1, , N -1)

考虑u 0=α, u N =β与⑩可写为

⑩

?2+h 2q 1?

?-1?????

?变系数：

-12+h 2q 2

-12+h 2q N -2

-1

??u 1??h 2f 1+α??????2

u h f 2??2????? ?=?? ?????2

u -1??N -2??h f N -2?

?h 2f +β???u 2+h 2q N -1?N -1????N -1?

-2(u i +1-2u i +u i -1) +hp i (u i +1-u i -1) +2h 2q i u i =2h 2f i (i =0,1, , N -1) -2u i +1+4u i -2u i -1+hp i u i +1-hp i u i -1+2h 2q i u i =2h 2f i (i =0,1, , N -1) (hp i -2) u i +1+(2h 2q i +4) u i -(2+hp i ) u i -1=2h 2f i (i =0,1, , N -1)

考虑u 0=α, u N =β与⑩可写为

⑩

?hp 1-24+2h 2q 1

hp 2-2?

? ???

?守恒形：

4+2h 2q 2

-hp 2-2 4+2h 2q N -1

??u 1??2h 2f 1+α??????2

u 2h f 22??? ?=??? ??? ?????2

u -hp N -1-2?2h f +β??N -1??N -1???

-a i +1/2(u i +1-u i ) +a i -1/2(u i -u i -1) +h 2q i u i =h 2f i

-a i +1/2u i +1+(a i +1/2+a i -1/2+h q i ) u i -a i -1/2u i -1=h f i

考虑u 0=α, u N =β与⑩可写为

⑩

?-a 3?2???????

a 3+a 1+h 2q 1

-a 5

a 5+a 3+h 2q 1

-a 3

1N -

3N -

+h 2q N -1

??u ??h 2f +α?

??1???2

??u 2??h f 2??? ?=?? ?????2

u ?N -1???h f N -1+β??-a 3?N -?2?

4. 数值例子

?最简形：

计算如下两点边值问题：

-u ''(x ) +u (x ) =e x (sinx -2cos x ),0<πu (0)="0," u="" (π)="">

其精确解为u (x ) =e x sin x 。

将区间[0,π]做N 等分，记h =π/N , x i =ih ,0≤i ≤N . 绝对误差最大值记为E ∞(h ) =

max u (x ) -u

0≤i ≤N

表1列出了4个结点处的精确解和取不同步长时所得的数值解。表2给出不同步长时在这4个结点处所得数值解的绝对误差。表3给出不同步长时在这4个结点处所得数值解的最大误差图1给出了取不同步长时所得的数值解曲线。

图2给出了精确解曲线和取h=π/10时所得数值解曲线。图3给出了取不同步长时所得数值解的误差曲线。

表1 部分结点处的精确解和取不同步长时所得的数值解

表2 取不同步长时部分结点处数值解的误差绝对值

表3 取不同步长时部分结点处数值解的最大误差

图1

图

图3

?变系数：

计算如下两点边值问题：

-u ''(x ) +u '(x ) +u (x ) =e x (2sinx -cos x ),0<πu (0)="0," u="" (π)="">

其精确解为u (x ) =e sin x 。

将区间[0,π]做N 等分，记h =π/N , x i =ih ,0≤i ≤N . 绝对误差最大值记为E ∞(h ) =

max u (x ) -u

0≤i ≤N

图2给出了精确解曲线和取h=π/10时所得数值解曲线。

图3给出了取不同步长时所得数值解的误差曲线。

表1 部分结点处的精确解和取不同步长时所得的数值解

表2 取不同步长时部分结点处数值解的误差绝对值

表3 取不同步长时部分结点处数值解的最大误差

图1

图2

图3

?守恒形：

计算如下两点边值问题：

d du

[2]+u (x ) =e x (sinx -4cos x ),0<>

dx dx

u (0)=0, u (π) =0. -

其精确解为u (x ) =e sin x 。

将区间[0,π]做N 等分，记h =π/N , x i =ih ,0≤i ≤N . 绝对误差最大值记为E ∞(h ) =

max u (x ) -u

0≤i ≤N

图2给出了精确解曲线和取h=π/10时所得数值解曲线。

图3给出了取不同步长时所得数值解的误差曲线。

表1 部分结点处的精确解和取不同步长时所得的数值解

表2 取不同步长时部分结点处数值解的误差绝对值

表3 取不同步长时部分结点处数值解的最大误差

图1

图

图3

5. 结论

由最简形和守恒形的图1，图2对比可发现，当步长越接近于0，其得到的解距离精确解的误差就越小。观察图3可发现误差的变化随着步长的变小而趋于小且稳定。因此，该差分方法对于二阶常微分方程的边值问题，解决方案十分稳定可靠。而变系数形应该在分析的时候产生了机制上的误差，导致不稳定，由目前看来并不是十分可靠。

参考文献：

[1] 孙志忠. 偏微分方程数值解法[M]. 北京：科学出版社， 2005 [2] 李荣华. 偏微分方程数值解法[M]. 北京：高等教育出版社， 2004

[3] 易昆南. 基于数学建模的数学实验[M]. 长沙：中国铁道出版社， 2014

附录：

Matlab 源代码： M 文件：

%计算结点函数，n 为选择步长 function f=buchang(n)

a=0;b=pi;alpha=0;beta=0;%输入边值条件 h=(b-a)/n;%步长 x=a+h:h:b-h;%结点 f=zeros(1,n-1);%初始化f for i=1:n-1 format long

f(i)=exp(x(i))*(sin(x(i))-2*cos(x(i)))*h^2; end

f(1)=alpha+f(1); f(n-1)=beta+f(n-1); %三角矩阵

P=(h^2+2)*ones(n-1,1); Q=-1*ones(n-2,1);

A=diag(P,0)+diag(Q,1)+diag(Q,-1); format long

u=A\f';%步长数值解曲线 for i=1:n-1

if h*i==pi/5||h*i==2*pi/5||h*i==3*pi/5||h*i==4*pi/5 %选择结点pi/5,2pi/5,3pi/5,4pi/5输出 disp(u(i)) end end

%计算误差和最大误差函数 function f=wucha(n) a=0;b=pi;alpha=0;beta=0; h=(b-a)/n; x=a+h:h:b-h;

f=zeros(1,n-1);

U=zeros(1,n-1); %初始化精确解 format long for i=1:n-1

U(i)=exp(x(i))*sin(x(i)); %精确解求解 end

for i=1:n-1 format long

f(i)=exp(x(i))*(sin(x(i))-2*cos(x(i)))*h^2; end

f(1)=alpha+f(1); f(n-1)=beta+f(n-1); P=(h^2+2)*ones(n-1,1); Q=-1*ones(n-2,1);

A=diag(P,0)+diag(Q,1)+diag(Q,-1); format long u=A\f'; for i=1:n-1 format long

if h*i==pi/5||h*i==2*pi/5||h*i==3*pi/5||h*i==4*pi/5 disp(abs(U(i)-u(i)')) %计算结点误差并输出 end end

disp(max(abs(U-u'))); %输出最大误差

%取不同步长时数值解曲线函数作图函数 function f=tu(n)

a=0;b=pi;alpha=0;beta=0; h=(b-a)/n; x=a+h:h:b-h; f=zeros(1,n-1); for i=1:n-1 format long

f(i)=exp(x(i))*(sin(x(i))-2*cos(x(i)))*h^2;

f(1)=alpha+f(1); f(n-1)=beta+f(n-1); P=(h^2+2)*ones(n-1,1); Q=-1*ones(n-2,1);

A=diag(P,0)+diag(Q,1)+diag(Q,-1); format long u=A\f'; switch n

%不同步长时选择不同样式作图 case 10

plot(x,u,'g' );hold on ; case 20

plot(x,u,'--y' );hold on ; case 40

plot(x,u,'-.r' );hold on ; case 80

plot(x,u,':m');hold on ; case 160

plot(x,u,'xb' );hold on ; end

%精确解曲线和取h=π/10时所得数值解曲线作图函数 function f=tu3(n) a=0;b=pi;alpha=0;beta=0; h1=(b-a)/n;

h2=(b-a)/160; %使精确解U(x)图形标准的步长选择 x1=a+h1:h1:b-h1;

x2=a+h2:h2:b-h2; %U(x)的x 轴 f=zeros(1,n-1); U=zeros(1,159); format long for i=1:159

U(i)=exp(x2(i))*sin(x2(i));

U=[0,U,0]; %补上边值条件 for i=1:n-1 format long

f(i)=exp(x1(i))*(sin(x1(i))-2*cos(x1(i)))*h1^2; end

f(1)=alpha+f(1); f(n-1)=beta+f(n-1); P=(h1^2+2)*ones(n-1,1); Q=-1*ones(n-2,1);

A=diag(P,0)+diag(Q,1)+diag(Q,-1); format long u=A\f'; u=u';

u=[0,u,0]; %补上边值条件 x1=[0,x1,pi]; %补上边值条件 x2=[0,x2,pi]; %补上边值条件

plot(x1,u,'-or' ,x2,U, '-.b' );hold on ; %精确解曲线和取h=π/10数值解曲线作图

%不同步长时数值解与精确解误差作图函数 function f=tu2(n) a=0;b=pi;alpha=0;beta=0; h=(b-a)/n; x=a+h:h:b-h; f=zeros(1,n-1); U=zeros(1,n-1); format long for i=1:n-1

U(i)=exp(x(i))*sin(x(i)); end U=[0,U,0]; for i=1:n-1 format long

f(i)=exp(x(i))*(sin(x(i))-2*cos(x(i)))*h^2;

f(1)=alpha+f(1); f(n-1)=beta+f(n-1); P=(h^2+2)*ones(n-1,1); Q=-1*ones(n-2,1);

A=diag(P,0)+diag(Q,1)+diag(Q,-1); format long u=A\f'; u=u'; u=[0,u,0];

sss=abs(U-u); %误差曲线 x=[0,x,pi]; switch n

%不同步长时的误差曲线作图 case 10

plot(x,sss,'g' );hold on ; case 20

plot(x,sss,'--y' );hold on ; case 40

plot(x,sss,'-.r' );hold on ; case 80

plot(x,sss,':m');hold on ; case 160

plot(x,sss,'xb' );hold on ; end

主函数： clear buchang(10) buchang(20) buchang(40) buchang(80) buchang(160) clear

wucha(10) wucha(20) wucha(40) wucha(80) wucha(160) clear tu(10) tu(20) tu(40) tu(80) tu(160) clear tu3(10) tu3(20) tu3(40) tu3(80) tu3(160) clear tu2(10) tu2(20) tu2(40) tu2(80) tu2(160)

范文四：偏微分方程数值解

第一章概述

大家采用下面的方法求解Terzaghi 一维固结方程。

1．1 偏微分方程工具箱的功能

偏微分方程工具箱(PDE Toolbox) 提供了研究和求解空间二维偏微分方程问题的一个强大而又灵活实用的环境。PDE Toolbox的功能包括：

(1) 设置PDE (偏微分方程) 定解问题，即设置二维定解区域、边界条件以及方程的形式和系数；

(2) 用有限元法 (FEM) 求解PDE 数值解；

(3) 解的可视化。

无论是高级研究人员还是初学者，在使用PDE Too1box时都会感到非常方便。只要PDE 定解问题的提法正确，那么，启动MATLAB 后，在MATLAB 工作空间的命令行中键人pdetool ，系统立即产生偏微分方程工具箱(PDE Toolbox)的图形用户界面(Graphical User Interface，简记为GUI) ，即PDE 解的图形环境，这时就可以在它上面画出定解区域、设置方程和边界条件、作网格剖分、求解、作图等工作，详见1．4节中的例子。我们将在第二章详细介绍GUI 的使用，在第二章给出大量典型例子和应用实例。除了用GUI 求解PDE 外，也可以用M 文件的编程计算更为复杂的问题，详见第三章和第四章

的内容。

1．2 PDE Toolbox求解的问题及其背景

1．2．1 方程类型

PDE Toolbox 求解的基本方程有椭圆型方程、抛物型方程、双曲型方程、特征值方程、椭圆型方程组以及非线性椭圆型方程。

椭圆型方程： -??(c ?u ) +au =f , in Ω, ，

椭圆型方程：-??(c ?u ) +au =f , in Ω,

其中Ω是平面有界区域，c , a , f 以及未知数u 是定义在Ω上的实（或复）函数。抛物型方程：d ?u -??(c ?u ) +au =f , in Ω. ?t

?2u 双曲型方程：?2-??(c ?u ) +au =f , in Ω. ?t

特征值方程：-?(c ?u ) +au =λdu , in Ω,

其中d 是定义在Ω上的复函数，λ是待求特征值。在抛物型方程和双曲型方程中，系数c , a , f 和d 可以依赖于时间t 。

可以求解非线性椭圆型方程：

-??(c (u ) ?u )+a (u ) =f (u ), in Ω,

其中c , a , f 可以是未知函数u 的函数。

还可以求解如下PDE 方程组；

??-??(c 11(u ) ?u 1)-??(c 12(u ) ?u 2)+a 11u 1+a 12u 2=f 1, ? ??-??(c 21(u ) ?u 1)-??(c 22(u ) ?u 2)+a 21u 1+a 22u 2=f 1

利用命令行可以求解高阶方程组。对于椭圆型方程，可以用自适应网格算法，还能与非线性解结合起

来使用。

另外，对于Poission 方程还有一个矩形网格的快速求解器。

1．2．2 边界条件

（1）Dirichlet 条件： h u =r

( 2 ) Neumann 条件： n ?(c ?u ) +qu =g

其中n 是Ω的边界?Ω上的单位外法向量，g , q , h 和r 是定义在?Ω上的函数。对于特征值问题仅限于齐次条件：g =0, 和r =0。对于非线性情形．系数g , q , h 和r 可以依赖于u ；对于抛物型方程和双曲型方程，系数可以依赖于时间t 。

对于方程组情形，边界条件为

( 1 ) Dirichlet 条件： h 11u 1+h 12u =2 r h 21u 1+h 22u =2 r

( 2 ) Neumann 条件： n ?(c 11?u 1) +n ?(c 12?u 2) +q 11u 1+q 12u 2=g 1

n ?(c 21?u 1) +n ?(c 22?u 2) +q 21u 1+q 22u 2=g 2

( 3 ) 混合边界条件为: h 11u 1+h 12u 2=r 1

n ?(c 11?u 1) +n ?(c 12?u 2) +q 11u 1+q 12u 2=g 1+h 11μ

n ?(c 21?u 1) +n ?(c 22?u 2) +q 21u 1+q 22u 2=g 2+h 12μ

其中μ的计算要使得Dirichlet 条件满足。在有限元法中，Dirichlet 条件也称为本质边界条件，Neumann 条件称为自然边界条件。

1.3 如何使用FDE Toolbox

1.3.1 定解问题的设置

员简单的办法是在PDE Tool上直接使用图形用户界面(GUl)。设置定解问题包括三个步骤：

(1)Draw模式：使用CSG(几何结构实体模型) 对话框画几何区域，包括矩形、圆、椭圆和多边形，也可以将它们组合使用。

(2)Boundary模式：在各个边界段上给出边界条件，

(3)PDE模式：确定方程的类型、系数c ，a ，f 和d c。也能够在不同子区域上设置不同的系数(反映材料的性质) 。

1.3.2 解PDE 问题

用GUI 解PDE 问题主要经过下面两个过程（模式）

(1)Mesh模式；生成网格．自动控制网格参数。

(2)Solve模式：对于椭圆型方程还能求非线性和自适应解。对于抛物型和双曲型力程．设置初始边值条件后能求出给定t 时刻的解。对于特征值问题，能求出给定区间内的特征值；求解后可以加密网格再求解。

1.3.3 使用Toolbox 求解非标准的问题

对于非标准的问题。可以用PDE Too1box的函数。或者用FEM(有限元法) 求解更为复杂的问题。

1.3.4 计算结果的可视化

从GUI 能够使用Plot 模式实现可视化。可以使用Color, Height 和Vector 等作图。对于抛物型和双曲型方程，还可以生成解的动画。这些操作通过命令行都很容易实现。

1.3.5 应用领域

在应用界面提供了丁如下应用领域

．结构力学——平面应力问题

．结构力学——平面应变问题

．静电场问题

．静磁场问题

．交流电磁场问题

．直流导体介质问题

．热传导问题

．9‘散问题

这些界面都有对话框，它包括PDE 的系数、边界条件、解的性质等。

1. 4 解偏微分方程的一个例子

解Poisson 方程-?u =f ，边界条件为齐次Dirichlet 类型。

第一步：启动MATLABl, 键入pdetool ，按回车键确定便可启动GUI ，然后在Options 菜单下选择Grid 命令，打开栅格, 栅格的使用，能使用户容易确定所绘图形的大小，如图1—

1--1

第二步：分步完成平面几何造型：R1-C1-E1+R2+C2。用菜单或快捷工具，分别画矩形R1、矩形R2、椭圆E1、圆C1、圆C2。画圆时，首先选中椭圆工具，按鼠标右键并拖动即可、或者在按ctrI 的同时，拖动鼠标也可绘制圆。然后在Set formula栏，进行编辑并用算术运算将将图形对象名称连接起来，删除默认的表达式键入R1-C1-E1+R2+C2，按等号健得到所需图形。若需要，还可进行储存．

形成M 文件。

选择Boundary 菜单中Boundary Mode 命令，进入边界模式。单击Boundary 菜单中Remove A11 Subdomain Borders 选项，去除子域边界。如果想将几何信息和边界信息进行存储，应选择Boundary 菜单中的ExPort Decomposed Geometry．Boundary Cond’s …命令，将它们分别储存于g,b 变量中, 通过MA TLAB 形成M 文件。

第三步：选取边界．单击Boundary 菜单中Specify Bounddy Conditions…选项，打开Boundary conditlons对话框，输入边界条件，如图1—4。本例取缺省条件。即将全部边界设为齐次Dirichlet 条件，边界颜色显示为红色。

第四步：选择PDE 菜单中PDE Mode命令，进入PDE 模式。单击PDE 菜单中PDE Specification…选项，打开PDE 对话框，设置方程类型。本例取缺省设置，类型为椭圆型，参数c , a , f 分别为1，0，10。

第五步：选择Mesh 菜单中Initialize Mesh命令，进行网格剖分。

第六步：选择Mesh 菜单中Refine Mesh命令，对网格加密。

第七步：选择Solve 菜单中So1ve PDE命令，解偏微分方程并显示图形解。

第八步：单击Plot 菜单中Parameters ?选项，打开Plot selection对话框，选中Color, Height (3—D Plot)和Show mesh三项。然后单击Plot 按钮，显示三维图形解。

第九步：如果要画等值线图和矢量场图，单击Plot 菜单中Parameters ?选项，打开Plot Selection对话框．选中Contour 和Arrows 两项。然后单击P1ot 按钮，可显示解的等值线图和矢量场图。

第二章 PDE 图形用户界面

2.1 PDE Toolbox菜单

File 菜单（如图1-1）

图1-1

New 新建一个几何结构实体模型（Constructive Solid Geomery, 简记为CSG ），默认文件名为

“Untitled ”。

Open ? 从硬盘装载M 文件

Save 将在GUI 内完成的成果储存到一个M 文件中。

Save As? 将在GUI 内完成的成果储存到另外一个M 文件中。

Print ? 将PDE 工具箱完成的图形送到打印机内进行硬拷贝。

Exit 退出PDE 工具图形用户界面。

2 Edit 菜单（如图1-2）

图1-2

Undo 在绘制多边形时退回到上一步操作。

Cut 将已选实体剪切到剪贴板上。

Copy 将已选实体拷贝到剪贴板上。

Paste ? 将剪贴板上的实体粘贴到当前几何结构实体模型中。

Clear 删除已选的实体。

Select All 选择当前几何结构实体造型CSG 中的所有实体及其边界和字域。

3 Options 菜单（如图1-3）

图1-3

Grid 绘图时打开或关闭栅格。

Grid Spacing? 调整栅格的大小。

Snap 打开或关闭捕捉栅格功能。

Axes Limits? 设置绘图轴的坐标范围。

Axes Equal 打开或关闭绘图方轴。

Turn off Toolbar Help 关闭工具栏按钮的帮助信息。

Zoom 打开或关闭图形缩放功能。

Application 选择应用的模式。

Refresh 重新显示PDE 工具箱中的图形实体。

4 Draw 菜单（如图1-4）

图1-4

Draw Mode 进入绘图模式。

Rectangle/square 以角点方式画矩形/方行（Ctrl+鼠标）。

Rectangle/square（centered ）

以中心方式画矩形/方行（Ctrl+鼠标）。

Ellipse/circle 以矩阵角点方式画椭圆/圆（Ctrl+鼠标）。

Ellipse/circle（centered ）以中心方式画椭圆/圆（Ctrl+鼠标）。

Polygon 画多边形，单击鼠标右键可封闭多边形。

Rotate ? 旋转已选的图形。

Export Geometry Description,Set Formula，Labels ?

将几何描述矩阵gd 、公式设置字符sf 和标识空间矩阵ns 输出到主工作空间

去。

单击Draw 菜单中Rotate …选项，可打开Rotate 比对活框，通过输入旋转的角度，可使选择的物体按输入的角度逆时针旋转。旋转中心的选择如果缺省，则为图形的质心，也可以输入旋转中心坐标。 5 Boundary 菜单（如图1-5）

Boundary Mode 进入边界模式。图1-5

Specify Boundary Conditions? 对于已选的边界输入条件，如果没有选择边界，则边界条件适用于所有的

边界。

Show Edge Labels 显示边界区域标识开关，其数据是分解几何矩阵的列数。

Show Subdomain Labels 显示子区域标识开关，其数据是分解几何矩阵中的子域数值。

Remove Subdomain Border 当图形进行布尔运算时，删除已选取的子域边界。

Remove All Subdomain Borders 当图形进行布尔运算时，删除所有的子域边界。

Export Decomposed Geometry，Boundary Cond’s ?

将分解几何矩阵g 、边界条件矩阵b 输出到主工作空间。

选择Boundary 菜单中Specify Boundary Conditions．命令可定义边界条件。在打开的Boundary condition对话框，可对已选的边界输入边界条件。共有如下三种不同的条件类型：

NeMmann 条件这里边界条件是由方程系数q 和g 确定的，在方程组的情况下(换成方程组模式) ，q 是2ⅹ2矩阵，g 是2x1矢量。

Dirichlet 条件 u 定义在边界上，边界条件方程是价h*u=r，这里h 是可以选样的权因子(通常为1) 。在方程组情况下，h 是2x2矩阵，r 是2x l矢量，

混合边界条件(仅适合于方程组情形) 它是Dirichlet 和Neumann 的混合边界条件，q 是2x 2矩阵，g 是2x1矢量，h 是1x 2矢量，r 是一个标量。

6 PDE 菜单（如图1-6）

图1-6

PDE Mode 进入偏微分方程模式。

Show Subdomain Labels 显示子区域标识开关。

PDE Specification? 调整PDE 参数和类型。

Export PDE Coefficients? 将当前PDE 参数c,a,f,d 输出到主工作空间，其参数变量为字符类型。 7 Mesh 菜单（如图1-7）

图1-7

Mesh Mode 输入网格模式。

Initialize Mesh 建立和显示初始化三角形网格。

Refine Mesh 加密当前三角型网格。

Jiggle Mesh 优化网格。

Undo Mesh Change 退回上一次网格操作。

Display Triangle Quality 用0~1之间数字化的颜色显示三角形网格的质量，大于0.6的网格可接受的。 Show Node Labels 显示网格节点标识开关，节点标识数据是点矩阵p 的列。

Show Triangle Labels 显示三角形网格标识开关，三角形网格标识数据是三角形矩阵t 的列。 Parameters ? 修改网格生成参数。

Export Mesh 输出节点矩阵p 、边界矩阵e 和三角形矩阵t 到主工作空间。

8 Solve 菜单（如图1-8）

图1-8

Solve PDE 对当前的几何结构实体CSG 、三角形网格和图形解偏微分方程。

Parameters ? 调整PDE 的参数。

Export Solution? 输出PDE 的解矢量u 。如果可行，将计算的特征值1输出到主工作空间。

1.1.9 Plot 菜单（如图1-9）

图1-9

Plot Solution 显示图形解。

Parameters ? 打开绘图方式对话框。

Export Movie? 如果动画被录制了，则动画矩阵M 将输出到主工作空间。

10 Window 菜单

从Window 菜单项，可选择当前打开的所有的MA TLAB 图形窗口，被选择的窗台移至前台。

11 Help 菜单

Help ? 显示帮助信息

About ? 显示版本信息

1．2 PDE 工具栏

主菜单下是工具栏，工具栏中喊有许多工具图标按钮，可提供快速、便捷的操作方式。从左到右5个按钮为绘图模式按钮，紧接着的6个为边界、网格、解方程和图形显示控制功能按钮，最右边的为图形缩放功能键。（如图1-10）

图

1-10

以角点方式画矩形/方行（Ctrl+鼠标）。

以中心方式画矩形/方行（Ctrl+鼠标）。

以矩形角点长轴方式画椭圆/圆（Ctrl+鼠标）。

以中心方式画椭圆/圆（Ctrl+鼠标）。

画多边形，按右键可封闭多边形。

进入边界模式。

打开PDE Specification（偏微分方程类型）对话框。

初始化三角形网格。

加密三角形网格。

解偏微分方程。

打开Plot Selection对话框，确定后给出解的三维图形。

为显示缩放切换按钮。

第三章典型方程及其应用实例

求解PDE 问题主要有两种方法，一种是使用图形用户界面，另一种是采用命令行编程。前者直观简便，而后者更为灵活。

2．1 求解椭圆方程的例子

例：单位圆上的Poisson 方程边值问题：

22??-?u =1, Ω={(x , y )|x +y <1}??u |?ω="0">

这一问题的精确解为：

u (x , y )1-x (=2-y 2)

若使用图形用户界面（Graphical User Interface, 简记为GUI ），则首先在MA TLAB 的工作窗口中键入pdetool ，按回车键确定，于是出现PDE Toolbox窗口。如果需要坐标网格，单击Options 菜单下的Grid 选项即可。下面分步进行操作。

（i ）画区域圆单击工具，大致在（0，0）位置单击鼠标右键同时拖拉鼠标到适当位置松开，绘制圆。为了保证所绘制的圆是标准的单位圆，在所绘圆上双击，打开Object Dialog对话框，精确地输入圆心坐标X-center 为0、Y-cebter 为0及半径Radius 为1，然后单击OK 按钮，这样单位远已画好。

（ii ）设置边界条件单击工具，图形边界变红，逐段双击边界，打开Boundary Condition对话框，输入边界条件。对于同一类型的边界，可以按Shift 键，将多个边界同时选择，统一设置边界条件。本题选择Dirichlet 条件，输入h 为1，r 为0，然后单击OK 按钮。也可以单击Boundary 菜单中Specify Boundary Conditions ?选项，打开Boundary Condition对话框输入边界条件，如图2-1。

（iii ）设置方程单击PDE 菜单中PDE Specification?选项，打开PDE Specification对话框，选项方程类型。本题单击Elliptic ，输入c 为1，a 为0，f 为1，然后单击OK 按钮，如图2-2。

图

2-1

图2-2

（iv ）网格剖分单击工具，或者单击Mesh 菜单中Initialize Mesh选项，可进行初始网格剖分，这时在PDE Toolbox窗口下方的状态栏内显示初始问网格的节点数和三角形单元数。本题节点数为144个，三角形单元数为254个。如果需要网格加密，再单击，或者单击Mesh 菜单中Refine Mesh 选项，这时节点数变为541个，三角形单元数为1016个，如此还可继续加密。

（v ）解方程单击工具，或者单击Solve 菜单中Solve 菜单中Solve PDE 选项，可显示方程色彩解。如果单击Plot 菜单中Parameters ?选项，出现Plot Selection对话框，如图2-3，从中可以选择Color ，Contour ，Arrows,Deformed mesh，Height （3-D polt），还可以设置等值线的数目等。本例中选择Color ，Contour ，Height （3-D polt）和Show mesh四项，然后单击Plot 按钮，方程的图形解如图2-4所示。除了作定解问题解u 的图形外，也可以作|grad u|，|cgrad u|等图形。

图

2-3

图2-4

（vi ）与精确解作比较单击Plot 菜单中Parameters ?选项，打开Plot Selection对话框，在Height （3-D plot ）行Property 下拉框中选user entry，且在该行的User entry输入框中键入u-（1-x.^2-y.^2）/4,单击Ｐlot 按钮就可以看到解的绝对误差图形，如图２-５．可见在边界处误差为０。

图2-5

（vii ）输出网格节点的编号、单元编号以及节点坐标单击Mesh 菜单中Show Node Labels选项，再单击工具

或

，即可显示节点编号。若要输出节点坐标，只需单击Mesh 菜单中Export Mesh?选项，

这时打开的Export 对话框中默认值为p ，e ，t ，这里p ，e ，t 分别表示points （点）、edges （边）、triangles （三角形）数据的变量，单击OK 按钮。然后在MA TLAB 命令窗口键入p ，按回车键确定，即可显示出节点按编号排列的坐标（二维数组）；键入e ，按回车键，则显示边界线段数据矩阵（7维数组）；输入t ，按回车键，则显示三角形单元数据矩阵（4维数组）。

（viii ）输出解的数值单击Solve 菜单中Export Solution?选项，在打开的Export 对话框中输入u ，单击OK 按钮确定。再在MATLAB 命令窗口中输入u ，按回车确定，即显示按节点编号排列的解的数值。

我们也可以用MATLAB 程序求解PDE 问题，同时显示解的图形： [p,e,t]=initmesh(‘circleg’,’hmax’,1); Error=[];err=1; While err>0.001,

[p,e,t]=refinemesh(‘circleg’,p,e,t); U=assempde(‘circleb1’,p,e,t,,1,0,1); Exact=(1-p(1,:).^2-p(2,:).^2)’/4; Err=norm(u-exact,’inf’); Error=[error err]; End

Pdemesh(p,e,t) Pdesurf(p,t,u) Pdesurf(p,t,u-exact)

通过命令行键入help+命令函数，如help pdemesh ，按回车键，可以调入有关命令函数的定义、参数格式等帮助信息。

2．2 求解抛物型方程的例子

例：考虑一个带有矩形孔的金属板上的热传导问题。板的左边保持在100c ，板的右边热量从板向环境空气定常流动，其他边及内孔边界保持绝缘。初始

t t 0时板的温度为0，于是概括为如下定解问题：

??d ???????????

-?u =0, ?t

u =100, ?u

=-1, ?n ?u =0, ?n

u |t =t 0=0.

域Ω的外边界顶点坐标为（-0.5，-0.8），（0.5，-0.8），（0.5，0.8），（-0.5，0.8）。内边界顶点坐标为（-0.005，-0.4），（0.05，-0.4），（0.05，0.4），（-0.05，0.4）。

使用GUI 求解这一问题。在PDE Toolbox窗口的工具栏中选择Generic Scalar模式。（i ）区域设置单击

工具，在窗口拖拉出一个矩形，双击矩形区域，在Object Dialog对话框中输

入Left 为-0.5，Bottom 为-0.8，Width 为1，Height 为1.6，单击OK 按钮，显示矩形区域R1。用同样方法作内孔R2，只要设置Left=-0.05，Bottom=-0.4，Width=0.1，Height=0.8即可。然后在Set formula栏中键入R1-R2。

（ii ）设置边界条件单击

，使边界变红色，然后分别双击每段边界，打开Boundary Cvondition

对话框，设置边界条件。在左边界条件。在左边界上，选择Dirichlet 条件，输入h 为1，r 为100；右边界上，选择Neumann 条件，输入g 为-1，q 为0；其他边界上，选择Neumann 条件，输入g 为0，q 为0。

（iii ）设置方程类型单击

，打开PDE Specification对话框，设置方程类型为Parabolic （抛物型），

d=1，c=1，a=0，f=0，单击OK 按钮。

（iv ）网格剖分单击

，或者加密网格，单击

。

（v ）初值和误差的设置单击Solve 菜单中Parameters ?选项，打开Solve Parameters 对话框，输入Time 为0:5，u （t0）为0，Relative tolerance为0.01，Absolute tolerance为0.001，然后单击OK 按钮。

（vi ）数值解的输出单击Solve 菜单中Export Solution?选项，在Export 对话框中输入u ，单击OK 按钮。再在MA TLAB 命令窗口中键入u ，按回车键，这时显示出按节点编号的数值解。

（vii ）解的图形单击Plot 菜单中Parameters ?选项，打开Plot Selection对话框，选Color ，Contour ，Arrows ，单击Plot 按钮，窗口显示出Time=5时解的彩色图形，如图2-6。

图2-6

MATLAB 程序：

[p,e,t]=initmesh(‘crackg ’);

u=parabolic(0,0:0.5:5,‘crackb ’,p,e,t,1,0,0,1);

pdeplot(p,e,t,‘xydata ’,u(:,11),‘mesh ’, ‘off ’, ‘colormap ’, ‘hot ’) 2．3 求解双曲型方程的例子

例：考虑如下二维波动方程的定解问题：

??2u

?2-?u =0, (x , y ) ∈Ω={(x , y ) |-1<1,><1},>

u |x =±1=0, ?

?u ?|y =±1=0, ?

?n ?

π?

t =0, u =(x )), ?2

?π

sin(x ). ?u ?

=3sin(πx ) e 2

??t ?

用GUI 求解。类似前面的例子，首先作正方形区域：设置断点坐标为（-1，1），（-1，-1），（1，-1），（1，1）。在Object Dialog对话框中输入Left 为-1，Bottom 为-1，Width 为2，Height2，单击OK 按钮。设置边界条件：左、右边界用Dirichlet 条件，输入h 为1，r 为0；上、下边界用Neumann 条件，输入q 为0，g 为0，作网格剖分。设置方程类型为Hyperbolic （双曲型），键入c=1，a=0，f=0，d=1。

单击Solve 菜单中Parameters ?选项，打开Solve Parameters对话框，在Time 栏中键入linspace （0，5，31），设置u 的初始值u （t0）为atan(cos(pi/2*x)),u’的初始值u ’(t0)为3*sin(pi*x).*exp(sin(pi/2*y)),Relative tolerance 为0.01，Absolute tolerance为0.001。

最后输出图形解：单击Plot 菜单中Parameters ?选项，打开Plot Selection对话框，选Color ，Contour ，Arrows ，Height （3-D Plot）和Show mesh五项，单击Plot 按钮，三维彩色图形的解如图2-7所示。

图2-7

MATLAB 的动画程序： [p,e,t]=initmesh(‘squareg ’); x=p(1,:)’;y=p(2,:)’; u0=atan(cos(pi/2*x));

ut0=3*sin(pi*x).*exp(sin(pi/2*y)); n=3l;

tlist=linspace(0,5,n);

uu=hyperbolic(u0,ut0,tlist,‘squareb3’,p,e,t,1,0,0,1); delta=-1:0.1:1;

[uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta); gp=[tn;a2;a3]; umax=max(max(uu)); umin=min(min(uu)); newplot M=moviein(n); For I=1:n,

pdeplot(p,e,t,‘xydata ’,uu(:,I),‘zdata ’,uu(:,I), ‘gridparam ’,gp, ‘colorbar ’, ‘off ’, ‘zstyle ’, ‘continuous ’);

axis([-1 1 –1 1 umin umax]);caxis([umin umax]); M(:,I)=getframe; End

Movie(M,10);

2．4 求解特征值问题的例子

‘mesh ’, ‘off ’, ‘xygrid ’,

‘on ’,

例：混合边界条件的特征值问题：

?-?u =λu , (x , y ) ∈Ω={(x , y ) |-1<1},?u |x="-1=0,">

|y =±1=0, ?

?n ?

??u 3??

-u ?|x =1=0, ???n 4? ?

首先作区域：方形的4个顶点为（-1，1），（-1，-1），（1，-1），（1，1）。在Object Dialog对话框中输入Left 为-1，Bottom 为-1，Width 为2，Height2。

单击

，逐段双击边界设置边界条件：左边界选Dirichlet 条件，输入h 为1，r 为0；右边界选Neumann

条件，输入q 为-3/4，g 为0；上、下边界选Neumann 条件，输入q 为0，g 为0。

在PDE Specification对话框中设置方程类型为Eigenmodes （特征值模式），PDE 的系数设置为c=1，a=0，d=1。

由于右边界上Neumann 条件，特征值可以出现负值，因此在求小于10的特征值时应在Solve Parameters对话框中键入[-inf 10]。

作网格剖分，然后单击

，得到第一个特征值对应的特征函数图形，如图2-8。

图2-8

再在Plot Selection对话框的Eigenvalue 项中选取不同的特征值，比如第二个特征值“2.06，并选Color ，Contour ，Height （3-D Plot）和Show mesh四项，如图2-9，可作出第二个特征值对应的特征函数的三维图形，如图2-10。

图

2-9

图2-10

MATLAB 程序：

[p,e,t]=initmesh(‘squareg ’);

[p,e,t]=refinemesh（‘squareg ’,p,e,t ）; [v,l]=pdeeig(‘squareb2’,p,e,t,1,0,1,[-inf 10]); pdesurf(p,t,v(:,4))

2．5 偏微分方程在一维空间的应用

已知偏微分方程初始条件时间t 和一维空间变量x ，MATLAB 自身存在对偏微分方程的解函数pdepe 。单一方程

假设我们要求解热传导方程

u =u

u （t , 0）=0，u （t , 1）=1

u （0, x ）=

2x 1+x

（1.1）

MATLAB 指定抛物型的偏微分方程形式如下：

c (x , t , u , u x ) u t =

边界条件为：

-m

x b (x , t , u , u x ) ）+s (x , t , u , u x )

p (x l , t , u ) +q (x l , t ) ?b (x l , t , u , u x ) =0p (x r , t , u ) +q (x r , t ) ?b (x r , t , u , u x ) =0

其中

x 为边界左端点，x 为边界右端点，初始条件为 u (0，x ) =f (x ) 。（我们注意到同一函数b 在方程

和边界条件中出现。）通常情况下，每个函数都会被指定到不同的MATLAB 文档中。也就是说，与方程有关的函数c,b 和s 都会被指定保存到一个MATLAB 文件中，与边界条件有关的函数p 和q 保存到令一个文档中（此外，我们需要注意到函数b 是相同的只需保存一次即可），最后将初始函数f (x ) 保存在第三个文档中。执行命令函数pdepe 将把m 个文档结合起来进行运算，返回方程的解。在我们的例子中，有：

c (x , t , u , u x ) =1

b (x , t , u , u x ) =u x s (x , t , u , u x ) =0

我们将它们保存到MATLAB 文档中，记作eqn1.m (m=0以后在加详述) 。 Function[c,b,s]=eqn1(x,t,u,DuDx)

%EQN1:MATLAB function M-file that specifies %a PDE in time and one space dimension. C=1; B=DuDx; S=0;

对于边界条件，我们有：

P(0,t,u) = u ； q(0,t) = 0 P(1,t,u) = u – 1 ； q(1,t) = 0 将它们保存到MATLAB 文档中，记作bc1.m. Function [p1,q1,pr,qr]=bc1(x1,u1,xr,ur,t)

%BC1:MATLAB function M-file that specifies boundary conditions %for a PDE in time and one space dimension. P1=u1;

Q1=0; Pr=ur-1; Qr=0;

对初始条件，有：

f （x ）=

2x 1+x

将它保存到MA TLAB 文档中，记作initial1.m. Function value=initial1(x)

%INITIAL1:MATLAB function M-file that specifies the initial condition %for a PDE in time and one space dimension. Value=2*x/(1+x^2);

最终我们准备好了用pdepe 函数处理偏微分方程。在如下的MA TLAB 文档中，我们选择坐标值x 和t ，求解偏微分方程并绘制解的表面图（如图2-11） %PDE1：MA TLAB script M-file that solves and plots %solutions to the PDE stored in eqn1.m M=0;

%NOTE:m=0 specifies no symmetry in the problem.Taking %m=1 specifies cylindrical symmetry,while m=2 specifies %spherical symmetry. %

%Define the solution mesh x=linspace(0,1,20); t=linspace(0,2,10); %Solve the PDE

U=pdepe(m,@epn1,@initial1,@bc1,x,t); %Plot solution Surf(x,t,u);

Title(‘Surface plot of solution’); Xlabel(‘Distance x’); Ylabel(‘Time t’);

图2-11

通常，我们会发现绘制解的剖面图是十分有用的。此时，t 是固定的，u 是x 的函数。解u （t,x ）作为t 和x 的向量指标被保存为指标矩阵。例如，u(1,5)返回在点（t(1),x(5)）出u 的值。我们可以使用命令plot(x,u(1,:))绘制u 的最初值的图像（t=0）（见图2-12）。

其实，我们有一种比较快的方法去绘制剖面图，如下的MATLAB 序列： Fig=plot(x,u(1,:),’erase ’, ’xor ’) For k=2:length(t)

Set(fig,’xdata ’,x, ’ydata ’,u(k,:)) Pause(.5) end

试过之后，你会发现它们是如何快速地接近热传导方程的平衡位置的。（平衡位置就是及时停止并要发生改变的位置）。

图2-12 3 PDE Toolboxz中的命令简介

本节将介绍PDE Toolbox中的命令函数。利用这些命令函数，可通过命令行而不是用PDE Tool图形用户界面来解偏微分方程。

3．1 PDE Toolbox中的函数及其分类

表3-1 PDE 数值计算函数

表3-2 用户界面算法

表3-3 几何算法

表3-4

表3-5 通用算法

表3-6 用户定义算法

3．2 PDE数值计算函数简介表3-7

3．3用户界面算法函数简介 pdecirc

画圆，修改几何图形矩阵。格式：pdecirc （xc,yc,radius,label ）

说明画一个以（xc,yc ）为圆心、radius 为半径的圆。Label 为标识符，可任选，也可缺省。 pdeellip

画椭圆，修改几何图形矩阵。

格式：pdeellip （xc,yc ）为中心，以radiusx 为x 半轴、radiuy 为y 半轴的椭圆。Angle 表示半轴与坐标轴之间的夹角，label 为标识符，可任选，也可缺省。 pdepoly

画多边形，修改几何图形矩阵。格式：pdepoly （x,y,label ）

说明以向量x,y 所对应的元素作为顶点画多边形。Label 为标识符，可任选，也可缺省。 Pderect

画矩形，修改几何图形矩阵

格式：pderect(xmin,xmax,ymin,ymax,label)

说明以（xmin,ymin ）,(xmax,ymin),(xmax,ymax),(xmin,ymax)为顶点画矩形。Label 为标识符，可任选，也可缺省。 Pdetool

进入偏微分方程工具用户界面（PDE Tool GUI）格式：pdetool

说明进入PDE Tool GUI不需要其他参数，只需在主窗口中键入pdetool 即可。

使用PDE Tool GUI可以画二维区域，定义边界条件，选择偏微分方程类型，创建和加密三角形网格，计算和显示解等等。用PDE Tool GUI可以使解偏微分方程非常简单和直观。 4

有限元法和有限差分法

MATLAB 的PDE Toolbox 采用有限元法(FEM ) 求解椭圆型(elliptic)、双曲线型(parabolic)、抛物线型(hyperbolic)、特征值型(eigenvalue) 和非线性椭圆型(nonlinear elliptic) 等偏微分方程。应用PDE Toolbox 求解上述偏微分方程有两种方式。其一是使用图形用户界面(graphical user in terface, 简称GU I)。在图形用户界面上设置了绘图模式(Draw mode) , 用以绘制求解问题的定义域; 边界模式(Boundary mode) , 用来定义边界条件; 偏微分方程模式(PDE mode) , 指定方程类型、系数; 网格模式(Mesh mode) , 自动剖分区域, 产生网格; 求解模式(Solve mode) , 按有限元算法计算网格节点上的未知函数; 图形显示模式(P lo t mode) , 将计算结果以二维、三维图形的方式显示出来。另一种方法是应用命令行函数(comm and line functions) , 即利用PDE Toolbox 所提供的函数编程求解。如利用in it mesh () 函数进行剖分区域、assempde () 函数求解偏微分方程、pdesurf () 函数图形显示计算结果等。相比之下, 前一种方法更为简单, 研究人员只需在图形用户界面的各种模式下根据实际问题定义相应的参数, 而无需编程, 即可得到计算结果。4．1 椭圆型方程

PDE Tool求解的基本椭圆型方程为

-??(c ?u ) +au =f , in Ω, （4-1）

其中Ω是平面有界区域,c,a,f 以及未知函数u 可以是定义在Ω复函数。边界条件可以是以下的情形之一：

[15]

?Dirichlet 条件：hu =r , ?Ω. （4-2）

?一般（广义）Neumann 条件：n ?(c ?u ) +qu =g , ?Ω, （4-3）

当q=0时成为Neumann 条件。

?混合条件：在?Ω上部分是Dirichlet 条件，部分是Neumann 条件。

其中n 是边界?Ω的单位法向量，g,q,h 和r 是定义在?Ω上的函数。在有限元法中，Dirichlet 条件也称为本质边界条件，Neumann 条件称为自然边界条件。

有限元解实际上是微分方程弱形式的解（弱解或广义解）在有限维空间的投影导出定解问题对应的弱形式

我们讨论方程（4-1）在一般Neumann 条件（4-3）下的解。取任意试验函数

v ∈V ，乘（4-1）式的两边，并在Ω上积分：（?（?c ?u )) v +auv ]dx =?fvdx . ?[-

利用Green 公式及条件（4-3），可得

[c ?u ) ??u +auv ]dx -?(-qu +g ) vds =?fvdx . ?（

?Ω

于是，问题（4-1），（4-3）的弱形式（虚功方程）为：

求u ∈V , 满足

?[(c ?u ) ??v +auv -fv ]dx -?(-qu +g ) vds =0, ?v ∈V .

（4-4）

（4-4）的解也称为问题（4-1），（4-3）的广义解（或弱解）。在一定意义下，问题（4-1），（4-3）与它的弱形式（4-4）是等价的。

如果问题是自共轭的，并且满足所谓椭圆型条件，那么问题（4-1）（4-3）也可以按最小势能原理导出泛函极小问题，这时它与虚功方程是等价的。

区域剖分

由于三角形剖分在几何上有很大的灵活性，对边界?Ω的逼近较好，因此在PDE Toolbox中区域Ω一般作三角形网格剖分。一般而言，在作三角形剖分和节点编号时要注意以下几点：

? 单元顶点不能是相邻单元边上的内点。 ? 尽量避免出现大的钝角、大的边长。

? 在u （x,y ）的梯度变化比较剧烈的地方，网格要加密。

? 单元编号可以任意，但节点编号的好坏，直接影响总刚度矩阵的带宽，应尽量使所有两个相邻节

点编号之差的绝对值中的最大者愈小愈好。

插值多项式的选取

在PDE Toolbox中采用的是单元上的线性函数：u(x,y)=ax+by+c. （4-5）

设节点

p i 上u 的值为u i , 即u (x i , y i ) =u i , i =1, 2, ???, N p , N p

为节点数。任取单元e ，三个顶点为

p i , p j , p m , 记e =p i p j p m 。

它们的顺序是逆时针的。为了使插值函数（4-5）在这三个顶点上分别取值

u i , u j , u m ,

那么a,b,c 应该满足

ax i +by i +c =u i , ax j +by j +c =u j , ax m +by m +c =u m ,

解得a,b,c ，再代入（4-5）式，可得查值函数为

u (x , y ) =N i (x , y ) u i +N j (x , y ) u j +N m (x , y ) u m ,

N i =

(a i x +b i y +c i ), 2?e

y i

, y m

其中

y i 1x i x i

a i =, b i =-, c i

y m 1x m x m

可由

为

p i p j p m

的面积。

N j ，N m 及a j , b j , c j , a m , b m , c m N j

及

a j , b j , c j

的表达式按i →j →m →i 轮换得到。

{δ}e =[u i , u j , u m ]T , [N ]=[N i , N j , N m ]T , 则在e 上有u =[u ]{δ}e =[B ]{δ}e ,

设

a i ??u ?u ?1? ?u = ?x , ?y ??可以表示为?u =2? b

??e ?i

u 的梯度向量

2X3常数矩阵。

单元刚度矩阵、单元荷载向量的形成设

a j b j

a m ?

?{δ}e =[B ]{δ}e b m ??，其中[B]使

{Φi }(i =1, 2, ???, N p ) 为V 的N p 维子空间的基函数，按所剖分的单元将（4-4）式改为

[c ?u ) ???+au ?-f ?]dxdy =∑?(g ?-qu ?) ds . ∑??（

n =1e n

n =1r n

（4-6）

这里NE 是单元数，

r n =?e n ??Ω。

总刚度矩阵和总荷载向量的组装

将单元刚度矩阵和单元荷载向量表达式代入（4-6），为了便于叠加，先

{δ}e {F }e {K }e 进行扩充，扩充为N p 维向量和N p ?N p 阶方阵，然后进行叠加。在PDE Toolbox 中

这一组装过程通过命令行eassempde 来实现。

约束处理，求解方程组

对于Neumann 条件，由于是自然边界条件，?Ω上不需要满足任何约束条件，因此可立即得到线性代数方程组：KU=F 一旦形成K 和F ，在MA TLAB 环境下立即可解出节点近似解的向量u 。如果是Dirichlet 边界条件，还需要对边界节点进行约束处理，然后再解方程组。 4．2 抛物型方程

下面说明抛物型方程如何简化成椭圆型方程来求解。这是通过PDE Toolbox函数parabolic 完成的。

-??(c ?u )+au =f , in Ω,

考虑?t （4-7）

初值为

u (x , 0) =u 0(x ), x ∈Ω. （4-8）

边界条件类似椭圆边值问题的提法，这里仅讨论Neumann 条件：

n ?(c ?u ) +qu =g , ?Ω. （4-9）

对于热传导方程（4-7）可以写成

ρC

-??(k ?u ) +h (u -u ∞) =f , ?t （4-10）

表示热量向环境扩散，其中是热源。

ρ是密度，C 是比热，k 是导热系数，h 是薄层传热系数，u ∞是环境温度，f

如果系数与时间无关，方程就是标准的椭圆型方程： -??(c ?u ) +au =f .

对Ω作三角形网格剖分，对于任意给定t ≥0, PDE 的解按有限元法的基底可以展开成

u (x , t ) =∑u i (t ) ?i (x ).

（4-11）

将展开式带入方程（4-7），两边乘以试验函数可得

?i ，

然后在Ω上积分，利用Green 公式和边界条件（4-9），

∑?d ?j ?i dx

i ΩΩ

du i (t )

+∑(???j ?(c ??i ) +a ?j ?i dx +?q ?j ?i ds ) v i (t ) dx dt i Ω?Ω

（4-12）

=?f ?j dx +?g ?j ds , ?j .

?Ω

上式可写成大型的线性稀疏的常微分方程组：这就是所谓的线性半离散化方法。

再解ODE （4-13）的初值，初值为

+KU =F . dt （4-13）

U i (0) =u 0(x i ).

得每一个节点

x i 、任一时刻t 的ODE 解。这里K 和F 是远边界条件下椭圆型方程

-??(c ?u ) +au =f , in Ω

的刚度矩阵和荷载向量。M 是质量矩阵。

当边界条件是与时间有关的Dirichlet 条件hu=r时，F 项则包括h 和r 的时间倒数，它可以用有限差分求解。常微分方程组是病态的，这时需作显式时间积分。

由于稳定性要求时间间隔很小，而隐式解由于每一时间段都要求解椭圆型方程，从而求解非常缓慢。常微分方程组的数值积分，可以由MATLAB 中suite 函数完成。对于这类问题它是有效的。 4．3 双曲型方程

类似于抛物型方程的有限元法，可以解双曲型问题：

?2u

d 2-??(c ?u ) +au =f , in Ω. ?t u (x , 0) =u 0(x ),

(x , 0) =v 0(x ), x ∈Ω, ?t 边界条件同上。

初值为

?2u

-c ?u =02?t 特别地，对于波动方程，波速为c 。

对区域Ω作三角形网格剖分，与抛物型方程处理方法一样，可以得到二阶常微分方程组：

d 2V M 2+KV =F , dt V i (0) =u 0(x i ), d V i (0) =V 0(x i ), ?i . dt 初值为

K 和M 是刚度矩阵和质量矩阵。

PDE Toolbox中提供的求解双曲型方程的命令函数是hyperbolic 。

4．4 特征值方程

在PDE Tooltox中求解的基本特征值问题是

-??(c ?u )+au =λdu ，

其中λ是未知复数。在固体力学中方程描述薄膜震动的问题，在量子力学中λ描述势阱a(x)中有界定态的能级。数值解包括方程的离散和代数特征值问题的求解。首先考虑离散化。按有限元基底将u 展开，两边同乘基函数，再在Ω上作积分，可以得到广义特征值方程：

KU=λMU ，

其中对应于右端项的质量矩阵的元素为

M i , j =?d (x ) ?j (x ) ?i (x ) dx .

利用命令函数assema （单元积分贡献的装配）生成。在通常情况下，当函数d （x ）为正时，质量矩阵M 为正定对称矩阵。同样，当c （x ）是正的而且在Dirichlet 边界条件下，刚度矩阵K 也是正定的。

广义特征值问题：

KU=λMU.

利用Arnoldi 算法进行移位和求逆矩阵, 直到所有的特征值都落在用户事先确定的区间. 在此, 求解的详细过程略去.

PDE Toolbox中提供的求解特征值问题的命令函数是pdeeig 。

范文五：偏微分方程数值解

抛物型方程的向前 Euler 格式

1. 问题介绍

考虑一维热传导方程: (1)

, 0), (22

T t x f x

u a

u ≤<>

其中 a 是正常数, ) (x f 是给定的连续函数。按照定解条件的不同给法,可将方程(1)的定解问题分为两类:

第一类、初值问题 (也称 Cauthy 问题) :求具有所需次数偏微商的函数 ) , (t x u , 满足方

程(1) (∞<∞-x )和初始条件:="">

), () 0, (x x u ?= ∞<>

第二类、初边值问题(也称混合问题) :求具有所需次数偏微商的函数 ) , (t x u ,满足方

程(1) (l x <0)和初始条件:>

), () 0, (x x u ?=

l x <>

及边值条件 (4)

. 0) , () , 0(==t l u t u T t ≤≤0

假定 ) (x ?在相应区域光滑, 并且在 l x , 0=满足相容条件, 使上述问题有唯一充分光滑

的解。

2. 网格剖分及差分格式

2.1网格剖分

考虑边值问题(1) , (4) 的差分逼近。取空间步长 N l h /=和时间步长 M T /=τ,其中 N,M 都是正整数。用两族平行直线:

)

, , 1, 0(N j jh x x j ===

) , , 1, 0(M k k t t k ===τ

将矩形域 }0; 0{T t l x G ≤≤≤≤=分割成矩形网格,网格节点为 ) , (k j t x 。以 h G 表示网格内点集合,即位于开矩形 G 的网点集合; h G 表示所有位于闭矩形 G 的网点集合; h Γ=h G -h G 是网格界点集合。差分格式

用差商代替原热传导方程中的导数,就

可

以

得到差分格式。如下图所示,向前差分格式联系到第(k+1)层的点 ) , (1+k j t x 和第 k 层的点 ) , (), , (1k j k j t x t x -及 ) , (1k j t x +。其差分格式如下:

, 22

j k

j k j k j k

k j

f h

u u u a

u u ++-=--++τ

) (j j x f f =,

) (0

j j j x u ??==, 00==k

N k

u u ,

其中 j = 1,2,…, N-1, k = 1,2,… ,M-1。以 2/h a r τ=表示网比。则方程(5)可以改写为:

j k j k

j k j k j

f ru

u r ru

u τ++-+=-++1

) 21(

将其写成矩阵形式:

??????

? ??=+-+++1112111

k N k k k u u u u

??????

? ??=-k N k

k k

u u u u 121

??????

??++=-k N

N k

i ru

f ru f f 12

01τττ

?????

?---=r r

r r r r r r C 212121

i k

k f Cu

u +=+1

3. 截断误差

在建立差分格式时,我们用差商代替导数,其实是略去了泰勒展开中的小量项,略去的小量项是

) , (12) , (24

) 1(k ik u ik i ik

t x

x t

u R

ξητ??

??=

我们称这个略去的小量项为截断误差。

4. 格式的稳定性

通过误差估计方程

111

) 2(2--++++-=k j

k j k j k j k j

e e e e r e

当 2

10≤

5. 数值例子

应用向前欧拉格式解决如下问题:

已知:a=1, f(x)=0, u(x,0)=1-x, u(0,t)=1, u(1,t)=0, u(x,t)=1-x,区间 [0,1] 求数值解。数值解:

第一层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第二层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第三层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第四层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第五层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第六层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第七层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第八层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 第九层 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 误差: 第一层 0 0 0 0 0 0 0 0 0 第二层 0 0 0 0 0 0 0 0 0 第三层 0 0 0 0 0 0 0 0 0 第四层 0 0 0 0 0 0 0 0 0 第五层 0 0 0 0 0 0 0 0 0 第六层 0 0 0 0 0 0 0 0 0 第七层 0 0 0 0 0 0 0 0 0 第八层 0 0 0 0 0 0 0 0 0 第九层 0.001 0.0014 0.0009 0 0.0001 0.0014 0.0013 0.0008 0.0004 误差中最大数为:0.0014

6. 参考文献

1. 李荣华 . 偏微分方程数值解 . 北京:高等教育出版社, 2005 2. 孙志忠 . 偏微分方程数值解 . 北京:科学出版社, 2005

3. 刘卫国 . MATLAB 程序设计教程(第二版) . 北京:中国水利水电出版社, 2010

附:

程序流程图:

程序代码:

%a=1

%f(x)=0

%u(x,0)=1-x u(0,t)=1 u(1,t)=0 %u(x,t)=1-x

%分段数

c=10;

a=1;

%区间各端点

g=0;

b=1;

t=b;

%步长

h=(b-g)/c

d=h;

r=a*d/h^2

%构造结点矩阵

M=zeros(c-1,c);

%第一列

for i=1:c-1

M(i,1)=1-i*d;

end

%构造系数矩阵

v1=(1-2*r)*ones(1,c-1);

v2=r*ones(1,c-2);

A=diag(v1,0)+diag(v2,1)+diag(v2,-1) %构造 f 向量

f=zeros(c-1,1);

for i=1:c-1

if i==1

f(i,1)=r;

else

f(i,1)=0;

end

%循环求解

for i=1:c-1

u=A*M(:,i)+f;

for j=1:c-1

M(j,i+1)=u(j,1);

end

N=M';

%输出结果

N(2:c,:)

转载请注明出处范文大全网 » 二阶偏微分方程的分类