大连理工大学高等数值分析偏微分方程数值解(双曲方程书稿)

双曲型方程的有限差分法

线性双曲型方程定解问题： （a ）一阶线性双曲型方程

()0=??+??x

u

x a t u （b ）一阶常系数线性双曲型方程组

0=??+??x

t u

A u 其中A ，s 阶常数方程方阵，u 为未知向量函数。 （c ）二阶线性双曲型方程（波动方程）

()022=??

?

??????-??x u x a x t u

()x a 为非负函数

（d ）二维，三维空间变量的波动方程

0222222=????

????+??-??y u x u t u 022222222=???? ????+??+??-??z u y u x

u t u §1 波动方程的差分逼近 1.1 波动方程及其特征

线性双曲型偏微方程的最简单模型是一维波动方程：

（1.1） 22

222x

u a t u ??=?? 其中0>a 是常数。

（1.1）可表示为：022

222=??-??x

u a t u ，进一步有

0=??? ????+?????? ????-??

u x a t x a t

由于

x

a

t ??

±??当a dt dx ±=时为()t x u ,的全导数

（=dt du dt dx x u t u ???+??x

u

a t u ??±??=）,故由此定出两个方向 （1.3） a

dx dt 1

±=

解常微分方程（1.3）得到两族直线 （1．4） 1C t a x =?+ 和 2C t a x =?- 称其为特征。

特征在研究波动方程的各种定解问题时，起着非常重要的作用。 比如，我们可通过特征给出（1.1）的通解。（行波法、特征线法） 将（1.4）视为),(t x 与),(21C C 之间的变量替换。由复合函数的微分法则

2

12211C u

C u x C C u x C C u x u ??+??=?????+?????=?? x C C u C u C x C C u C u C x u ?????

? ????+????+?????? ????+????=??2

212121122 2221222122

12C u C C u C C u C u ??+???+???+??= 2

2

22122122C u

C C u C u ??+???+??= 同理可得

a t t a t C -=??-=??1，a t

C

=??2 ????

????-??=?????+?????=??21

2211C u C u a t C C u t C C u t u

t

C C u C u a C u t C C u C u a C t u ?????? ?????? ????-?????+?????? ?????? ????-?????=??2122112122 ?????????-??+????????-???-=21222

2

22221222

C C u C u a C u C C u a ???

?????+???-??=22221221

22

2C u C C u C u a 将22x u ??和22t

u

??代入（1.1）可得： ????????+???-??2222122122

2C u C C u C u a ????????+???+??=222212212

22C u C C u C u a 即有

02

12=???C C u

求其对2C 的积分得：

()11

C f C u

=?? 其中()1C f 是1C 的任意可微函数。 再求其对1C 的积分得：

（1.5） ()()11,dC C f t x u ?= ()()()()at x f at x f C f C f ++-=+=212211 其中()?1f 和()?2f 均为任意的二次连续可微函数。 （1.5）为（1.1）的通解，即包含两个任意函数的解。

为了确定函数()at x f -1和()at x f -2的具体形式，给定u 在x 轴的初值

（1.5） ()

()+∞<<∞-?????=??===x x t

u x u t t 10

00?? 将（1.5）式代入上式，则有 （ⅰ）()()()x x f x f 021?=+

注意()=t x u t ,()()()a at x f a at x f ?+'+--'21；()=0,x u t ()()()()x a x f x f 112?='-'，有

（ⅱ）()()()x a

x f x f 1121

?='-' 并对x 积分一次，得

()()()C d a

x f x f x

+=

-?ξξ?10121 与（ⅰ）式联立求解，得

()()()2

21211002C

d a x x f x ++=?ξξ?? ()()()2

21211001C

d a x x f x --=?ξξ?? 将其回代到通解中，即得（1.1）在（1.5）条件下的解：

（1.6） ()t x u , ()()[]at x at x ++-=0021??()ξξ?d a

at

x at x 121?+-

即为法国数学家Jean Le Rond d ’Alembert (1717-1783)提出的著名的D ’Alembert 公式。

由D ’Alembert 公式还可以导出解的稳定性，即当初始条件（1.5）仅有微小的误差时，其解也只有微小的改变。如有两组初始条件：

()()()()()()()()()()+∞<<∞-??

?====x x x u x x u x x u x x u t t 1

2120

101~0,,

0,~0,,

0,????

满足 δ??<-00~，δ??<-1

1~，则 ()()≤-t x u t x u ,,21

()()at x at x +++00~21??+()()at x at x -+-1

1~21?? ()()ξξ?ξ?d a at x at

x 11~21-+?+- 即

()()≤-t x u t x u ,,21=?++at a

221

2121δδδ()δt +1

显然，当t 有限时，解是稳定的。

此外，由D ’Alembert 公式可以看出，解在()00,t x 点，()00>t 的值仅依赖于x 轴上区间[]0000,at x at x +-内的初始值()x 0?，()x 1?，与其他点上的初始条件无关。故称区间[]0000,at x at x +-为点()00,t x 的依存域。它

是过点()00,t x 的两条斜率分别为a

1±的直线在x 轴上截得的区间。

对于初始轴0=t 上的区间[]21,x x ，过1x 点作斜率为

a

1

的直线at x x +=1；过1x 点作斜率为a

1

-

的直线at x x -=2。它们和区间[]21,x x 一起构成一个三角区域。此三角区域中任意点()t x ,的依存区间都落在

[]21,x x 内部。所以解在此三角形区域中的数值完全由区间[]21,x x 上的初

始条件确定，而与区间外的初始条件无关。这个三角形区域称为区间

[]21,x x 的决定域。在[]21,x x 上给定初始条件，就可以在其决定域中确定

初值问题的解。

1.2显格式

现在构造（1.1）的差分逼近。取空间步长h 和时间步长τ，用两族平行直线

jh x x j ==， ,2,1,0±±=j ，τn t t n ==， ,2,1,0=n

作矩形网络。于网点()n j t x ,处Taylor 展开成

()()()

()()

22

11,,,2,h O t x u h

t x u t x u t x u n j xx n j n j n j +=+--+ ()()()

()()

22

11,,,2,ττ

O t x u t x u t x u t x u n j tt n j n j n j +=+--+

代入（1.1）,并略去截断误差，则得差分格式： （1.7）

=+--+2

1

12τn j

n j n j u u u 2

1

12

2h u u u a

n

j n j n j -++-

,2,1,0±±=j ， ,2,1,0=n

这里n j u 表示u 于网点()n j t x ,处的近似值。初值条件（1.5）用下列差分

方程近似：

（1.8） ()j j x u 00?=

（1.9）

()j j

j x u u 101?τ

=-

注意：（1.7）的截断误差阶是()22h O +τ，而（1.9）的截断误差阶仅是

()τO 。为此需要提高（1.9）的精度，可用中心差商代替t u ，即

（1.10）

()j j

j x u u 1112?τ

=--

为了处理1

-j u ，在（1.7）中令0=n ，得

=+--2

1

012τ

j

j j u u u 2

1

0012

2h

u u u a

j j j -++-

进一步，

=+--1012j j j u u u ()0

10012-++-j j j u u u r

其中h

a r τ=

。并用（1.10）式的()

j j j x u u 11

12τ?-=-代入上式得 ()()=-+-j j j j x u x u 110122τ??()()()()1001022-++-j j j x x x r ???

即

（1.11） =1j

u ()()()()

()()j j j j x x r x x r 10210102

12

τ????+-+--+

这样，利用（1.8） (1.11)，可以由初始层()0=n 的已知值，算出第一层()1=n 各网格节点上的值。然后利用（1.7）或显式三层格式

（1.12） =+1n j u ()()1

211212--+--++n j

n j n j n j u u r u u r 可以逐层求出任意网点值。

以上显式三层格式也可用于求解混合问题：

（1.13）()()()()()()()()???

????====??=??t t l u t t u x x u x x u x

u a t u t βα??,,00,0,102

2

222

取J L h =

，N

T

=τ。除（1.7）～（1.9）外。再补充边值条件 （1.14） ()παn u N =0，()πβn u N J =

1．

3稳定性分析

下面我们要讨论（1.7）的稳定性。为引用Fourier 方法，我们把波动方程（1.1）化成一阶偏微分方程组，相应地把显式三层格式（1.7）化成二层格式。

一种简单的做法是引进变量t

u

v ??=

，于是（1.1）化为 v t

u =??，22

2x u a t v ??=?? 这样会使得初值()0,x u 与()0,x v 不适定（不唯一），更合理的方法

是再引进一个变量x

u

a

??=ω，将(1.1)化为 （1.15） t u v ??=，x a t v ??=??ω，x

v

a t ??=??ω

注意到：

x a x u a x a x u a t u t v ??=???? ????? ?

?????=??=??=??ω22

222； x

v

a x t u a t x u a x u a t t ??=???=???=???? ????? ??????=??22ω 若令???

?

??=ωv U ，???

?

??=00a a A ，则(1.5)可写成 (1.16)

0=??-??x

t U

A U 相应地，将（1.7）写成等价的双层格式：

(1.17) ???

????-=--=-+-+-+--++h v v a h a v v n j n j n j n j n j n j n j n j 11111212

12

121τωωωωτ 即

()'17.1 (

)

