当前位置:首页 >> 理学 >>

微分方程数值解法(戴嘉尊 第二版)习题讲解


成都信息工程学院>>精品课程>>微分方程数值解

微分方程数值解 习题解答
杨韧 吴世良(编)

成都信息工程学院 数学学院
二 O 一 O 年四月编写

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解



录?

第一章 常微分方程数值解 ......................................................................3 第二章 抛物型方程的差分方法 ..............................................................8 第三章 椭圆型方程的差分方法 ............................................................16 第四章 双曲型方程的差分方法 ............................................................25

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

第一章 常微分方程数值解
1.解: 由欧拉公式得
2 2 0.1 yn +1 = yn + hf ( xn , yn ) = yn + h( 1+1x2 ? 2 yn ) = yn ? 0.2 yn + 1+ x2
n n

由梯形公式得

yn +1 = yn + 1 2 h[ f ( xn , yn ) + f ( xn +1 , yn +1 )]
2 2 1 1 = yn + 1 2 h[( 1+ x 2 ? 2 yn ) + ( 1+ x 2 ? 2 yn +1 )]
n n+1

= yn ? hy ? hy
2 n

2 n +1

+ h(
1 2

1 2 1+ xn

1 + 1+ x 2 )
n+1

2 2 1 1 1 hyn +1 + yn +1 = yn ? hyn + 2 h( 1+ x 2 + 1+ x 2 )
n n+1

yn +1 =
欧拉公式计算结果
xn

2 1 1 ?1 + 1 + 4h( yn ? hyn +1 2 h( 1+ x 2 + 1+ x 2 ))
n n+1

2h

yn

y ( xn )

y ( xn ) ? yn

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.1000 0.1970 0.2854 0.3609 0.4210 0.4656 0.4957 0.5137 0.5219 0.5227

0 0.0990 0.1923 0.2752 0.3448 0.4000 0.4412 0.4698 0.4878 0.4972 0.5000

0 0.0010 0.0047 0.0102 0.0160 0.0210 0.0244 0.0259 0.0259 0.0247 0.0227

梯形公式计算结果
xn
yn
y ( xn )
y ( xn ) ? yn

0

0

0

0

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.0985 0.1915 0.2742 0.3439 0.3992 0.4406 0.4695 0.4877 0.4973 0.5002

0.0990 0.1923 0.2752 0.3448 0.4000 0.4412 0.4698 0.4878 0.4972 0.5000

0.4758*1.0e-3 0.8291*1.0e-3 0.9894*1.0e-3 0.9580*1.0e-3 0.7886*1.0e-3 0.5523*1.0e-3 0.3097*1.0e-3 0.0988*1.0e-3 0.0640*1.0e-3 0.1773*1.0e-3

2.解:显然, y = e ? x 是原初值问题的准确解。 由梯形公式得
h yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , yn +1 )] 2 h = yn + ( ? yn ? yn +1 ) 2

整理可得: 于是:
y n +1

yn +1 =

2?h yn 2+h
n +1 n +1

2?h ?2?h? ) yn = ? =( ? y n ?1 = 2 + h) ?2+h?
n

2

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

?2?h? y0= ? ? ?2+h?

?2?h? 亦即: y n = ? ? ?2+h?

因为 x = 0 + nh = nh , n =
x

x 2h 1 1 1 ,令 t = ? , =? ? 有 2+h h t 2 h

x x x x ? ? ? ? 2h ? h ? t 2 t 2 = (1 + t ) (1 + t ) y n = ?1 ? ? = (1 + t ) ? 2+h?

从而 lim y n = lim(1 + t )
h→0 t →0

?

x t

? lim(1 + t )
t →0

?

x 2

= e? x

同理可以证明预报-校正法收敛到微分方程的解.

3.解: 局部截断误差:
电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

? Rn +1 = y ( xn +1 ) ? yn +1

= y ( xn ) + ∫ =∫ =∫ =∫ =∫
xn+1 xn xn+1

xn+1

xn

f ( x, y ( x))dx ? [ y ( xn ) + hf ( xn +1 , y ( xn +1 )) ]

f ( x, y ( x))dx ? hf ( xn +1 , y ( xn +1 )) f ( x, y ( x)) ? f ( xn +1 , y ( xn +1 ))dx y′( x) ? y′( xn +1 )dx y′′( xn +1 + θ ( x ? xn +1 ))( x ? xn +1 )dx (0 < θ < 1) ( lagrange中值定理 )
xn+1 xn

xn xn+1

xn xn+1

xn

= y′′( xn +1 + θ ( x ? xn +1 )) ∫

( x ? xn +1 )dx ( xn < x < xn +1 )(积分中值定理 )

2 ′′ =?1 2 h y ( xn +1 + θ ( x ? xn +1 ))

若 M = max ? y′′( xn +1 + θ ( x ? xn +1 )) ,则有
x0 ≤ x ≤ X

2 Rn +1 ≤ 1 2h M

整体截断误差:

ε n +1 ≤ ε n + Rn + ∫

xn+1

xn

f ( xn +1 , y ( xn +1 )) ? f ( xn +1 , yn +1 ) dx

≤ ε n + R + Lh y ( xn +1 ) ? yn +1 , 这里 Rn ≤ R ≤ ε n + R + Lh ε n +1

即有

(1 ? Lh) ε n +1 ≤ ε n + R.
若 Lh < 1 ,则有

εn ≤

R 1 ? Lh 1 ? Lh +

ε n ?1



ε0
(1 ? Lh)
L ( X ? x0 ) 2 n

+

R ? 1 + 1+ ? 1 ? Lh ? 1 ? Lh
L ( X ? x0 ) 2

+

? 1 n ?1 ? (1 ? Lh) ?

≤e

R ? ε0 + ?e Lh ?

? ? 1? ?

4.试分析中点格式
h yn +1 = yn + hf ( xn + h 2 , yn + 2 f ( xn , yn ))

的局部截断误差及整体截断误差。 解: 中点公式的局部截断误差:
电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解
? y ( xn +1 ) ? yn +1 = ∫ xn+1

xn xn+1

h h ? ? f ( x, y( x)) ? f ( xn + 2 , y( xn ) + 2 f ( xn , yn ))? ?dx h h h [ f ( x, y ( x)) ? f ( xn + h 2 , y ( xn + 2 )) + f ( xn + 2 , y ( xn + 2 ))

=∫ =∫

xn

h ? f ( xn + h 2 , y ( xn ) + 2 f ( xn , yn ))]dx xn+1 h h h h [ y′( x) ? y′( x + h 2 )] + [ f ( xn + 2 , y ( xn + 2 )) ? f ( xn + 2 , y ( xn ) + 2 f ( xn , yn ))]dx

xn

所以上式为
en +1 = y′′( xn + θ ) ∫
xn+1 xn

( x ? xn ? h 2 )dx +



xn+1

xn

h h h f ( xn + h 2 , y ( xn + 2 )) ? f ( xn + 2 , y ( xn ) + 2 f ( xn , yn ))dx
2

en +1 ≤ Lh y′′( xn + ? ) ? h8 ≤ LMh3

中点公式的整体截断误差:

ε n +1 = y ( xn +1 ) ? yn +1
= y ( xn ) ? yn + ∫ = y ( xn ) ? yn + ∫
xn+1 xn xn+1 h h ? ? f ( x, y ( x)) ? f ( xn + 2 , yn + 2 f ( xn , yn ))? ?dx h [ f ( x, y ( x)) ? f ( xn + h 2 , y ( xn ) + 2 f ( xn , y ( xn )))

xn

h h h + f ( xn + h 2 , y ( xn ) + 2 f ( xn , y ( xn ))) ? f ( xn + 2 , yn + 2 f ( xn , yn ))]dx

记 F =∫

xn+1

xn

h h h f ( xn + h 2 , y ( xn ) + 2 f ( xn , y ( xn ))) ? f ( xn + 2 , yn + 2 f ( xn , yn ))]dx

注这里 f 满足 Lipschitz 条件, 所以有
h F ≤ Lh | y ( xn ) + h 2 f ( xn , y ( xn )) ? yn ? 2 f ( xn , yn )) | h ≤ Lh[| y ( xn ) ? yn | + h 2 | f ( xn , y ( xn )) ? 2 f ( xn , yn )) |]

≤ Lh(| ε n | + h 2 L | y ( xn ) ? yn |] ≤ Lh(1 + h 2 L) | ε n |

综上有:
| ε n +1 |≤| ε n | + R + Lh(1 + h 2 L) | ε n | ,

即有
2 2 | ε n |≤ (1 + Lh + 1 2 L h ) | ε n ?1 | + R


