您所在位置：主题：我们的实习题目,将Basic转化为PASCAL,请大家帮忙啊!急啊!
gao503
[专家分：0] 发布于 2005-09-27 16:01:00
111 方阵A的列主元Doolittle分解
1. 功能 对方阵A 作LR分解, L为单位下三角阵, R为上三角阵.
2. 数学方法简介[2]
要选列主元,美不先将列的元酸楚在选,选出后作行交换,之后再按先列后行计算.每步的方法步驺为(i=1,2,3…,n-1):
1) aji=aji-∑ljkkki, (1≤k≤i-1) j=i,i+1,…,n
2) 选主元,|Sm|=max|aji|,(i≤j≤n)
3) 交换地I行与地m行,并将SM→Rji
4) lji=aji/sm, j=i+1,i+2,…,n
rij=aij-∑ljkrki, (1≤k≤i-1) j=i+1,i+2,…,n
按上述选主元步驺A作LR分解,必有
|lij|≤1, (i,j=1,2,…,n)
说明L中的元素都是可控的,有利于控制舍入误差的增长,使计算稳定.
我们称列主元Doolittle分解为你LR分解,因为在选列主元时如果作过行交换,则LR≠A,设LR=B,对B作选列主元时行交换的相反的行交换才是A.若是用拟LR分解解线性方程组,对其增广矩阵施行拟LR分解,则最后不必再施行行交换,就可求得线性方程组的解,并且减少了舍入误差,算法稳定.
3. 使用说明
与第110节相同,只是增加了记录选列主元时行交换的工作单元C(N),起先按顺序存放前N个自然数的一维数组.将行交换情况记录下来,最后输出.如例题中C(1)=1,C(2)=4,C(3)=3,C(4)=2,说明选主元时将第二行与第四行作了对调.
10 '* * * * * * * * * * * * * * * * * * * * * * *
20 '* 111 方阵A的列主元doolittle分解 *
30 '* * * * * * * * * * * * * * * * * * * * * * *
40 INPUT "N,E="; N, E
50 DIM A(N, N), C(N)
60 PRINT TAB(3); "EXAMPLE"; TAB(8); "矩阵 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 410
160 PRINT TAB(10); "矩阵 R"
170 PRINT
180 FOR I = 1 TO N
190 FOR J = I TO N
200 PRINT TAB(11 * (J - 1)); USING "###.#####"; A(I, J);
210 NEXT J
220 PRINT
230 NEXT I
240 PRINT
250 FOR I = 1 TO N
260 A(I, I) = 1
270 NEXT I
280 PRINT : PRINT TAB(10); "矩阵 L"
290 FOR I = 1 TO N
300 FOR J = 1 TO I
310 PRINT USING "###.#####"; A(I, J); : PRINT " ";
320 NEXT J
330 PRINT
340 NEXT I
350 PRINT : PRINT
360 FOR I = 1 TO N
370 PRINT "C("; I; ")="; C(I); : PRINT " ";
380 NEXT I
390 DATA 1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256
400 DATA
410 END
500 '子程序 '
510 FOR I = 1 TO N
520 C(I) = I
530 NEXT I
540 Y = 0
550 FOR I = 1 TO N
560 IF ABS(A(I, 1)) <= ABS(Y) THEN 590
570 Y = A(I, 1)
580 I1 = I
590 NEXT I
600 IF ABS(Y) - E > 0 THEN 620
610 GOTO 1090
620 IF I1 = 1 THEN 660
630 FOR J = 1 TO N
640 T = A(I1, J): A(I1, J) = A(1, J): A(1, J) = T
645 Q = C(I1): C(I1) = C(1): C(1) = Q
650 NEXT J
660 FOR I = 2 TO N
670 A(I, 1) = A(I, 1) / A(1, 1)
680 NEXT I
690 P = 0
700 FOR K = 2 TO N
710 FOR I = K TO N
720 FOR R = 1 TO K - 1
730 P = P + A(I,R) * A(R, K)
740 NEXT R
750 A(I, K) = A(I, K) - P
760 P = 0
770 NEXT I
780 '选主元'
790 T = A(K, K): I0 = K
800 FOR I = K + 1 TO N
810 IF ABS(T) >= ABS(A(I, K)) THEN 840
820 T = A(I, K)
830 I0 = I
840 NEXT I
850 IF I0 = K THEN 960
860 '行交换'
870 FOR J = 1 TO N
880 Y = A(I0, K)
890 A(I0, J) = A(K, J)
900 A(K, J) = Y
910 NEXT J
920 Q = C(I0)
930 C(I0) = C(K)
940 C(K) = Q
950 IF ABS (A(K,K))-E<0 THEN 1090
960 FOR I = K + 1 TO N
970 A(I, K) = A(I, K) / A(K, K)
980 NEXT I
990 P = 0
1000 FOR J = K + 1 TO N
1010 FOR R = 1 TO K - 1
1020 P = P + A(K, R) * A(R, J)
1030 NEXT R
1040 A(K, J) = A(K, J) - P
1050 P = 0
1060 NEXT J
1070 NEXT K
1080 W = 0: RETURN
1090 W=1:RETURN
1. 功能 对方阵A 作LR分解, L为单位下三角阵, R为上三角阵.
2. 数学方法简介[2]
要选列主元,美不先将列的元酸楚在选,选出后作行交换,之后再按先列后行计算.每步的方法步驺为(i=1,2,3…,n-1):
1) aji=aji-∑ljkkki, (1≤k≤i-1) j=i,i+1,…,n
2) 选主元,|Sm|=max|aji|,(i≤j≤n)
3) 交换地I行与地m行,并将SM→Rji
4) lji=aji/sm, j=i+1,i+2,…,n
rij=aij-∑ljkrki, (1≤k≤i-1) j=i+1,i+2,…,n
按上述选主元步驺A作LR分解,必有
|lij|≤1, (i,j=1,2,…,n)
说明L中的元素都是可控的,有利于控制舍入误差的增长,使计算稳定.
我们称列主元Doolittle分解为你LR分解,因为在选列主元时如果作过行交换,则LR≠A,设LR=B,对B作选列主元时行交换的相反的行交换才是A.若是用拟LR分解解线性方程组,对其增广矩阵施行拟LR分解,则最后不必再施行行交换,就可求得线性方程组的解,并且减少了舍入误差,算法稳定.
3. 使用说明
与第110节相同,只是增加了记录选列主元时行交换的工作单元C(N),起先按顺序存放前N个自然数的一维数组.将行交换情况记录下来,最后输出.如例题中C(1)=1,C(2)=4,C(3)=3,C(4)=2,说明选主元时将第二行与第四行作了对调.
10 '* * * * * * * * * * * * * * * * * * * * * * *
20 '* 111 方阵A的列主元doolittle分解 *
30 '* * * * * * * * * * * * * * * * * * * * * * *
40 INPUT "N,E="; N, E
50 DIM A(N, N), C(N)
60 PRINT TAB(3); "EXAMPLE"; TAB(8); "矩阵 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 410
160 PRINT TAB(10); "矩阵 R"
170 PRINT
180 FOR I = 1 TO N
190 FOR J = I TO N
200 PRINT TAB(11 * (J - 1)); USING "###.#####"; A(I, J);
210 NEXT J
220 PRINT
230 NEXT I
240 PRINT
250 FOR I = 1 TO N
260 A(I, I) = 1
270 NEXT I
280 PRINT : PRINT TAB(10); "矩阵 L"
290 FOR I = 1 TO N
300 FOR J = 1 TO I
310 PRINT USING "###.#####"; A(I, J); : PRINT " ";
320 NEXT J
330 PRINT
340 NEXT I
350 PRINT : PRINT
360 FOR I = 1 TO N
370 PRINT "C("; I; ")="; C(I); : PRINT " ";
380 NEXT I
390 DATA 1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256
400 DATA
410 END
500 '子程序 '
510 FOR I = 1 TO N
520 C(I) = I
530 NEXT I
540 Y = 0
550 FOR I = 1 TO N
560 IF ABS(A(I, 1)) <= ABS(Y) THEN 590
570 Y = A(I, 1)
580 I1 = I
590 NEXT I
600 IF ABS(Y) - E > 0 THEN 620
610 GOTO 1090
620 IF I1 = 1 THEN 660
630 FOR J = 1 TO N
640 T = A(I1, J): A(I1, J) = A(1, J): A(1, J) = T
645 Q = C(I1): C(I1) = C(1): C(1) = Q
650 NEXT J
660 FOR I = 2 TO N
670 A(I, 1) = A(I, 1) / A(1, 1)
680 NEXT I
690 P = 0
700 FOR K = 2 TO N
710 FOR I = K TO N
720 FOR R = 1 TO K - 1
730 P = P + A(I,R) * A(R, K)
740 NEXT R
750 A(I, K) = A(I, K) - P
760 P = 0
770 NEXT I
780 '选主元'
790 T = A(K, K): I0 = K
800 FOR I = K + 1 TO N
810 IF ABS(T) >= ABS(A(I, K)) THEN 840
820 T = A(I, K)
830 I0 = I
840 NEXT I
850 IF I0 = K THEN 960
860 '行交换'
870 FOR J = 1 TO N
880 Y = A(I0, K)
890 A(I0, J) = A(K, J)
900 A(K, J) = Y
910 NEXT J
920 Q = C(I0)
930 C(I0) = C(K)
940 C(K) = Q
950 IF ABS (A(K,K))-E<0 THEN 1090
960 FOR I = K + 1 TO N
970 A(I, K) = A(I, K) / A(K, K)
980 NEXT I
990 P = 0
1000 FOR J = K + 1 TO N
1010 FOR R = 1 TO K - 1
1020 P = P + A(K, R) * A(R, J)
1030 NEXT R
1040 A(K, J) = A(K, J) - P
1050 P = 0
1060 NEXT J
1070 NEXT K
1080 W = 0: RETURN
1090 W=1:RETURN