(

)

?????-+=-+=+-+-+--++111112121

2

1

21n j n j n j n j n j n j n j n j v v r r v v ωωωω 其中τ

1

--=

n j

n j n j u u v ，n

j 2

1-

ωh

u u a n j n j 11

1+-+-=。可直接验证之。

记h

a r τ

=

为网比。用Fourier 方法可以证明，差分方程(1.17)稳定的必要条件是网比 （1.19） h

a r τ

=1≤。 充分条件是网比

（1.19） h a r τ

=

1<。 Courant 等证明，1=r 时，差分解仍稳定，收敛。但是要求有更光滑的初值。习惯上也称1≤r 为Courant 条件或C-F-L （Courant-Fridrichs-Lewy ）条件。

稳定性条件（1.19）有直观的几何解释。从方程（1.12）

=+1n j u ()()1

211212--+--++n j

n j n j n j u u r u u r 可看出，n j u 依赖于前两层的值：11--n j u ，11-+n j u ，1-n j u ，2

-n j

u ，而这四个值由依赖于，2-n j u 依赖于：31--n j u ，31-+n j u ，3-n j u ，4

-n j

u 1-n j u 依赖于：21--n j u ，21-+n j u ，3-n j u ，2

-n j u 11-+n j u 依赖于：22-+n j u ，2-n j u ，21-+n j u ，3

1-+n j u 11--n j u 依赖于：22--n j u 2-n j u ，21-+n j u ，31-+n j u

以此类推，可知，n j u 最终依赖于初始层0=n 上的下列值：

0n j u -，01+-n j u ，… ，0j u ，… ，01-+n j u ，0

n j u +

因此，称x 轴上含于区间[]n j n j x x +-,的网点为差分解n j u 的依存域，

它是x

轴上被过()n j t x ,和()0,n j x -以及()n j t x ,和()0,n j x +的两条直线所切割下来的区间所覆盖的网域。而过()n j t x ,的两条特征线为：()n j t t a x x -±=-。

差分格式稳定的必要条件为：h a r τ

=

1≤或a

h 1≤τ，并且进而a

h 1

-≥-

τ

。 可见差分格式稳定的必要条件是：

差分解的依存域必须包含微分方程解的依存域，否则差分格式不稳定。

用依存域的概念容易证明：当1>r 时，差分解不收敛。 1.4 隐式

为了得到绝对稳定的差分格式，用第1-n 层、n 层、1+n 层的中心差商的加权平均去逼近xx u 得到下列差分格式：

=+--+2

112τn j

n j n j u u u ()????????+-++--++-----+-++-+++2

1

1

1112112111112