2 2 n 1 2 2 n ?1 R 1 2 2 ≤ (1 + Lh + 1 ] 2 L h ) | ε 0 | + ? hL ? 1 L2 h 2 [1 + (1 + Lh + 2 L h ) + (1 + Lh + 2 L h )
2

11、解:

≤ e L ( X ? x0 ) | ε 0 | +

R L ( X ? x0 ) (e ? 1) Lh

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

11、解:令 f(x,y)=-y+x+1

y

n +1

=

y

n

+ h(? y + xn + 1) = (1 ? h) y + h( xn + 1) = 0.9 y + xn + 0.1
n n n

(n=0,1,2,…………9)

y

n +1

=

y

n

h + [? y + xn + 1 + (? y + xn +1 + 1)] n n +1 2

=(1-

h h h ) y + ( x n + 1) + (- y + xn +1 + 1 ) n +1 n 2 2 2
n n +1

=0.95 y + 0.05( xn + 1) + 0.05(? y

+ xn+1 + 1)

由初值 y(0)=1 出发,按上述公式计算,结果如下表所示:

x
y

n

0 1

0.1 1.005

0.2 1.019025

0.3 1.041217

0.4 1.070801

0.5 1.107075

0.6 1.149403

0.7 1.197210

0.8 1.249975

0.9 1.307227

1 1.368541

n

12、计算结果与精确解比较。
解:由标准四阶 P-K 方程可得:
0.2 ? ? y n +1 = y n + 6 (k 1 + 2 k 2 + 2 k 3 + k 4) ? k 1 = xn + y n ? ? 0.2 ? = xn + + y + 0.1 k 1 k 2 n ? 2 ? k 3 = xn + 0.1 + y n + 0.1k 2 ? ? k 4 = xn + 0.2 + y n + 0.2 k 3 ?

计算求解得:
n xn 精确解 y(xn) P-K 解 yn 误差 ε 0 0 1 1 0 1 0.2 1.242806 1.242806 10
?6

2 0.4 1.583649 1.583649 1.3 × 10
?5

3 0.6 2.044238 2.044238 2.5 × 10
?5

4 0.8 2.651082 2.651082 4 × 10
?5

5 1 3.436564 3.436564 6.2 × 10
?5

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

第二章 抛物型方程的差分方法
1.解:由
n+1 n n n ?Um = (1? 2r)Um + r(Um ?1 + Um+1 ) ? 0 ?Um = sin(mhπ ) m = 0,1,...M ?U n = U n = 0 n = 1,2,... M ? 0

m = 1,2,..., M ?1

h = 0.1, r =

k 1 = 0.1, k = 0.001, M = = 10 2 h h

得古典显式差分格式
n +1 n n n ?U m = 0.8U m + 0.1(U m ?1 + U m + 1 ) ? 0 m = 0,1, ...10 ?U m = sin(0.1m π ) ?U n = U n = 0 n = 1, 2, ... M ? 0

U

x

0

0.1

0.2

0.3

0.4

t
0 1 0 0 0.3090 0.3060 0.5876 0.5820 0.8090 0.8011 0.9511 0.9417

U

x

0.5 1.0

0.6

0.7

0.8

0.9

t
0 1

1.0000 0 0.9902 0

0.9511

0.8090

0.5878

0.3090

0.9417

0.8011

0.5820

0.3060

2. 解: k = h = 0.1, r =

k = 10 Crank-Nicolson 格式为 h2

n+1 n+1 n+1 n n n ?11Um ? 5(Um ?1 + Um+1 ) = ?9Um + 5(Um?1 + Um+1 ) m = 1,2,...M ?1 ? 0 ?Um = sin(0.1mπ ) m = 0,1,...10 ? Un =Un = 0 n = 1,2,... 10 ? 0

即有

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

? 1 1U 1n + 1 ? 1 1U 2n + 1 ? ? ? ? 1 1U n + 1 8 ? n +1 ? ? 1 1U 9
矩阵形式

? 5 (U 0n + 1 + U 2n + 1 ) = ? 9 U 1n + 5(U 0n + U 2n ) ? 5 (U 1n + 1 + U 3n + 1 ) = ? 9 U 2n + 5(U 1n + U 3n ) ? 5 (U 7n + 1 # + U 9n + 1 ) = ? 9 U 8n + 5(U 7n + U 9n )

? 5 (U 8n + + U 1n0+ 1 ) = ? 9 U 9n + 5(U 8n + U 1n0 )

? 11 ? ?5 ? ? ? ? ? ? ?


?5 11

?5

?5

11 ?5

? ?U 1n +1 ? ? ? 9 ? ? n+` ? ? 5 ? ?U 2 ? ? ? ? ?? ? ?=? ? ? ? ?? n +1 ? ? ? ? ?5 U 8 ? ? n +1 ? ? 11 ? ? ?U 9 ? ? ?

5 ?9

5

5

?9 5

? ?U 1n ? ?? n? ? ?U 2 ? ? ?? ? ? ? ? ?? n? ? ? 5 U8 ?? n? ?9 ? ? ?U 9 ? ?

?? 9 5 ? ?11 ?5 ? ? 5 ?9 5 ? ??5 11 ?5 ? ? ? ? ? ? ? ? ? 5 ?9 5 ?5 11 ?5 ? ? ? ? 5 ?9 5 ?5 11 ?5 ? ? ? ? ? ? B=? A= ? 5 ?9 5 ?5 11 ?5 ? ? ? ? 5 ?9 5 ?5 11 ?5 ? ? ? ? ? ? ? ? 5 ?9 5 ?5 11 ?5 ? ? ? ? 5 ?9 5 ? ?5 11 ?5? ? ? ? ? 5 ? 9? ?5 11? ? ? ? ?
?0.1283 0.0823 0.0528 0.0338 0.0215 0.0136 0.0083 0.0048 0.0022 ? ? 0.0823 0.1811 0.1161 0.0743 0.0473 0.0299 0.0183 0.0105 0.0048 ? ? ? ? 0.0528 0.1161 0.2026 0.1296 0.0826 0.0521 0.0320 0.0183 0.0083 ? ? ? ? 0.0338 0.0743 0.1296 0.2109 0.1344 0.0848 0.0521 0.0299 0.0136? A ?1 = ? 0.0215 0.0473 0.0826 0.1344 0.2131 0.1344 0.0826 0.0473 0.0215? ? ? ? 0.0136 0.0299 0.0521 0.0848 0.1344 0.2109 0.1296 0.0743 0.0338? ? 0.0083 0.0183 0.0320 0.0521 0.0826 0.1296 0.2026 0.1161 0.0528 ? ? ? ? 0.0048 0.0105 0.0183 0.0299 0.0473 0.0743 0.1161 0.1811 0.0823 ? ? 0.0022 0.0048 0.0083 0.0136 0.0215 0.0338 0.0528 0.0823 0.1283? ? ?

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

? sin 0.1π ? ?0.3090? ?sin 0.2π ? ?0.5878? ? ? ? ? ?sin 0.3π ? ?0.8090? ? ? ? ? ?sin 0.4π ? ?0.9511? U 0 = ?sin 0.5π ? = ?1.0000 ? ? ? ? ? ?sin 0.6π ? ?0.9511? ?sin 0.7π ? ?0.8090? ? ? ? ? ?sin 0.8π ? ?0.5878? ?sin 0.9π ? ?0.3090? ? ? ? ?

U1 = A?1BU0 ? -0.7434 ? 0.1646 ? ? 0.1055 ? ? 0.0675 = ? 0.0430 ? ? 0.0271 ? 0.0167 ? ? 0.0096 ? 0.0043 ? 0.1646 0.1055 0.0675 0.0430 0.0271 0.0167 0.0096 0.0043 ? ?0.3090? ?0.1059? ? ? ? ? -0.6378 0.2321 0.1485 0.0947 0.0597 0.0367 0.0210 0.0096? ? ?0.5878? ?0.2015? 0.2321 -0.5948 0.2593 0.1652 0.1042 0.0641 0.0367 0.0167? ?0.8090? ?0.2773? ? ? ?? ? 0.1485 0.2593 -0.5781 0.2688 0.1696 0.1042 0.0597 0.0271? ?0.9511? ?0.3260? 0.0947 0.1652 0.2688 -0.5738 0.2688 0.1652 0.0947 0.0430?i?1.0000? = ?0.3428? ?? ? ? ? 0.0597 0.1042 0.1696 0.2688 -0.5781 0.2593 0.1485 0.0675? ?0.9511? ?0.3260? 0.0367 0.0641 0.1042 0.1652 0.2593 -0.5948 0.2321 0.1055? ?0.8090? ?0.2773? ? ? ?? ? 0.0210 0.0367 0.0597 0.0947 0.1485 0.2321 -0.6378 0.1646? ?0.5878? ?0.2015? ? ?0.3090? ?0.1059? 0.0096 0.0167 0.0271 0.0430 0.0675 0.1055 0.1646 -0.7434? ? ? ? ?

U2 = A?1BU1 ? -0.7434 ? 0.1646 ? ? 0.1055 ? ? 0.0675 = ? 0.0430 ? ? 0.0271 ? 0.0167 ? ? 0.0096 ? 0.0043 ? 0.1646 -0.6378 0.2321 0.1485 0.0947 0.0597 0.0367 0.0043 ? ?0.1059? ?0.0363? ? ? ? ? 0.0096? ? ?0.2015? ?0.0691? 0.0167? ?0.2773? ?0.0951? ? ? ?? ? 0.0271? ?0.3260? ?0.1118? 0.0430?i?0.3428? = ?0.1175? ? ? ?? ? 0.0675? ?0.3260? ?0.1118? 0.1055? ?0.2773? ?0.0951? ? ? ?? ? 0.0210 0.0367 0.0597 0.0947 0.1485 0.2321 -0.6378 0.1646? ?0.2015? ?0.0691? ? ? ? ? 0.0096 0.0167 0.0271 0.0430 0.0675 0.1055 0.1646 -0.7434? ? ?0.1059? ?0.0363? 0.1055 0.2321 -0.5948 0.2593 0.1652 0.1042 0.0641 0.0675 0.1485 0.2593 -0.5781 0.2688 0.1696 0.1042 0.0430 0.0947 0.1652 0.2688 -0.5738 0.2688 0.1652 0.0271 0.0597 0.1042 0.1696 0.2688 -0.5781 0.2593 0.0167 0.0367 0.0641 0.1042 0.1652 0.2593 -0.5948 0.0096 0.0210 0.0367 0.0597 0.0947 0.1485 0.2321

求解结果:
Ux 0 t

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0

0 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0

0.1 0 0.1059 0.2015 0.2773 0.3260 0.3428 0.3260 0.2773 0.2015 0.1059 0 0.2 0 0.0363 0.0691 0.0951 0.1118 0.1175 0.1118 0.0951 0.0691 0.0363 0

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

3.解:古典显格式
n +1 n n n Um ?Um U n ? 2U m + Um +1 = m?1 k h2

右边界

?u ?x

x =1

n n UM + 1 ? U M ?1 ≈ =0 2h

O (h 2 )

n n UM +1 = U M ?1 n 为基础建立古典显格式 在右边界,以 U M

n+1 n n n n n UM = (1? 2r)UM + r(UM ?1 +UM +1) = (1? 2r)UM + 2rUM ?1

古典显式差分格式
n +1 n n n ? Um = (1 ? 2r )U m + r (U m ?1 + U m +1 ) ? 0 ? U m = 1 m = 0,1,..., M ? n n = 1, 2,... ? U0 = 0 n +1 n n ? n = 0,1,... ?U M = (1 ? 2r )U M + 2rU M ?1
0 由 0 层的边界点 U10 = 1算 U
1 10

,U

2 10

,

右边界
m, n+1

其中 r =

k , Mh = 1 h2
m-1, n

截断误差 O (k + h 2 )
m, n m+1, n

4. 解:显式差分格式

(U

n +1 m

?U k

n m

)



U

n m ?1

n ? 2U m +U 2 h

n m +1

? βU

n m

n +1 n n n ? Um = (1 ? 2α r ? β r )U m + α r (U m ?1 + U m + 1 ) ? 0 U = f (m h) m = 0,1, ..., M ? ? mn k = 1, 2... 即 ? U 0 = g 1 ( nk ) ? k ? U n = g ( nk ) k = 1, 2... M h = 1, r = 2 M 2 ? h ?

截断误差 O (k + h 2 ) 方法 I:矩阵方法讨论稳定性
电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

矩阵形式

?U1n+1 ? ?1? 2αr ? βr αr ? ?U1n ? ?αrg1(nk)? ? n+1 ? ? ?? n? ? 0 ? 1? 2αr ? βr αr ?U2 ? ? αr ? ? ?U2 ? ? ? ? ? ? ?? ? ? ? ? =? ? ?? ?+? ? ? ? ? ?? ? ? n? ? ?Un+1 ? ? ? ? 0 ? αr 1? 2αr ? βr αr U ? 8n+1 ? ? ? ? ? 8n ? ? αr 1? 2αr ? βr? ? ?U9 ? ? ? ?U9 ? ? ? ?αrg2 (nk)? ? ?
即 IU n +1 = BU n + e n C=I-1B=B 的特征值

λi = 1 ? 2α r ? β r + 2α r cos(

jπ jπ ) = 1 ? β k ? 4α r sin 2 ( ) M 2M jπ )| ρ ( A) = max | λi |=|1 ? β k ? 4α r sin 2 ( i 2M

即 |1? β k ? 4αr sin2 (

1 jπ 时稳定。 ) |≤ 1时稳定, r ≤ 2α 2M

方法 II:Fourier 级数(Von Neumann 方法)法讨论稳定性

N1 = {0}

a0 = 1 b0 = 1 ? 2α r ? β r b?1 = b1 = α r
βh
2

N 0 = {?1, 0,1}

G(β , k ) = (e0 )?1[(1 ? 2α r ? β r )e0 + α r (eiβ h + e?iβ h )] = 1 ? 2α r ? β r + 2α r cos(β h) = 1 ? 4α r sin 2
| G( β , k ) |=|1 ? 4α r sin 2

? βk

βh
2

? β k |≤|1 ? 4α r sin 2

βh
2

|≤ 1

要求 ?1 ≤ 1 ? 4α r sin 2 β h ≤ 1
2
? 2α r sin 2

βh
2

≤1? αr ≤

1 时稳定。 1即 r≤ 2α 2

5.解:显式差分格式
n +1 n n +1 n +1 (U m ?Um ) U n +1 ? 2U m +Um n +1 +1 = α m ?1 ? βU m k h2

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解
n +1 n +1 n +1 n ? (1 + 2α r + β r )U m ? α r (U m ?1 + U m +1 ) = U m ? 0 m = 0,1, ..., M ?U m = f ( m h ) ? k 即 ? U 0n = g 1 ( nk ) k = 1, 2... r = 2 h ? ? n 1 1 k = 1, 2... h= ,k = ? U M = g 2 ( nk ) M N ?

(1) 差分格式截断误差 O (k + h 2 ) (2) 方法 I. 矩阵方法讨论稳定性 差分格式矩阵形式
? n+1? ? n? ?αrg ((n+1)k) ? U1 ? ? ? ?U1 ? ? 1 ?1+2αr+βr ?? ?αr ? ? ? ? ? ? ? ? n+1? n ? ? ?U ? 1+2αr+βr ?αr ? ?αr ? ?U2 ? 2? ? ? ? ? ?? ? ? ? ? ? ? ?? ? = I ? ?+? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? n 1 n ?αr 1+2αr+βr ?αr ? U + ? U ? ? ? ? ? ? 8 ? ? ?? 8 ? ? α 1 2 α β r r r ? + + ? ? ? ? n+1 ?U ? ?Un? ?αrg ((n+1)k)? ? 9 ? ? 9? ? ? ? ? 2

0

0

即 AUn+1 = IUn + en

C = A?1
jπ ? 1 1 <1 ) |= jπ M 1 + 2α r + β r ? 2α r cos M

ρ ( A ) = | (1 + 2α r + β r ? 2α r cos

无条件稳定。 方法 II. Fourier 级数(Von Neumann 方法)法分析稳定性

N 1 = {? 1, 0,1} a 0 = 1 + 2α r + β r

N 0 = {0} a ? 1 = a1 = ? α r b0 = 1

G ( β , k ) = [(1 + 2α r + β r )e0 ? α r (eiβ h + e? iβ h )]?1 ? e0 = (1 + 2α r + β r ? 2α r cos β h) ?1 = (1 + β r + 2α r (1 ? cos β h)]?1 βh = (1 + 4α r sin 2 + β k ) ?1 2
1 1 + 4α r s in
2

| G ( β , k ) |=

βh
2

+ βk

<1

无条件稳定。
电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

6. 解:差分格式化为
n +1 n n n Um = (1 ? 2iwr )U m + iwr (U m ?1 + U m +1 )

方法 1:矩阵方法
λi = 1 ? 2iwr + 2iw cos
jπ jπ = 1 ? 4iwr sin 2 ( ) M 2M jπ ρ ( A) = max | λi |= 1 + 16w2 r 2 sin 4 ≥1 i 2M

显式差分格式不稳定。 方法 2:Fourier 级数法

G ( β , k ) = ( e 0 ) ?1 [(1 ? 2 iwr ) e 0 + iwr ( e i β h + e ? i β h )] = 1 ? 2 iwr + 2 iwr cos β h = 1 ? 4 iwr sin 2
| G ( β , k ) |2 = 1 + 16 w 2 r 2 sin 4 jπ ≥1 2M

βh
2

显式差分格式不稳定 7. 解:先将差分方程写成方程组形式
+1 n +1 +1 ? U nj +1 ? U nj U nj ? V jn U nj + + U nj ? 1 ? 2U j 1 θ θ (1 ) + ? = ? 2 k k h ? n ? n +1 ?V j = U j

用 Von-Neuman 可以得到增长矩阵的特征值为
λ1,2 =
1 + 2θ ± 1 ? 16θ r sin 2 2(1 + θ + 4r sin 2
βh
2

βh
2

)

=

1 + 2θ ± 1 ? 16θ r sin 2 1 + 2θ + (1 + 8r sin 2
βh
2

βh
2

)

由此可得到稳定性结论如下 a, 当 1 + 2θ < 0 时,格式恒不稳定 c, 当 θ > 0 时, b, 当 ? 1 2 < θ ≤ 0 时, 格式无条件稳定

1 格式条件稳定, 其条件为 r ≤ 16 θ

8. 解: 由 Von-Neuman 方法得增长矩阵
? 1 G=? 2 ? 2 arcsin
βh
2

?2 arcsin 2 1

β h ?1
2

? ? 1 ? i? 2 ? ? ?2 arcsin

βh
2

2 arcsin 2 1

β h ?1
2

? ? ?

求得特征值为

λ1,2 =

1 ? α 2 ± 2iα ,其中 α 2 = arcsin 2 1+ α 2

βh
2

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

即得, λ1,2 = 1 ,故格式无条件稳定. 9. 解: 由 Von-Meuman 方法可以得到(1) (2)的增长因子分别为
(1) G1 ( β , h) = 1 + cτ ? 4r (sin 2 (2) G2 ( β , h) =
β1h
2

+ sin 2

β2h
2

) )

1 1 ? cτ + 4r (sin 2

β1h
2

+ sin 2

β2h
2

由此可知(1)条件稳定 条件为 r ≤ 1 (2)无条件稳定。 4 ;
+1 10. 解:消去 U l ,m 2 便可得到 U ln,m 与 U ln,m 的关系式, 即为
2 ? r 2 τ ?? r 2 τ ? n +1 r 2 2 n ?1 ? δ x ? c ??1 ? δ y ? c ? U lm = δ x δ y U lm 2 ?? 2 2 ? 4 ? 2

?, n + 1

由 Von-Meuman 方法可以得到增长因子

G ( β , h) =
显然无条件稳定。

1h 2h 4r 2 sin 2 β2 sin 2 β2 1h ? τ2 c)(1 + 2r sin 2 (1 + 2r sin 2 β2

β2 h
2

? τ2 c)

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

第三章 椭圆型方程的差分方法
1、解: (1)分别将 U l +1,m +1 ,U l +1,m ?1 , U l ?1,m +1 , U l ?1,m?1 在(l,m)结点处做 Taylor 展开,并
代入
1 (U l +1,m +1 + U l +1,m ?1 + U l ?1, m +1 + U l ?1, m ?1 ? 4U l ,m ) = 0 整理可得其截断误差为 2h 2

R = O(h2 )

(2)方法类似于(1), 同理可得其截断误差为 R = O ( h 4 ) 2、解:作网格
1 25 ln[( + 1) 2 + 12 ] = ln 3 9

2 34 ln[( + 1) 2 + 12 ] = ln 3 9

ln[(0 + 1) 2 +
ln[(0 + 1)2 +

2 13 ] = ln 3 9
2

2

ln[(1 + 1) 2 +
ln[(1 + 1)2 +

2 40 ] = ln 3 9
2

2

1 10 ] = ln 3 9

1 37 ] = ln 3 9

1 16 ln[( +1)2 + 02 ] = ln 3 9

2 25 ln[( + 1)2 + 02 ] = ln 3 9

五点差分格式
10 16 ? ? 4U 1 ? U 2 ? U 3 = ln 9 + ln 9 ? ? ?U + 4U ? U = ln 25 + ln 37 1 2 4 ? 9 9 ? ? ? U + 4U ? U = ln 13 + ln 25 1 3 4 ? 9 9 ? 40 34 ? ?U 2 ? U 3 + 4U 4 = ln + ln 9 9 ? = ln 160 = 0.680725 81 25 × 37 = ln = 2.435345 81 13 × 25 = ln = 1.389376 81 40 × 34 = ln = 2.820791 81

矩阵方程

?4 ? ?1 ? ? ?1 ? ?0

?1 4 0 ?1

?1 0 4 ?1

0 ? ? U 1 ? ? 0.680725 ? ?U ? ? ? ? 1? ? ? 2 ? = ? 2.435345 ? ? AU = K ? 1? ?U 3 ? ?1.389376 ? ?? ? ? ? 4 ? ?U 4 ? ? 2.820791 ?

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

? 0.2917 ? 0.0833 ?1 A =? ? 0.0833 ? ? 0.0417

0.0833 0.2917 0.0417 0.0833

0.0833 0.0417 0.2917 0.0833

0.0417 ? 0.0833 ? ? 0.0833 ? ? 0.2917 ?

? 0.634804 ? ?1.059992 ? ? U = A ?1 K = ? ? 0.798500 ? ? ? ? 1.169821 ?

3、解:首先将网格进行剖分。由于 h=0.5,一共有 9 个网格结点。如图所示的实
红点,其中虚红点是在计算过程中所需要的过度点。

对 U11 建立差分格式得:
?4U11 + U 01 + U 21 + U10 + U12 = 0

(1)
对 U 00 建立差分格式。首先对 U 00 直接应用五点差分格式
U ?10 + U10 + U 0?1 + U 01 ? 4U 00 = 0

(2)
显然, U ?10 与 U 0 ?1 是网格虚点,需要消去,为此,在利用 x=0,y=0 的边界有如下 差分格式
U10 ? U ?10 ? 2hU 00 = 2h

(3)
U 01 ? U 0?1 ? 2hU 00 = ?2h

(4)
电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

联立(2)-(4)得
2U10 + 2U 01 ? (4 + 2h + 2h)U 00 = 0

(5)
类似,对 U 02 U 20 U 22 U10 U 01 U12 U 21 建立差分格式分别为
2U12 + 2U 01 ? (4 + 2h + 2h)U 02 = ?6h

(6)
2U10 + 2U 21 ? (4 ? 2h ? 2h)U 20 = 8h

(7)
2U12 + 2U 21 ? (4 + 2h + 2h)U 22 = 0

(8)
2U11 + U 20 + U 00 ? (4 + 2h)U10 = ?3h

(9)
2U11 + U 02 + U 00 ? (4 + 2h)U 01 = 3h

(10)
2U11 + U 20 + U 22 ? (4 + 2h)U 21 = ?4h

(11)
2U11 + U 02 + U 22 ? (4 + 2h)U12 = 3h

(12)
联立(1), (5), (6)-(12)得
2 0 2 0 0 0 0 0 ? ?(4 + 2h + 2h) ? ?U 00 ? ? 0 ? ? ? ?U ? ? ?3h ? 1 ? (4 + 2 ) 1 0 2 0 0 0 0 h ? ? ? ? 10 ? ? ? ? ?U 20 ? ? 8h ? 0 2 ?(4 + 2h + 2h) 0 0 2 0 0 0 ? ? ?? ? ? 1 0 0 ?(4 + 2h) 2 0 1 0 0 ? ? ?U 01 ? ? 3h ? ? ? ?U11 ? = ? 0 ? 0 1 0 1 ?4 1 0 1 0 ? ? ?? ? ? ? + h 0 0 1 0 2 (4 2 ) 0 0 1 ? ? ?U 21 ? ? 4h ? ? ? ?U ? ? ?6h ? ?(4 + 2h + 2h) 0 0 0 2 0 0 2 0 ? ? ? 02 ? ? ? ? + h 0 0 0 0 2 0 1 (4 2 ) 1 ? ? ?U12 ? ? 3h ? ? ? ?U ? ? 0 ? ?(4 + 2h + 2h) ? 0 0 0 0 0 2 0 2 ? ? ? 11 ? ?

h=0.5, 所以有
电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

? ? ? ? ? ? ? ? ? ? ? ? ? ?

?6 1 0 1 0 0 0 0 0

2 ?5 2 0 1 0 0 0 0

0 1 ?2 0 0 1 0 0 0

2 0 0 ?5 1 0 2 0 0

0 2 0 2 ?4 2 0 2 0

0 0 2 0 1 ?5 0 0 2

0 0 0 1 0 0 ?6 1 0

0 0 0 0 1 0 2 ?5 2

0 ? 0 ? ? 0 ? ? 0 ? 0 ? ? 1 ? 0 ? ? 1 ? ?6? ?

?U ?U ? ?U ? ?U ?U ? ?U ?U ? ?U ?U ?

00 10 20 01 11 21 02 12 11

? ? 0 ? ? ? 1 .5 ? ? ? ? 4 ? ? ? ? 1 .5 ? = ? 0 ? ? ? ? 2 ? ? ?3 ? ? ? ? 1 .5 ? ? ? 0 ?

? ? ? ? ? ? ? ? ? ? ? ? ? ?

解得 U 00 = ?1.1844 , U10 = ?2.1317 , U 20 = ?7.0172 , U 01 = ?1.421 , U11 = ?1.9784
U 21 = ?2.8855 , U 02 = ?0.4655 , U12 = ?1.4752 , U 22 = ?1.4536

T T T 4 、 解 :(1) 将 节 点 以 自 动 顺 序 排 列 U = (U 1T , U 2 , "U q ) , 其 中 ?1

U iT = (U1i ,U 2i ,"U q ?1,i )
则 Dirichlet 问题的差分格式可以写成方程组 AU=b 其中 A,B 即为题中已给形式,b 为由边界条件离散化后已确定的已知向量.

(2)考虑 AU = λAU ,即有
? BU1 + U 2 = λAU1 ? ?" ? ?U i ?1 + BU i + U i +1 = λAU i ?" ? ? ?U q ? 2 + BU q ?1 = λAU q ?1

因 B 为对称阵,故存在正交阵 Q 使
B QT BQ = Λ B = diag (λ1B , λ2B ,", λP ?1 )

用 Q T 左乘以上各式
?QT BU1 + QT U 2 = λ AQT U1 ? ?" ? T T T T ?Q U i ?1 + Q BU i + Q U i +1 = λ AQ U i ?" ? T T T ? ?Q U q ? 2 + Q BU q ?1 = λ AQ U q ?1

记 Vr = QT U r 则

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

?Λ BV1 + V2 = λAV1 ? ?" ? ?Vi ?1 + Λ BVi + Vi +1 = λAVi ?" ? ? ?Vq ? 2 + Λ BVq ?1 = λAVq ?1

选取以上方程组中的每一个小方程组的第 j 个方程得到
?λ jBV j ,1 + V j ,2 = λ AV j ,1 ? ?" ? B ?V j ,i ?1 + λ j V j ,i + V j ,i +1 = λAV j ,i ? ?" ?V j ,q ? 2 + λ jBV j ,q ?1 = λAV j , q ?1 ?

j=1,2…..p-1

这是一齐次三对角线性方程组,若要使其有非零解,则系数行列式值为零,

λ jB ? λA
1

1 λ ? λA 1 % %
B j

%

=0

1

λ jB ? λA

根据三对角阵特征值公式有
π λlA = λ jB + 2cos lq , l = 1, 2,", q ? 1 ,
B λm = ?4 + 2 cos mpπ , m = 1, 2,", p ? 1

于是证得
mπ λl ,m = ?4 + 2(cos lπ p + cos q ) (l = 1, 2, ", p ? 1; m = 1, 2, ", q ? 1)

(3)求解 AU=b 的 Jacobi 迭代矩阵为 G = D ?1 ( L + U ) . 因为 A = D ? L ? U ,所以 G = D ?1 ( D ? A) = I ? D ?1 A .

G 设 x 为 A 的对应于 λlm 的特征向量,于是 A Gx = ( I ? D ?1 A) x = (1 + 1 4 λlm ) x

这说明 x 也是 G 的特征向量,对应的特征值为
G lπ mπ lπ mπ 1 λlm = 1+ 1 4 [ ?4 + 2(cos p + cos q )] = 2 (cos p + cos q ) (l = 1, 2, ", p ? 1; m = 1, 2, ", q ? 1)

因在 (0, π ) 上, cos θ 为减函数,于是
G π π ρ (G ) = max | λlm |= 1 2 (cos p + cos q )

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

5、解:原方程组的系数矩阵为
?1 2 ?2 ? ?. A=? ?1 1 1 ? ? ?2 2 1 ? ?

将矩阵 A 分裂如下形式

A = D ? L ?U ? 1 0 0 ? ? 0 0 0 ? ? 0 ?2 2 ? ? ? ? ? ? =? ?0 1 0 ? ? ? ?1 0 0 ? ? ?0 0 ?1? ? ?0 0 1 ? ? ? ? ?2 ?2 0 ? ? ? ?0 0 0 ? ?

Jacobi 迭代矩阵记为 J ; Guass-Seidel 的迭代矩阵记为 G
? 0 ?2 2 ? ? J = D (L + U ) = ? ? ?1 0 ?1? ? ? ?2 ?2 0 ? ?
?1

,

?1 0 0 ? ?0 ?2 2 ? ?0 ?2 2 ? ? ? ? ? ? ?1 G = ( D ? L) U = ?1 1 0 ? × ? ?0 0 ?1? = ?0 2 ?3? ? ? ?2 2 1? ? ?0 0 0 ? ? ? ?0 0 2 ? ?

?1

通过计算得 ρ ( J ) = 0.5811× 10 ?5 < 1 , ρ (G ) = 2 > 1 。 故原结论成立。

6、解:方法同第 5 题,此处省略。 7、解:原方程组的系数矩阵为
?4 1 0? ?. A=? ?1 6 2 ? ? ?0 2 4? ?

将矩阵 A 分裂如下形式

A = D ? L ?U ? 4 0 0 ? ? 0 0 0 ? ?0 ?1 0 ? ? ? ? ? ? =? ? 0 6 0 ? ? ? ?1 0 0 ? ? ?0 0 ?2 ? ? ?0 0 4? ? ? ? 0 ?2 0 ? ? ? ?0 0 0 ? ?

Jacobi 迭代矩阵记为 J ; Guass-Seidel 的迭代矩阵记为 G , SOR 迭代矩阵记为 S
?0 ?1 1 J = D (L + U ) = ? ?? 6 ? ?0 ?1 4 0 ?1 2 0? ? , ρ ( J ) = 0.4564 ?1 3? 0? ?

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解
?1

?1/ 4 0 ? ?4 0 0? ?0 ?1 0 ? ?0 ? ? ? ? ? ?1 , ?1/ 3 ? G = ( D ? L) U = ? 1 6 0 ? × ?0 0 ?2 ? = ?0 0.0417 ? ? ? ?0 2 4? ? ?0 0 0 ? ? ? ?0 ?0.0208 0.1667 ? ?

ρ (G ) = 0.2083
S = ( D ? wL) ?1 ((1 ? w) D + wU ) 0 0? 0 ? ? ?0.8 ?0.45 0 ? ?4 ? ?3.2 ?1.8 ? ? ? ? ? = ?1.8 6 0 ? × ? 0 ?4.8 ?3.6 ? = ? 0.24 ?0.665 ?0.6 ? ? ? ? 0 ?3.2 ? ? 0 3.6 4 ? ? ? 0 ? ? ? ?0.216 0.5985 ?0.26 ? ?
即有 ρ ( S ) > ρ ( J ) > ρ (G ) 。 故 Guass-Seidel 迭代法的收敛最快,SOR 迭代法收敛最慢
?1

8、解:

五点差分格式

? ?U l ?1,m ? U l +1,m ? U l ,m +1 ? U l ?1,m ? U ? ?U = m ( 3 ? m ) m = 0 , 1, 2 , 3 ? 0 ,m ?U ? 4 , m = 0 m = 0 , 1, 2 , 3 ? π l l = 0 , 1, 2 , 3 , 4 ?U l ,0 = s i n 4 ? ? ? U l , 3 = 0 l = 0 , 1, 2 , 3 , 4
矩阵形式

l ,m ?1

+ 4U

l ,m

= 0

l =1, 2 ,3 m =1,2

AU = K

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

其中

? 4 ?1 0 ? ? B = ? ?1 4 ?1? ? ? ? 0 ?1 4 ? ?

?B A=? ?? I

?I ? B? ?

3π 2 0 0]′ =[2.7071 1 0.7071 2 0 0]′ 4 4 U =[1.0835 0.7405 0.4168 0.8863 0.4616 0.2196]′ K =[2 +sin 1 sin
? 0.2948 ? 0.0932 ? ? 0.0282 A?1 = ? ? 0.0861 ?0.0497 ? ? 0.0195
Jacobi 迭代公式
+ 1) U l(,n = m

π

0.0932 0.0282 0.0861 0.0497 0.0195 ? 0.3230 0.0932 0.0497 0.1056 0.0497 ? ? 0.0932 0.2948 0.0195 0.0497 0.0861? ? 0.0497 0.0195 0.2948 0.0932 0.0282 ? 0.1056 0.0497 0.0932 0.3230 0.0932 ? ? 0.0497 0.0861 0.0282 0.0932 0.2948 ?

1 ) (n) (n) (n) (U l(?n1, m + U l + 1, m + U l , m ? 1 + U l , m + 1 ) l = 1, 2, 3; m = 1, 2 4 1 + 1) (n) ( n + 1) (n) (U l(?n1, m + U l + 1, m + U l , m ? 1 + U l , m + 1 ) l = 1, 2, 3; m = 1, 2 4

Gauss-Seidel 迭代公式
+ 1) U l(,n = m

9、解:

设 f ( y ) = ay + b a ?0 + b = 0 b=0 4 a + b = 64 a = 16

五点差分格式

? ? U l ?1, m ? U l + 1, m ? U l , m ?1 ? U l , m + 1 + 4U ? U 0 ,m = 0 ? ? U l ,0 = 0 ? 3 ? U l = 0, 1,,, 2 3 4 l ,4 = l ? ? ? U 4 , m = 1 6 m m = 1, 2 , 3

l ,m

= 0

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

矩阵形式 其中

AU = K
?B A=? ?? I ? ?0 ?I B ?I 0? ?I ? ? B? ?

? 4 ?1 0 ? ? B=? ? ?1 4 ?1? ? ? 0 ?1 4 ? ?

K =[2.8214 6.3929 11.0000 4.8929 11.0000 21.6071 5.0000 14.1071 31.6786]′

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

第四章 双曲型方程的差分方法
1、解: a = 1 b = 3 x 2 c = x + y代入(4.2)得
dy ? 3 x 2 dx = 0 ? y = x 3 + A 19 = 27 + A A = ?8

特征线 y = x 3 ? 8 特征线通过Γ上点 ( xR , 0) ,即 xR = 2 特征关系 a ( x, y )
du du = c ( x, y ) ? = x + y ? du = ( x + x 3 ? 8) dx dx dx

u=

x2 x4 + ? 8x + B 2 4

代初始条件及Γ上点(2,0), 22 = 2 + 22 ? 8 × 2 + B B = 14

x2 x4 59 或 u = + ? 8 x + 14 , 在点(3,19) u = = 14.75 2 4 4
2、解: a = 1 b =

x c = 2x u

x ? dx = 0 (1) ? ady ? bdx = 0 ? dy ? ?? u ? ?adu ? cdx = 0 ?du ? 2 xdx = 0 (2) ?
由(2)特征关系 u = x 2 + A

i) 在边界 x = 0 y ≥ 0上u |x =0 = 0,特征方程的解为u = x 2

特征线由(1)解出, y = x + B ,沿 (0, yR ) 出发的特征线为 y = x + yR 过点(2,5)的特征线 5 = 2 + yR即y = x + 3
ii) 在边界 y>0,x>0 上, u | y =0 = 0代入u = x 2 + A 。
2 沿点 ( xR , 0) 上,特征方程的解 u = x 2 ? xR

由特征线(1) dy =

2x 2 x ?x
2 2 R

2 dx ? y = x 2 ? xR

2 过点(5,4)的特征线 4 = 52 ? xR xR = 3即y = x 2 ? 9

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

iii) 求初始条件 u | y =0 = x 上的解

在 y=0 上,由 u = x 2 + A 在点 ( xR , 0) 出发的解,
xR = xR 2 + A u = x 2 + ( xR ? xR 2 )

代入特征线 dy =

x x + xR ? x
2 2 R

2 dx ,特征线 y = x 2 + xR ? xR +B

2 通过 ( xR , 0) , 0 = xR + B y = x 2 + xR ? xR ? xR

通过 R(4,0)的特征线上 P(4.05,y)处的 y 及 u 的解析解为

y =

4.05 2 + 4 ? 4 2 ?

4 = 0.0982

u = 4 .0 5 2 + 4 ? 4 2 = 4 .4 0 2 5
iv) 计算近似值,由(4.9.1)和(4.9.2)
(1) ? ? a 0 ( y p ? y 0 ) = b0 ( x p ? x 0 ) ? (1) ? ? a 0 (u p ? u 0 ) = c0 ( x p ? x0 )

由 x0 = 4 y0 = 0 u0 = 4

x p = 4.05 得

?1 ? ( y (1) ? p ? 0) = 2 ? (4.05 ? 4) ? (1) ? ?1 ? (u p ? 4) = 8 ? (4.05 ? 4) y (1) p = 2 × 0.05 = 0.1
又由(4.10)
(2) 1 1 ? ? 2 ( a 0 + a p )( y p ? y 0 ) = 2 ( b0 + b p )( x p ? x 0 ) ?1 (2) 1 ? ? 2 ( a 0 + a p )( u p ? u 0 ) = 2 ( c 0 + c p )( x p ? x 0 ) xp 4.05 a p = 1 bp = = = 1.9308 c p = 2 x p = 8.1 (1) 4.4 up

u (1) p = 8 × 0.05 + 4 = 4.4

? (1 + 1)( y (p2 ) ? 0) = (2 + 1.9308)(4.05 ? 4) ? ? ? (2) ? ? (1 + 1)( u p ? 4) = (8 + 8.1)(4.05 ? 4) y (p2 ) = 1 2 × 3.9308 × 0.05 = 0.0983 u (p2 ) =
1 2

× 16.1 × 0.05 + 4 = 4.4025

3、解: a = x ? y b = u c = x + y 代入 4.71 得

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

?( x ? y )dy ? udx = 0 三个未知量两个方程,不能求解析解。 ? ?( x ? y )du ? ( x + y )dx = 0
求近似解,由 4.9

x0 = 1

y 0 = 0 u 0 = 1 x p = 1 .1 代 入 4 .9 c0 = 0

a 0 = 1 b0 = 0

(1 ) ) ? y (1 = 0 .1 p ? 1 ? ( y p ? 0 ) = 1 ? (1 .1 ? 1) ? ? (1 ) ) u (1 = 1 .1 ? p ? 1 ? ( u p ? 1) = 1 ? (1 .1 ? 1)
) ) = 1 .1 ? 0 .1 = 1 b 0 = u (1 = 1 .1 又 a p = x p ? y (1 p p ) c p = x p + y (1 = 1 .1 + 0 .1 = 1 .2 p

代 入 (4.10)得

(2) 1 1 ? ? 2 ( a 0 + a p )( y p ? y 0 ) = 2 ( b 0 + b p )( x p ? x 0 ) ?1 (2) 1 ? ? 2 ( a 0 + a p )( u p ? u 0 ) = 2 ( c 0 + c p )( x p ? x 0 ) 又 (2) ? ? (1 + 1)( y p ? 0 ) = (1 + 1 .1)(1 .1 ? 1) ? ? (2) ? ? (1 + 1)( u p ? 4 ) = (1 + 1 .2 )(1 .1 ? 1) y (p2 ) = 1 2 × 2 .1 × 0 .1 = 0 .1 0 5

u (p2 ) =

1 2

× 2 .2 × 0 .1 + 1 = 1 .1 1

4、解: 原方程化为 ? ?f? ? 0 ? ?+? 1 ?y ? g ? ? ( x ? 2)( x +1) ? ? ?f? ? ? ? =0, 1? 2 x ( x ? 2)( x +1) ? ?x ? g ?
1 2? x

?1

可以求出 A 的特征值为: λ1 = 对应的左特征向量为: s1 = (1 ?
x +1 3

, λ2 = ? x1 +1 .

x +1) x +1) +1 , ( x ? 2)( ), s2 = ( x3 , ? ( x ? 2)( ), 3 3

因此可以得到方程的正规形式并求解. 5、解: (1) 由 Von Neumann 方法,可以得到增长因子为:

G ( β , Δt ) =

(b cos(

βh
2

)) 2 ? (ar sin

βh
2

)) 2 ? iabr sin( βh)

(b cos(

βh
2

)) 2 + (ar sin

βh
2

)) 2

所以 G ( β , Δ t ) = 1 ,格式无条件稳定。 (2) 利用格式的泰勒展式近似
? ?u ? 2 ? ?u ? 2 a? +b? =c, ? ? ? ? x ? m + 12 ? ? t ? m + 12
电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月
n+ 1 n+ 1

成都信息工程学院>>精品课程>>微分方程数值解

可以得到。
n ? um ? ? ? ln ? i β m h 6、解: (1) 由 Von Neumann 方法,令 ? n ? = ? n ? e ,可以得到增长矩阵为: ? vm ? ?η l ?

? cos( βh) ir sin( βh) ? G ( β , Δt ) = ? ? ir sin( βh) cos(βh) ? ? ? ?,
特征值为: λ ± = cos β h ± ir sin β h ,所以

λ± = 1 ? (1 ? r 2 ) sin 2 βh ,
得到当 r ≤ 1 时,格式稳定。 (2) 同教材 193 页格式(4.81)的讨论。 7、解: 特征值为: λ ± = ± ? p ′ ( v ) ,对应的左特征向量为:
s1 = ( 1 2,2
?1 ? p′( v )

), s2 = ( 1 2,2

1 ? p′( v )

),

则可以相应地写出原方程的正规形式。 8、略 9、解: 取 k = h = 0.2 r = 显式差分格式
n +1 n n n ?1 ?U m =Um ?1 + U m + 1 ? U m ? 0 1 m = 1, 2, 3, 4 ?U m = 8 sin(0.2 mπ ) ? 1 1 ?U m = 18 {sin[0.2( m + 1)π ] + sin[0.2( m ? 1)π ]} ?U n = U n = 0 5 ? 0

k =1 h

U

x

0

0.2

0.4

0.6

0.8

1.0

t
0 0.2 0.4 0.6 0.8 1.0 0 0 0 0 0 0 0.0735 0.0594 0.0227 -0.0227 -0.0594 -0.0735 0.1189 0.0962 0.0367 -0.0367 -0.0962 -0.1189 0.1189 0.0962 0.0367 -0.0367 -0.0962 -0.1189 0.0735 0.0594 0.0227 -0.0227 -0.0594 -0.0735 0 0 0 0 0 0

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月

成都信息工程学院>>精品课程>>微分方程数值解

10、由 Von Neumann 方法可以得到,这里省略。

电子文档制作:成都信息工程学院 数学学院 杨韧 吴世良,2010 年 4 月


相关文章:
微分方程数值解法(戴嘉尊 第二版)习题讲解.pdf
微分方程数值解法(戴嘉尊 第二版)习题讲解 - 成都信息工程学院>>
微分方程数值解法(戴嘉尊_第二版)习题讲解.doc
微分方程数值解法(戴嘉尊_第二版)习题讲解 - 成都信息工程学院>>
微分方程数值解习题课.doc
微分方程数值解习题课_理学_高等教育_教育专区。微分方程 初值问题数值解 习题课 1 一、应用向前欧拉法和改进欧拉法求由如下积分 y ? ? e dt ?t 2 0 x ...
微分方程数值方法习题二.doc
微分方程数值方法习题二 - 微分方程数值方法 常微分方程初值问题习题一 1. 对
微分方程数值解习题(李立康)_图文.doc
微分方程数值解习题(李立康) - 常微分方程习题 习题 1.用 Euler 方法求初值问题 ?u ? ? 1 ? 2 tu ? ?u(0) ? 0 《李立康》 在t ? 1 时的近似解(...
《微分方程数值解法》复习、练习题.doc
微分方程数值解法》复习、练习题 - 第一章 复习题 1、建立差分格式的三个主要步骤(三个离散化) 。 2、差分格式的相容性、收敛性概念。 3、Poisson 方程的...
微分方程数值解法答案.doc
微分方程数值解法答案_理学_高等教育_教育专区。微分...这次考试题目的类型: 20 分的选择题,主要是基本...( )(2) 第二个差分方程的截断误差为 O( ) 2...
偏微分方程的数值解法课后习题答案.doc
微分方程数值解法课后习题答案 - 第二习题答案 q1 q2 a2 q1 a1 q2 a3 q4 a4 q5 a5 q6 第二章 第三章 第四...
偏微分方程的数值解法课后习题答案.doc
微分方程数值解法课后习题答案 - 第二习题答案 q1 a1 q2 a2 q
常微分方程教程_丁同仁(第二版)_习题解答.pdf
微分方程教程_丁同仁(第二版)_习题解答 - 常微分方程教程(第二版)-丁同仁等编-高等教育出版社-参考答案 习题 2-1 判断下列方程是否为恰当方程,并且对恰当...
matlab微分方程例题.ppt
matlab微分方程例题_计算机软件及应用_IT/计算机_专业...求简单微分方程解析解 2、求微分方程数值解. ...
偏微分方程数值解习题解答案.doc
微分方程数值解习题解答案 - q1 q2 2 第二章 第三章 第四章 第五章
偏微分方程的数值解法课后习题答案.doc
微分方程数值解法课后习题答案 - 第二习题答案 六章 q1 q2 a2 六
偏微分方程数值解法答案.doc
微分方程数值解法答案_理学_高等教育_教育专区。偏微分方程数值解法的课后习题答案,是大家想要的!只要你掌握这几个习题,基本上整个书的大部分已经全部理解了 ...
常微分方程教程(丁同仁、李承治第二版)习题解答第6....pdf
微分方程教程(丁同仁、李承治第二版)习题解答第6章6-1 - 特别说明 此
《常微分方程教程》第二版(丁同仁 李承治)课后习题答案....pdf
《常微分方程教程》第二版(丁同仁 李承治)课后习题答案 高等教育出版社_理学_高等教育_教育专区。《常微分方程教程》第二版(丁同仁 李承治)课后习题答案 高等教育...
偏微分方程数值解习题.pdf
微分方程数值解习题_理学_高等教育_教育专区。矩阵、微积分、数值分析 1. 证明 n 阶三对角矩阵 ?? 2 1 ? ? 1 ?2 1 ? ? T=? ? O O O? ? ? ...
偏微分方程数值习题解答.doc
偏微分方程数值习题解答 - 李微分方程数值解习题解答 1-1 如果? '
微分方程数值解问题复习题.pdf
微分方程数值解问题复习题_理学_高等教育_教育专区。...u ? y =0 = x 的解析解在点 (3,19) 之值...微分方程数值解法(戴嘉尊... 29页 1下载券 微分...
偏微分方程数值解习题解答案.doc
微分方程数值解习题解答案_理学_高等教育_教育专区。偏微分方程的一点学习资料,可以看看 第二习题答案 第二章 第三章 第四章 第五章 第六章 q1 q2 a2 ...