109 方阵A的Crout分解

1 功能 对方阵A作LU分解,其中L为下三角阵, U为单位上三角阵.

2 数学方法简介 对奇异方阵A通常是消去主对角线以下的元素,使A变为下三角阵,记未变换的A为A⁽⁰⁾ .

第一步, 若aij⁽⁰⁾ <>0, 令uij =aij⁽⁰⁾ /a11⁽⁰⁾ , j=2,3,…..,n .构造初等方阵

1 u12 u13 . . .u1n

U1= 1

1

.

1

则

1 -u12 -u13 . . . –u1n

U1^(-1) = 1

1

1

用u1^(-1) 右乘A,

a11⁽⁰⁾ 0 0 0

AU1^(-1) = a21⁽⁰⁾ a22⁽¹⁾ a23⁽¹⁾ … a2n⁽¹⁾

… …

an1⁽⁰⁾ an2⁽¹⁾ an3⁽¹⁾ … a nn⁽¹⁾

第二步,若a22⁽¹⁾ <>0, 令u2j =a2j⁽¹⁾ /a22⁽¹⁾ , j=3,4,……,n. 构造初等方阵

1 0 0 … 0

u2^(-1)= 1 -u23 -u2n

1

1

用u2^(-1) 右乘A⁽¹⁾,则

a11⁽⁰⁾ 0 0 … 0

A⁽¹⁾U2^(-1)=A⁽⁰⁾ U1^(-1)U2^(-1)= a21⁽⁰⁾ a22⁽⁰⁾ 0 … 0

a31⁽⁰⁾ a32⁽¹⁾ a33⁽²⁾ a3n⁽²⁾

… … …

an1⁽⁰⁾ an2⁽¹⁾ an3⁽²⁾ ann⁽²⁾

若a33⁽²⁾<>0,可按上述步骤继续往下作,akk^(k-1) <>0,k=1,2,……,n,作完n-1步,A U1^(-1) U2^(-1)…Un-1^(-1)=A^(n-1) =L

则L为下三角阵,令U^(-1) =U1^(-1)U2^(-1)…Un-1^(-1) ,它为单位上三角阵,则

A=A^(n-1) U=LU

这里L为下三角阵,U为单位上三角阵.

为了省去一些工作单元,用紧凑格式编写程序,即将L存放在A的左下方(含对角线),将单位上三角阵U存放在A的右上方(不含对角线元,对角线元全为1不必存放).

上述变换过程的计算方法步骤为,对I=1,2,…,n, 有

lji =aji –∑ljk uki , j=i,i+1,…,n

uij=(aij –∑ljk ukj)/ lii , j=i +1,i+2,…,n

每步先算i列的元素,再算i行的元素.

3 使用说明

输入参数

N 整变量,方阵A的阶数.

A(N,N) N*N个元素的二维实数组,按行存放矩阵A.

输出参数

W 标志,W=0分解完毕,W=1分解进行不下去.

L(N,N) 下三角阵L;

U(N,N) 单位上三角阵U.

10 ‘*************************************’

20 ‘* 109 方阵的A的crout分解 *’

30 ‘*************************'************’

40 INPUT "N,E="; N, E

50 Dim A(N, N)

60 Print Tab(3); "EXAMPLE"; Tab(8); "JUZHEN A": Print

70 For I = 1 To N

80 For J = 1 To N

90 READ A(I, J): Print USING; "####.#"; A(I, J);

100 Next J

110 Print

120 Next I

130 Print: GoSub 500

140 Print Tab(5); "W="; W: Print

150 If (W = 1) Then GoTo 380

160 Print Tab(15); "JUZHEN L"

170 Print

180 For I = 1 To N

190 For J = 1 To I

200 Print USING; "##.#####"; A(I, J);: Print " ";

210 Next J

220 Print

230 Next I

240 For I = 1 To N

250 A(I, I) = 1

260 Next I

270 Print

280 Print Tab(15); "JUZHEN U"

290 Print

300 For I = 1 To N

310 For J = I To N

320 Print Tab(10 * (J - 1)); USING; "##.#####"; A(I, J);

330 Next J

340 Print

350 Next I

360 Data 1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256

370 Data

380 End

500 'SON'

510 If Abs(A(1, 1)) <= E Then GoTo 750

520 For J = 2 To N

530 A(1, J) = A(1, J) / A(1, 1)

540 Next J

550 P = 0

560 For K = 2 To N

570 For I = K To N

580 For R = 1 To K - 1

590 P = P + A(I, R) * A(R, K)

600 Next R

610 A(I, K) = A(I, K) - P

620 P = 0

630 Next I

640 P = 0

650 For J = K + 1 To N

660 For R = 1 To K - 1

670 P = P + A(K, R) * A(R, J)

680 Next R

690 If Abs(A(K, K)) <= E Then GoTo 750

700 A(K, J) = (A(K, J) - P) / A(K, K)

710 P = 0

720 Next J

730 Next K

740 W = 0: Return

750 W = 1: Return