22212h u u u h u u u h u u u a n j n j n j n j n j n j n j n j n j θθθ 或 ()[]

1

2212222

2211

-++-+=n j

x n j x n j x n j

t

u u u h

a u θδδθθδδτ 其中10≤≤θ是参数。

可以证明，对于41≥θ时，差分格式绝对稳定；4

10<≤θ时，差分格式的充要条件是：θ

τ411-<=

h a r 。 当0=θ就是显格式（1.7），一个常用的隐式格式是取4

1

=θ此时，差分格式为：

=+--+2

112τ

n j

n j n j u u u ()()()[]1

111111111112

222224----+-++-++++-++-++-n j n j n j n j n j n j n j n j n j u u u u u u u u u h

a

或 ()

1

12222

2241

-+++=n j

n j n j x n j

t

u u u h

a u δδτ

高维波动方程！

§3 一阶双曲方程

双曲方程与椭圆方程和抛物方程的一个重要区别是，双曲方程具有特征和特征关系，其解对初值有局部依赖性质。初值的函数性质（如间断、弱间断等）也沿着特征传播，因而其解一般没有光滑性质。我们在构造双曲方程的差分逼近时，应充分注意这些特性。

下面对于一阶双曲方程，介绍几种常见的差分格式 3.1 迎风格式

首先考虑一阶线性常系数双曲方程 （3.1）

0=??+??x

u a t u 此方程虽简单，但是对我们构造差分格式很有启发。我们的主要的目的是构造差分格式，因此只限于考虑纯初值问题。 设0a ≠, 定义特征线:

a dx dt 1

= 或 dx a dt

= 则在每一条这样的特征线上,

0du u u dx u u a dt t x dt t x

????=+=+=???? 因此,在特征线上,u 等于常数.

对于（3.1）按照用差商代替微商的方法，自然有如下三种格式：

()12.3 01

1

=-+--+h

u u a

u u n

j n j n j

n j τ

（左偏心格式）

()22.3 011=-+-++h

u u a

u u n j

n j n j

n j τ

（右偏心格式）

()32.3

021

11

=-+--++h

u u a

u u n

j n j n j

n j τ

（中心格式）

其中()12.3和()22.3的截断误差的阶为()h O +τ，

()32.3的截断误差的阶为()2h O +τ。

记

（3.3） h

a r τ= 将()12.3～()32.3式改写为：

()12.3' ()n j n j n j u r ru u -+=-+111

()22.3' ()n j n j n j ru u r u 111-+-+= ()32.3' n j n j n j n j u r u r u u 1112

2

+-+-+=

用Fourier 方法分析稳定性可知，()32.3绝对不稳定。0≥a 时，()22.3不稳定，而()12.3当1≤a h

τ稳定，；0≤a 时， ()12.3不稳定，而()22.3当1

≤a h

τ

稳定。这两个稳定条件意味着差分方程的依存域必须包含微分方程的依存域。

同样的思想可用于构造变系数方程

()0=??+??x

u

x a t u 的差分格式。此时a 可能变号，因此相应的格式为：

（3.6） ???????<=-+-≥=-+-++-+0

,00,0111

1j n

j n j j n j n j

j n j n j j n j n j a h u u a u u a h u u a u u 当当τ

τ 其中()j j x a a =。

稳定性条件为

（3.7）

1max ≤j j

a h

τ

由（3.7），并取j j a h r ??

? ??=τ，则知()12.3'和()22.3'

右端的系数非负。

当0≥j a 时，

()()∞

-+∞

+=-+≤-+≤=n

n j

j

j j n j j n j j n j j

n u r r u r u r u U U m ax 11m ax 111

当0

()()∞

-+∞

+=-+≤-+≤=n

n j j

j j n j j n j j n j j

n u r r u r u r u U

U m ax 11m ax 111

其中n U 是以n j u 为分量的的向量。总之，∞

∞+≤n n U U 1

。这说明（3.6）

稳定，按气体力学的含义（a(x)表示气流速度），称（3.6）为迎风格式。

初边值问题：边值条件应该在迎风方向给出!

3.2 积分守恒的差分格式

迎风格式是根据特征走向构造出来的向前或向后差分格式。现在以积分守恒方程出发构造差分格式。

所谓守恒方程是指如下散度型偏微分方程 （3.13）

()0,=??+??x

u x f t u

设G 是xt 平面中任意有界域，由Green 公式

()dxdt x u x f t u ??? ????+????

,G

()udx fdt -=?Γ

其中G ?=Γ。于是可将（3.13）写成积分守恒方程 （3.14） ()udx fdt -?Γ

=0

1. Lax 格式

首先，我们从（3.14）出发构造所谓Lax 格式。取G 为()n j A ,1+，

()1,1++n j B ，()1,1+-n j C 和()n j D ,1-为顶点的开矩形。ABCDA =Γ为其

边界，则

（3.15） ()udx fdt -?Γ

=

()dx u DA

-?+()dx u BC

-?+fdt AB ?

+

fdt CD

?

右端第一积分用梯形公式，第二积分用中矩形公式即 ()dx u DA

-?≈h u u n j n j 22

1

1?+-

+-，

()dx u BC

-?≈()h u

n j

21

-?-+

第三、第四积分用如下矩形公式计算：

fdt AB

?

≈τ?+n j f 1，

fdt CD

?

≈()τ-?-n j f 1

从而有

h u u u n

j n j n j 22111????

? ??+-+-+()

011=?-+-+τn j n j f f 两端同除以τ?h 2得Lax 格式 （3.16）

()

τ

n j n

j n j

u u u 111

21+-++-0211=-+-+h f f n j n j 其中()()n j j n j t x u x f f ,,=，此格式的截断误差为()2h O +τ。

特别地，au f =时，Lax 格式为关于

0=??+??x

u

a t u 的显格式： ()

τ

n j n

j n j

u u u 111

21+-++-0211=-+-+h u u a n j n j 即 ()()()n j n j n j

u r u r u 111

112

1-++++-= 其稳定性条件为

1≤h

a τ

。 现在回过头来看绝对不稳定格式()32.3

τ

n j

n j u u -+1021

1=-+-+h

u u a

n

j n j

Lax 格式实际是用()n j n j u u 112

1+-+取代n

j u 的结果，这样一个变化就使得绝

对不稳定格式成为条件稳定，并保持截断误差为()2h O +τ。

双曲方程组

112121111f x

u b x u

b x u b t u m m =??++??+??+??

m m mm m m m f x

u b x u

b x u b t u =??++??+??+?? 2211 即

f u

u +??=??x

B t 若矩阵][ij b B =相似于对角矩阵，则称为双曲方程组，可以化成m 个一阶双曲方程组，分别求解。

数值分析常微分方程的数值解法

《计算机数学基础》数值部分第五单元辅导 14 常微分方程的数值解法 一.重点内容 1. 欧拉公式： ）心知1）a 儿+1 =儿 + hfg ，儿） m 1、 伙=0丄2,…川一 1） I 无=x Q +kh 局部截断误差是0（*）。 2. 改进欧拉公式： 预报一校正公式： 预报值 _v*+1 =儿+ hf （x k ,儿） - h - 校正值 y M = y k +-［f （x kt y k ） + /（x A+1, y M ）］ 即 儿+1 =儿+ ￡ "（忑'儿）+心+「儿+ hfg ,儿））］ 或表成平均的形式： 儿=儿+ hfg ,儿） '儿=儿+"（无+】，儿） +K ） 改进欧拉法的局部截断误差是0（2） 3. 龙格一库塔法 二阶龙格一库塔法的局部截断误差是0（爪） 三阶龙格一库塔法的局部截断误差是0（护） 四阶龙格F 塔法公式：儿计=儿+ 2（匕+ 2心+ 2? + ?） 四阶龙格一库塔法的局部截断误差是0（爪）。 二实例 y' = — y — xv f2（0 < x < 0.6） 例1用欧拉法解初值问题｛ ' ? -取步长/匸02计算过程保留 b （o ）= 1 4位小数。 解/i=0.2. f （x ）= —y —xy 2<,首先建立欧拉迭代格式 y*+i =儿+ hf g,y k ） = y k -hy k -hx k y ； =0.2 儿（4 - x k y k ）（k = 0,1,2） K 2=f(x n +^h, yk+-hK\)t gg+舟人，>'n +y/?A3）；

当k=0, xi=0.2 时，已知x（）=0,y（）=l,有 y(0?2)今i=0?2X l(4-0X 1)=0.8000 当k=\. M=0?4时，已知“=0?2」尸0?8,有 y(0?4)今2=0.2 X 0.8X(4-0.2X0.8)=0.614 4 当k=2, xs=0.6 时，已知x2=0.4,y2=0.6144,有 y(0?6)今3=0.2 X0.6144X (4-0.4 X 0.4613)=0.8000 「J, ,2 ?_ ZX 例2用欧拉预报一校正公式求解初值问题\y + v +V sinx=,取步长/?=0.2,计算 .y ⑴=1 y(0.2),y(0.4)的近似值，计算过程保留5位小数。 解步长力=0.2,此时/(x,y)=—y—fsiiu 欧拉预报一校正公式为： 预报值兀I = y k + hfg y k) - I J_ 校正值)3=儿+尹(忑，儿)+ fg,儿+1)] 有迭代格式] 预报值儿+] = y k 4-h(-y k -y； sin x k) =y k (0?8-0?2儿sin x k) < h 、—— 2 校止值y如]=儿 +尸[(一片一力sinxJ + LN+i-yl sin.v I+1)] ——?> =儿(°?9一0?1儿sin心)一0?1(儿+| +y；j sin心利) 、"M=0.別=1」)=1 时，Xj=1.2> 有 儿=yo(°?8-O?2yo sinx0) = 1 x (0.8-02x lsin 1) = 0.63171 y(1.2) ?= lx(0.9-0.1xlxsinl)-0.1(0.63171+0.631712sinl.2) = 0.71549 当 T xi=1.2, yi=0.71549 时,x2=1.4,有 y2 =儿(0.8-0?2儿sinXj) = 0.71549x(0.8-02x0.71549sinl.2) =0.47697 y(14) z y2 = 0.71549x(0.9-0.1x0.71549xsin 1.2)-0.1(0.47697+ 0.476972 sin 1.4) =0.52608 V = 8 — 3y 例3写出用四阶龙格一库塔法求解初值问题^ ‘的计算公式，取步长/匸0.2计 b(0) = 2 算y(0.4)的近似值。讣算过程保留4位小数。 解此处.心,刃=8 —3”四阶龙格一库塔法公式为 艰=儿 + % + 2? + 2勺 + ?) 1 h, y n+ y/?A3)： 本例计算公式为: 0 2 呱严儿+三(32?+2?+心

偏微分方程数值解期末试题及标准答案

偏微分方程数值解试题(０6B ） 参考答案与评分标准 信息与计算科学专业 一(10分）、设矩阵A 对称，定义)(),(),(2 1)(n R x x b x Ax x J ∈-=，)()(0x x J λλ?+=.若0)0('=?，则称称0x 是)(x J 的驻点(或稳定点).矩阵A 对称(不必正定)，求证0x 是)(x J 的驻点的充要条件是:0x 是方程组 b Ax =的解 解: 设n R x ∈0是)(x J 的驻点，对于任意的n R x ∈,令 ),(2),()()()(2 000x Ax x b Ax x J x x J λλλλ?+-+=+=, （3分) 0)0('=?,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0，则有0||||),(2000=-=--b Ax b Ax b Ax ，得到b Ax =0． (3分) 反之,若n R x ∈0满足b Ax =0,则对于任意的x ,)(),(2 1)0()1()(00x J x Ax x x J >+==+??,因此0x 是)(x J 的最小值点. (４分） 评分标准：)(λ?的展开式3分, 每问3分，推理逻辑性１分 二（10分）、 对于两点边值问题:?????==∈=+-=0 )(,0)(),()('b u a u b a x f qu dx du p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ],[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈ 建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和Galerkin 形式的变分方程。 解: 设}0)(),,(|{11=∈=a u b a H u u H E 为求解函数空间,检验函数空间．取),(1b a H v E ∈,乘方程两端,积分应用分部积分得到 (３分) )().(),(v f fvdx dx quv dx dv dx du p v u a b a b a ==+=??,),(1 b a H v E ∈? 即变分问题的Galerkin 形式. （3分)

偏微分方程数值解

偏微分方程数值解 偏微分方程地构建科学、工程学和其他领域的数学模型的主要手段。一般情况下,这些模型都需要用数值方法去求解。本书提供了标准数值技术的简明介绍。借助抛物线型、双曲线型和椭圆型方程的一些简单例子介绍了常用的有限差分方法、有限元方法、有限体方法、修正方程分析、辛积分格式、对流扩散问题、多重网络、共轭梯度法。利用极大值原理、能量法和离散傅里叶分析清晰严格地处理了稳定性问题。本书全面讨论了这些方法的性质,并附有典型的图像结果,提供了不同难度的例子和练习。 本书可作为数学、工程学及计算机科学专业本科教材,也可供工程技术人员和应用工作者参考。 偏微分方程数值解---学习总结(2) 关于SobolveSobolve空间的几个重要定理 迹定理 : ΩΩ是 RdRd 的一个有界开子集，具有李普希茨连续边界?Ω?Ω, s>12s>12, 则 a.存在唯一的连续线性映射γ0:Hs(Ω)→Hs?12(?Ω),满足γ0v=v ∣∣?Ω,?v∈Hs(Ω)∩C0(Ωˉˉˉˉ), b.存在唯一的连续映射R0：Hs?12(?Ω)→Hs(Ω),满足γ0°R0°φ=φ,?φ∈Hs?12(?Ω).(1)(2)(1)a.存在唯一的连续线性映射γ0:Hs(Ω)→Hs?12(?Ω),满足γ0v=v|?Ω,?v∈

Hs(Ω)∩C0(Ωˉ),(2)b.存在唯一的连续映射R0：Hs?12(?Ω)→Hs(Ω),满足γ0°R0°φ=φ,?φ∈Hs?12(?Ω). 迹定理把区域内部与边界联系起来. 上面定理中边界?Ω?Ω当被它的一个子集ΣΣ代替时，结论依然成立. S=1时， γ0:H1(Ω)→H12(?Ω)?L2(?Ω)||γ0v||0,?Ω≤||γ0v||2,?Ω≤C||v||1=C(||v||0+||?v||0).γ0:H1(Ω)→H12(?Ω)? L2(?Ω)||γ0v||0,?Ω≤||γ0v||2,?Ω≤C||v||1=C(||v||0+||? v||0). 注意几个范数 ||?||k||?||0||?||1||??||0=||?||k,2=||?||L2=||?||1,2=(||?||20+||??||20)12=|?|1.(3)(4)(5)(6)(3)||?||k=||?||k,2(4)||? ||0=||?||L2(5)||?||1=||?||1,2=(||?||02+||??||02)12(6)||?? ||0=|?|1. 庞加莱不等式(Poincare inequality): 假设ΩΩ是 RdRd 的一个有界联通开子集，ΣΣ是边界?Ω?Ω的一个非空的李普希茨连续子集. 则存在一个常数 CΩ>0CΩ>0满足 ∫Ωv2(x)dx≤CΩ∫Ω|?v(x)|2dx,?v∈H1Σ(Ω),其中H1Σ(Ω)={v ∈H1(Ω),γΣv=v∣∣Σ=0}.∫Ωv2(x)dx≤CΩ∫Ω|?v(x)|2dx,?v∈HΣ1(Ω),其中HΣ1(Ω)={v∈H1(Ω),γΣv=v|Σ=0}.

偏微分方程数值解试题及答案

偏微分方程数值解试题(06B) 参考答案与评分标准 信息与计算科学专业 一（10分）、设矩阵A 对称，定义)(),(),(2 1 )(n R x x b x Ax x J ∈-= ，)()(0x x J λλ?+=.若0)0('=?,则称称0x 是)(x J 的驻点（或稳定点）.矩阵A 对称（不必正定），求证0x 是)(x J 的驻点的充要条件是：0x 是方程组 b Ax =的解 解: 设n R x ∈0是)(x J 的驻点,对于任意的n R x ∈,令 ),(2 ),()()()(2 000x Ax x b Ax x J x x J λλλλ?+ -+=+=, (3分) 0)0('=?,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0,则有 0||||),(2000=-=--b Ax b Ax b Ax ,得到b Ax =0. (3分) 反之,若 n R x ∈0满足 b Ax =0,则对于任意的 x ,)(),(2 1 )0()1()(00x J x Ax x x J >+ ==+??,因此0x 是)(x J 的最小值点. (4分) 评分标准:)(λ?的展开式3分, 每问3分,推理逻辑性1分 二（10分）、 对于两点边值问题：????? ==∈=+-=0 )(,0)() ,()(' b u a u b a x f qu dx du p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ] ,[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈ 建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和 Galerkin 形式的变分方程。 解: 设}0)(),,(|{11 =∈=a u b a H u u H E 为求解函数空间,检验函数空间.取),(1 b a H v E ∈,乘方程两端,积分应用分部积分得到 (3分) )().(),(v f fvdx dx quv dx dv dx du p v u a b a b a ==+=??,),(1 b a H v E ∈? 即变分问题的Galerkin 形式. (3分) 令?-+=-=b a dx fu qu dx du p u f u u a u J ])([21),(),(21)(22,则变分问题的Ritz 形式

常微分方程数值解法

i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。

偏微分方程数值解试题06B答案

专业班级 姓名 学号 开课系室数学与计算科学学院 考试日期

偏微分方程数值解试卷 一（15分）、（1）简述用差分方法求解抛物型方程初边值问题的数值解的一般步骤.(2)写出近似一阶偏导数 n m x u |??的三种有限差分逼近及其误差阶,写出近似 n m x u |22 ??的差分逼近及其误差阶. 评分标准： （1） 7分，三个离散4分，其他步骤3分 （2） 8分，每个格式及误差2分。 二（15分）、（1）以抛物型方程的差分格式为例，解释差分格式的相容性,稳定性和收敛性概念,分析相容性,稳定性和收敛性与误差的关系，简述 Lax 等价性定理。(2) 简述差分格式稳定性分析的Fourier 级数法（或称为Neumann Von 方法，分离变量法）的一般步骤。 （1）8分，解释概念6分，等价关系2分 （2）7分，典型波2分，放大因子与条件3分，其他2分 三（20分）、对于边值问题 ?? ???=?=∈=??+???0 |) 1,0()1,0(),(,92 222G u G y x y u x u （1）建立该边值问题的五点差分格式（五点棱形格式又称正五点格式），推导截 断误差的阶。 （2）取3/1=h ，求边值问题的数值解（写出对应的方程组的矩阵形式并求解） （3）就取5/1=h 的情况写出对应方程组的系数矩阵（用分块矩阵表示）。 解：（1）7分，离过程与格式

第二页（共五页） 四（20分）、对于初边值问题??? ????≤≤==<<=≤<<

数值分析第五版全答案chap1

第一章 绪 论 1．设0x >,x 的相对误差为δ，求ln x 的误差。 解：近似值*x 的相对误差为* **** r e x x e x x δ-=== 而ln x 的误差为()1ln *ln *ln **e x x x e x =-≈ 进而有(ln *)x εδ≈ 2．设x 的相对误差为2%，求n x 的相对误差。 解：设()n f x x =，则函数的条件数为'()||() p xf x C f x = 又1'()n f x nx -= , 1 ||n p x nx C n n -?∴== 又((*))(*)r p r x n C x εε≈? 且(*)r e x 为2 ((*))0.02n r x n ε∴≈ 3．下列各数都是经过四舍五入得到的近似数，即误差限不超过最后一位的半个单位，试指 出它们是几位有效数字：*1 1.1021x =,*20.031x =, *3385.6x =, *456.430x =,*5 7 1.0.x =? 解：*1 1.1021x =是五位有效数字； *20.031x =是二位有效数字； *3385.6x =是四位有效数字； *456.430x =是五位有效数字； *57 1.0.x =?是二位有效数字。 4．利用公式(2.3)求下列各近似值的误差限：(1) *** 124x x x ++,(2) ***123x x x ,(3) **24 /x x . 其中****1234,,,x x x x 均为第3题所给的数。 解：

*4 1*3 2*13*3 4*1 51 ()102 1()102 1()102 1()102 1()102x x x x x εεεεε-----=?=?=?=?=? ***124***1244333 (1)() ()()() 111101010222 1.0510x x x x x x εεεε----++=++=?+?+?=? ***123*********123231132143 (2)() ()()() 1111.10210.031100.031385.610 1.1021385.610222 0.215 x x x x x x x x x x x x εεεε---=++=???+???+???≈ **24****24422 *4 33 5 (3)(/)()() 110.0311056.430102256.43056.430 10x x x x x x x εεε---+≈??+??=?= 5计算球体积要使相对误差限为1，问度量半径R 时允许的相对误差限是多少？ 解：球体体积为343V R π= 则何种函数的条件数为 2 3 '4343 p R V R R C V R ππ=== (*)(*)3(*)r p r r V C R R εεε∴≈= 又(*)1r V ε=

数值分析讲义线性方程组的解法

数值分析讲义 第三章线性方程组的解法 §3.0 引言 §3.1 雅可比(Jacobi)迭代法 §3.2 高斯-塞德尔(Gauss-Seidel)迭代法 §3.3 超松驰迭代法§3.7 三角分解法 §3.4 迭代法的收敛性§3.8 追赶法 §3.5 高斯消去法§3.9 其它应用 §3.6 高斯主元素消去法§3.10 误差分析 §3 作业讲评3 §3.11 总结

§3.0 引言 重要性：解线性代数方程组的有效方法在计算数学和科学计算中具有特殊的地位和作用.如弹性力学、电路分析、热传导和振动、以及社会科学及定量分析商业经济中的各种问题. 分类：线性方程组的解法可分为直接法和迭代法两种方法. (a) 直接法：对于给定的方程组，在没有舍入误差的假设下，能在预定的运算次数内求得精确解.最基本的直接法是Gauss消去法，重要的直接法全都受到Gauss消去法的启发.计算代价高. (b) 迭代法：基于一定的递推格式，产生逼近方程组精确解的近似序列.收敛性是其为迭代法的前提，此外，存在收敛速度与误差估计问题.简单实用,诱人.

§3.1 雅可比Jacobi 迭代法 (AX =b ) 1 基本思想: 与解f (x )=0 的不动点迭代相类似,将AX =b 改写为X =BX +f 的形式,建立雅可比方法的迭代格式：X k +1=BX (k )+f ,其中，B 称为迭代矩阵.其计算精度可控,特别适用于求解系数为大型稀疏矩阵(sparse matrices)的方程组. 2 问题: (a) 如何建立迭代格式？ (b) 向量序列{X k }是否收敛以及收敛条件? 3 例题分析: 考虑解方程组??? ??=+--=-+-=--2.453.82102 .72103 21321321x x x x x x x x x (1) 其准确解为X *={1, 1.2, 1.3}. 建立与式(1)相等价的形式： ??? ??++=++=++=84.02.01.083.02.01.072 .02.01.02 13312321x x x x x x x x x (2) 据此建立迭代公式： ?????++=++=++=+++84 .02.01.083.02.01.072.02.01.0)(2)(1)1(3 )(3 )(1)1(23)(2)1(1k k k k k k k k k x x x x x x x x x (3) 取迭代初值0) 0(3 )0(2)0(1===x x x ,迭代结果如下表. JocabiMethodP31.cpp

偏微分方程数值解复习题(2011硕士)

偏微分方程数值解期末复习（2011硕士） 一、考题类型 本次试卷共六道题目，题型及其所占比例分别为： 填空题20%；计算题80% 二、按章节复习内容 第一章 知识点：Euler法、向前差商、向后差商、中心差商、局部截断误差、整体截断误差、相容性、收敛性、阶、稳定性、显格式、隐格式、线性多步法、第一特征多项式、第二特征多项式、稳定多项式、绝对稳定等； 要求: 会辨认差分格式, 判断线性多步法的误差和阶; 第二章 知识点：矩形网格、(正则,非正则)内点、边界点、偏向前(向后,中心)差商、五点差分格式、增设虚点法、积分插值法、线性椭圆型差分格式、极值原理、比较定理、五点差分格式的相容收敛和、稳定性等； 要求: 建立椭圆型方程边值问题的差分格式, 极值原理; 第四章 知识点：最简显格式、最简隐格式、CN格式、双层加权格式、Richardson 格式、网格比、传播因子法(分离变量法) 、传播因子、传播矩阵、谱半径、von Neumann条件、跳点格式、ADI格式、线性椭圆型差分格式、极值原理、比较定理、五点差分格式的相容收敛和稳定性等； 要求: 建立抛物型方程边值问题的差分格式, 计算局部截断误差; 第五章 知识点：左偏心格式、右偏心格式、中心格式、LF格式、LW格式、Wendroff 格式、跳蛙格式、特征线、CFL条件等； 要求: 建立双曲型方程边值问题的差分格式, 计算局部截断误差; 第七章 要求: 会用线性元(线性基)建立常微分方程边值问题的有限元格式

三 练习题 1、 已知显格式21131()22 n n n n u u h f f +++-=-，试证明格式是相容的，并求它的阶。 P39+P41 2、用Taylor 展开原理构造一元函数一阶导数和二阶导数的数值微分公式。 提示：向前、向后和中心差商与一阶导数间关系，二阶中心差商与二阶导数 之间的关系 课件 3、用数值微分方法或数值积分方法建立椭圆型方程 2222(,),(,),u u f x y x y x y ??--=?∈Ω?? :01,01x y Ω≤≤≤≤ 内点差分格式。 P75+课件 4、构造椭圆型方程边值问题的差分格式. P101 (4)题 5、构建一维热传导方程220,(0)u u Lu a a t x ??=-=>??的数值差分格式(显隐格式等)。 参考P132-135相关知识点 6、设有逼近热传导方程22(0)u u Lu a f a const t x ??≡-==>??的带权双层格式 ()()1111111122(1)2k k j j k k k k k k j j j j j j u u a u u u u u u h θθτ++++-+-+-??=-++--+?? 其中[0,1]θ∈，试求其截断误差。并证明当2 1212h a θτ=-时，截断误差的阶最 高阶为24()O h τ+。 P135+P165+课件 7、传播因子法证明抛物型方程22(0)u u Lu a f a const t x ??≡-==>??的最简显隐和六点CN 格式稳定性。 P156+课件 8、对一阶常系数双曲型方程的初边值问题 0,0,0,0,(,0)(),0,(0,)(),0, u u a t T x a t x u x x x u t t t T φψ???+=<≤<<∞>?????=≤<∞??=≤≤?

常微分方程的数值解

实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60;

x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27

最新数值分析课程第五版课后习题答案(李庆扬等)1

第一章 绪论（12） 1、设0>x ，x 的相对误差为δ，求x ln 的误差。 [解]设0*>x 为x 的近似值，则有相对误差为δε=)(*x r ，绝对误差为**)(x x δε=，从而x ln 的误差为δδεε=='=* ****1)()(ln )(ln x x x x x ， 相对误差为* * ** ln ln ) (ln )(ln x x x x r δ εε= = 。 2、设x 的相对误差为2%，求n x 的相对误差。 [解]设*x 为x 的近似值，则有相对误差为%2)(*=x r ε，绝对误差为**%2)(x x =ε，从而n x 的误差为n n x x n x n x x n x x x ** 1 *** %2%2) ()()()(ln * ?=='=-=εε， 相对误差为%2) () (ln )(ln *** n x x x n r == εε。 3、下列各数都是经过四舍五入得到的近似数，即误差不超过最后一位的半个单位，试指出它们是几位有效数字： 1021.1*1=x ，031.0*2=x ，6.385*3=x ，430.56*4=x ，0.17*5 ?=x 。 [解]1021.1*1 =x 有5位有效数字；0031.0* 2=x 有2位有效数字；6.385*3=x 有4位有效数字；430.56* 4 =x 有5位有效数字；0.17*5?=x 有2位有效数字。 4、利用公式（3.3）求下列各近似值的误差限，其中* 4*3*2*1,,,x x x x 均为第3题所给 的数。 （1）* 4*2*1x x x ++； [解]3 334* 4*2*11** *4*2*1*1005.1102 1 10211021)()()()()(----=?=?+?+?=++=? ??? ????=++∑x x x x x f x x x e n k k k εεεε； （2）* 3*2 *1x x x ；

偏微分方程数值解(试题)

偏微分方程数值解试题 1、考虑一维的抛物型方程： 2200, [0,], 0t T (,), (,)(,0)() x x u u x t x u x t u u x t u u x x ππνπ?==??=∈≤≤??=== （1）导出时间离散是一阶向前Euler 格式，空间离散是二阶精度的差分格式； （2）讨论（1）中导出的格式的稳定性； （3）若时间离散为二阶精度的蛙跳格式， 11 2n n n t t u u u t t +-=?-= ?? 空间离散是二阶精度的中心差分，问所导出的格式稳定吗？为什么？ 2、考虑Poission 方程 2(,)1, (,)0, in AB and AD (,)0, in BC and CD u x y x y u n u x y -?=∈Ω ?=?= 其中Ω是图1中的梯形。 使用差分方法来离散该方程。由于梯形的对称性，可以考虑梯形的一半，如图2， 图2 从物理空间到计算区域的几何变换 图1 梯形

为了求解本问题，采用如下方法：将Ω的一半投影到正方形区域?Ω ，然后在?Ω上使用差分方法来离散该方程。在计算区域?Ω 上用N N ?个网格点，空间步长为1/(1)N ξη?=?=-。 （1）引入一个映射T 将原区域Ω（带有坐标,x y ）变换到单位正方形?Ω（带有坐标,ξη）。 同时导出在新区域上的方程和边界条件。 （2）在变换区域，使用泰勒展开导出各导数项在区域内部和边界点上的差分格式。 3、对线性对流方程0 constant >0u u a a t x ??+=??，其一阶迎风有限体积法离散格式为 1?n j u +=?n j u a t x ?-?（?n j u 1?n j u --） （1）写出0a <时的一阶迎风有限体积法的离散格式； （2）写出a 为任意符号的常数的一阶迎风有限体积法的守恒形式。 （3）使用0 u u u t x ??+=??说明一阶迎风有限体积法不是熵保持的格式。 4、对一维Poission 方程 , (0,1) (0)(1)0 x xx u xe x u u ?-=∈? ==? 将[]01，分成(1)n +等分，写出用中心差分离散上述方程的差分格式，并问： （1）该差分格式与原微分方程相容吗？为什么？ （2）该差分格式稳定吗？为什么？ （3）该差分格式是否收敛到原微分方程的解？为什么？ （4）取(1)6n +=，写出该差分格式的矩阵表示。 5、叙述二重网格方法的执行过程，并对一维常微分方程边值问题 2 25, (0,1) (0)(1)0 xx u x x x u u πππ?-=∈? ==?（sin(5)+9sin(15)） 给出限制算子和延拓算子矩阵（以细网格h ：7n =，粗网格2h ：3n =为例）。 6、对一阶波动方程 01(,0)sin(), (0,1)2(0,)(1,)u u t x u x x x u t u t π???+=???? ? =∈?? =??? （1）写出用中心差分进行空间离散，用一阶向后Euler 进行时间离散的差分格式；

偏微分方程数值解法

一、 问题 用有限元方法求下面方程的数值解 2 u u u f t ?-?+=? in (]0,T Ω? 0u = on []0,T ?Ω? ()00,u x u = in Ω 二、 问题分析 第一步 利用Green 公式，求出方程的变分形式 变分形式为：求()()21 00,;u L T H ∈Ω，使得 ()())(2 ,,,,u v u v u v f v t ???+??+= ???? ()10v H ?∈Ω （*） 以及 ()00,u x u =. 第二步 对空间进行离散，得出半离散格式 对区域Ω进行剖分，构造节点基函数，得出有限元子空间：()12,,,h NG V span ???=???，则（*）的Galerkin 逼近为： []0,t T ?∈，求()()1 0,h h u t x V H ∈?Ω，使得 ()()()()() () )(2 ,,,,h h h h h h h d u t v u t v u t v f v dt +??+= h h v V ?∈ （**） 以及()0,0h h u u =，0,h u 为初始条件0u 在h V 中的逼近，设0,h u 为0u 在h V 中的插值. 则0t ?≥，有()()1 N G h i i i u t t ξ? == ∑，0,h u =01 N G i i i ξ?=∑，代人（**）即可得到一常微分方程组. 第三步 进一步对时间进行离散，得到全离散的逼近格式 对 du dt 用差分格式.为此把[]0,T 等分为n 个小区间[]1,i i t t -，其长度1i i T t t t n -?=-= ,n t T =. 这样把求i t 时刻的近似记为i h u ，0 h u 是0u 的近似.这里对（**）采用向后的欧拉格式，即 ()()() () )(2 11 11 1 ,,,,i i i i h h h h h h h i h u u v u v u v f v t ++++-+??+ = ? h h v V ?∈ (***) i=0,1,2…,n-1. 0 h u =0,h u 由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：

数值分析实验2_求解线性方程组直接法

一 实验目的 1．掌握求解线性方程组的高斯消元法及列主元素法； 2. 掌握求解线性方程组的克劳特法； 3. 掌握求解线性方程组的平方根法。 二 实验内容 1．用高斯消元法求解方程组（精度要求为610-=ε）： 1231231 233272212240x x x x x x x x x -+=??-+-=-??-+=? 2．用克劳特法求解上述方程组（精度要求为610-=ε）。 3. 用平方根法求解上述方程组（精度要求为610-=ε）。 4. 用列主元素法求解方程组（精度要求为610-=ε）： 1231231 233432222325x x x x x x x x x -+=??-+-=??--=-? 三 实验步骤（算法）与结果 1. 程序代码（Python3.6）： import numpy as np def Gauss(A,b): n=len(b) for i in range(n-1): if A[i,i]!=0: for j in range(i+1,n): m=-A[j,i]/A[i,i] A[j,i:n]=A[j,i:n]+m*A[i,i:n] b[j]=b[j]+m*b[i] for k in range(n-1,-1,-1): b[k]=(b[k]-sum(A[k,(k+1):n]*b[(k+1):n]))/A[k,k]

print(b) 运行函数： >>> A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=np.float) >>> b=np.array([7,-1,0],dtype=np.float) >>> x=Gauss(A,b) 输出： 结果：解得原方程的解为x1=3.5,x2=-1,x3=-2.25 2 程序代码（Python3.6）： import numpy as np A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=float) L=np.array([[1,0,0],[0,1,0],[0,0,1]],dtype=float) U=np.array([[0,0,0],[0,0,0],[0,0,0]],dtype=float) b=np.array([7,-1,0],dtype=float) y=np.array([0,0,0],dtype=float) x=np.array([0,0,0],dtype=float) def LU(A): n=len(A[0]) i=0 while i