- 注册时间
- 2004-9-23
- 最后登录
- 1970-1-1
|
发表于 2007-2-16 20:09:42
|
显示全部楼层
可以分析一下高斯消去法！
我以前用BAS写过一个，你看看我分析的原理吧！你可以移植成LAVA~！
以下是引用的部份原文：
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
下面我们来研究一下这个程序的算法吧.
如,我们有一个方程
x+2y+3z+4w=30①
2x+3y+4z+w=24②
3x+4y+z+2w=22③
4x+y+2z+3w=24④
传统的解法就是将x化成以其它三个未知数表示的代数式,再代入其它三个方程,这样一个一个地消元.这叫作代入法.或者先将每个方程中的某一个未知数的系数化成一样再作加减法来消元.这叫作加减法
现在我们就研究一下用加减法解线情方程(晕!这些方法我初二就会了,不过这里我们要找规律).
第一步
x+2y+3z+4w=30①
y+2z+7w=36 ⑤=①*2-②
2y+8z+10w=68 ⑥=①*3-③
7y+10z+13w=96⑦= ①*4-④
第二步
x+2y+3z+4w=30①
y+2z+7w=36 ⑤
-4z+4w=4 ⑧=⑤*2-⑥
4z+36w=156 ⑨= ⑤*7-⑦
第三步
x+2y+3z+4w=30①
y+2z+7w=36 ⑤
-4z+4w=4 ⑧
160w=640 ⑩=⑧*4-⑨*(-4)
我们现在把方程的编号改一下
x+2y+3z+4w=30⑴
y+2z+7w=36 ⑵
-4z+4w=4 ⑶
160w=640 ⑷
第四步
-160x-320y-480z=-224 ⑸=⑷*4-⑴*160
-160y-320z=-128 ⑹=⑷*7-⑵*160
640z=1920 ⑺=⑷*4-⑶*160
160w=640 ⑷
第五步
1024x-2048y=5120 ⑻=⑺*(-160)-⑸*640
1024y=2048 ⑼=⑺*(-160)-⑹*640
640z=1920 ⑺
160w=640 ⑷
第六步
-1048576x=-1048576 ⑽=⑼*1024-⑻*2048
1024y=2048 ⑼
640z=1920 ⑺
160w=640 ⑷
第七步
x,y,z,w的系数都不等于0,方程有且只有一组解
因此:原方程的解是:
x=1
y=2
z=3
w=4
若我们把上面的方程组的系数提出来进行计算就是
第一步
1 2 3 4 30 ①
0 1 2 7 36 ⑤=①*2-②
0 2 8 10 68 ⑥=①*3-③
0 7 10 13 96 ⑦= ①*4-④
第二步
1 2 3 4 30 ①
0 1 2 7 36 ⑤
0 0 -4 4 4 ⑧=⑤*2-⑥
0 0 4 36 156 ⑨= ⑤*7-⑦
第三步
1 2 3 4 30①
0 1 2 7 36 ⑤
0 0 -4 4 4 ⑧
0 0 0 160 640 ⑩=⑧*4-⑨*(-4)
我们现在把方程的编号改一下
1 2 3 4 30⑴
0 1 2 7 36 ⑵
0 0 -4 4 4 ⑶
0 0 0 160 640 ⑷
第四步
-160 -320 -480 0 -224 ⑸=⑷*4-⑴*160
0 -160 -320 0 -128 ⑹=⑷*7-⑵*160
0 0 640 0 1920 ⑺=⑷*4-⑶*160
0 0 0 160 640 ⑷
第五步
1024 -2048 0 0 5120 ⑻=⑺*(-160)-⑸*640
0 1024 0 0 2048 ⑼=⑺*(-160)-⑹*640
0 0 640 0 1920 ⑺
0 0 0 160 640 ⑷
第六步
-1048576 0 0 0 -1048576 ⑽=⑼*1024-⑻*2048
0 1024 0 0 2048 ⑼
0 0 640 0 1920 ⑺
0 0 0 160 640 ⑷
不难看出,第一至第三步是从上到下每一个方程都依次参加了加减,规律就是第一次循环是用第一个方程的第一个未知数系数乘以第一个后面的每个方程的第一个未知数系数后再一一对应地减去第一个后面的每个方程的第一个未知数系数乘以第一个方程的第一个未知数系数.第二个循环则是用第二个方程的第二个未知数系数乘以第二个后面的每个方程的第二个未知数系数后再一一对应地作减去第二个后面的每个方程的第二个未知数系数乘以第二个方程的第二个未知数系数.以此类推,直到最后一个方程只剩下一个未知数系数为止.这样就相当于把最后一个未知数解出来了.
下面再看第四至第六步:就相反.第一个循环中我们是用最后一个方程的最后一个未知数系数乘以最后一个方程以前所有方程的最后一个未知数系数一一对应地减去最后一个方程以前所有方程的最后一个未知数系数乘以最后一个方程的最后一个未知数系数.以此类推.直到第个方程都只剩下一个未知数系数为此.这样,只要这些未知数系数不为0那么方程就有且只有一组解.如果这些未知数系数当中只要有一个为0且等号右边出为0,则原方程有无数组解.反之则无解!
综上所述,就是解线性方程的数学程序.下面我们再来看一下解线性方程的BAS算法:
由于有多少未知数就要有多少个方程才能解出原方程的解.而且解方程只用到系数.那么这个方程组的所有系数就要用一个二维数组来表示.又因为在计算中要进行一个循环,为了方便,我们就把上面就的那所有的系数用一组三维数组来表示吧.第一步到第三步用一组,第四步到第六步用一组.当然我们也可以用一组四维数组A(1,2,9,10)(最多能计算九元一次方程)来计算.
我们先用下面这个程序输入系数
设M=4(解四元一次方程)
10 FOR I=1 TO M:CLS:A=0:FOR J=1 TO M+1
20 INPUT A(0,0,I,J)
30 FOR K=1 TO M:A=A+A(0,0,I,K):NEXT:IF A=0 THEN BEEP:CLS:20
40 CLS:NEXT:NEXT
(30句是为了防止输入的所有系数都为0)
根据上面的计算方法:下面我们来看 0011表示A(0,0,1,1))
第一步
0111=0011
0112=0012
0113=0013
0114=0014
0115=0015
0121=0021*0021-0011*0021
0122=0022*0021-0011*0022
0123=0023*0021-0011*0023
0124=0024*0021-0011*0024
0125=0025*0021-0011*0025
0131=0031*0031-0011*0031
0132=0032*0031-0011*0032
0133=0033*0031-0011*0033
0134=0034*0031-0011*0034
0135=0035*0031-0011*0035
0141=0041*0041-0011*0041
0142=0042*0041-0011*0042
0143=0043*0041-0011*0043
0144=0044*0041-0011*0044
0145=0045*0041-0011*0045
第二步
0211=0111
0212=0112
0213=0113
0214=0114
0215=0115
0221=0121
0222=0122
0223=0123
0224=0124
0225=0125
0231=0131*0132-0022*0131
0232=0132*0132-0022*0132
0233=0133*0132-0022*0133
0234=0134*0132-0022*0134
0235=0135*0132-0022*0135
0241=0141*0142-0022*0141
0242=0142*0142-0022*0142
0243=0143*0142-0022*0143
0244=0144*0142-0022*0144
0245=0145*0142-0022*0145
第三步
0311=0211
0312=0212
0313=0213
0314=0214
0315=0215
0321=0221
0322=0222
0323=0223
0324=0224
0325=0225
0331=0231
0332=0232
0333=0233
0334=0234
0335=0235
0341=0241*0243-0033*0241
0342=0242*0243-0033*0242
0343=0243*0243-0033*0243
0344=0244*0243-0033*0244
0345=0245*0243-0033*0245
下面我们来转换一下数据,以便进行后面的计算:
1311=0311
1312=0312
1313=0313
1314=0314
1315=0315
1321=0321
1322=0322
1323=0323
1324=0324
1325=0325
1334=0331
1332=0332
1333=0333
1334=0334
1335=0335
1341=0341
1342=0342
1343=0343
1344=0344
1345=0345
第四步
1241=1341
1242=1342
1243=1343
1244=1344
1245=1345
1231=1341*1334-1344*1331
1232=1342*1334-1344*1332
1333=1343*1334-1344*1333
1234=1344*1334-1344*1334
1235=1345*1334-1344*1335
1221=1341*1324-1344*1321
1222=1342*1324-1344*1322
1223=1343*1324-1344*1323
1224=1344*1324-1344*1324
1225=1345*1324-1344*1325
1211=1341*1314-1344*1311
1212=1342*1314-1344*1312
1213=1343*1314-1344*1313
1214=1344*1314-1344*1314
1215=1345*1314-1344*1315
第五步
1141=1241
1142=1242
1143=1243
1144=1244
1145=1245
1131=1231
1132=1232
1133=1233
1134=1234
1135=1235
1121=1231*1223-1233*1221
1122=1232*1223-1233*1222
1123=1233*1223-1233*1223
1124=1234*1223-1233*1224
1125=1235*1223-1233*1225
1111=1231*1213-1233*1211
1112=1232*1213-1233*1212
1113=1233*1213-1233*1213
1114=1234*1213-1233*1214
1115=1235*1213-1233*1215
第六步
1041=1141
1042=1142
1043=1143
1044=1144
1045=1145
1031=1131
1032=1132
1033=1133
1034=1134
1035=1135
1021=1121
1022=1122
1023=1123
1024=1124
1025=1125
1011=1121*1112-1122*1111
1012=1122*1112-1122*1112
1013=1123*1112-1122*1113
1014=1124*1112-1122*1114
1015=1125*1112-1122*1115
好了,
x=1115/1011
y=1125/1022
z=1135/1033
w=1145/1044
只要1011,1022,1033,1044都不等于0方程就有且只有一组解
根据上面的规律.我们就得出下面的程序:
现在我们用M表示未知数个数
I 表示计算的循环次数
J 表示一个方程组的方程序号
K 表示一个方程中从左到右未知数系数的序
100 GRAPH: FOR I=1 TO M-1
110 FOR K=1 TO I: FOR J=1 TO M+1:A(0,I,K,J)=A(0,I-1,K,J):NEXT :NEXT
120 FOR J=I+1 TO M
130 FOR K=1 TO M+1
140 A(0,I,J,K)=A(0,I-1,I,K)*A(0,I-1,J,I)-A(0,I-1,I,I)*A(0,I-1,J,K)
150 NEXT :NEXT :NEXT
180 FOR I=1 TO M: FOR J=1 TO M+1:A(1,M-1,I,J)=A(0,M-1,I,J):NEXT :NEXT
190 FOR I=M-2 TO 0 STEP -1
200 FOR K=I+2 TO M: FOR J=1 TO M+1:A(1,I,K,J)=A(1,I+1,K,J):NEXT :NEXT
210 FOR J=I+1 TO 1 STEP -1
220 FOR K=1 TO M+1
230 A(1,I,J,K)=A(1,I+1,I+2,K)*A(1,I+1,J,I+2)-A(1,I+1,I+2,I+2)*A(1,I+1,J,K)
240 NEXT :NEXT :NEXT /100~240句解方程主体程序/
得出结果后不要忘了对解进行判断,然后再打印到屏屏幕上:
270 A=0:B=0:FOR I=1 TO M:A=A+(A(1,0,I,M+1)<>0):B=B+(A(1,0,I,I)<>0):NEXT /270句判断解的情况/
280 IF A=0 AND B=0 THEN PRINT "原方程有无数个解！":END
290 IF B=0 THEN PRINT "原方程无解！":END
300 FOR I=1 TO M
310 X$(I)=LEFT$((STR$(A(1,0,I,M+1)/A(1,0,I,I))+" "),7):NEXT
315 FOR I=1 TO M:IF VAL(X$(I))<1 AND VAL(X$(I))>0 THEN X$(I)=LEFT$(("0"+X$(I)),7)
316 IF VAL(X$(I))<0 AND VAL(X$(I))>-1 THEN X$(I)=LEFT$(("-0"+RIGHT$(X$(I),6)),7)
320 NEXT :IF M>5 THEN 340
330 FOR I=1 TO M OCATE I,1 RINT A$(I);"=";X$(I);:NEXT :END
340 N=M-INT(M/2)
350 FOR I=1 TO M
360 LOCATE I-(I>N)*N,(I>N)*10+1
370 PRINT A$(I);"=";X$(I);:NEXT :END /280~370句显示解的三种结果,其是300~370句显示原方程有且只有一组时的解/
好了.结合上面所说的全部,我们得出了下面这个程序:
0 :
10 DIM A(1,8,9,10)
20 A$(1)="x":A$(2)="y":A$(3)="z":A$(4)="w":A$(5)="v"
30 A$(6)="u":A$(7)="t":A$(8)="s":A$(9)="r"
50 CLSRINT "请输入未知数个数 (2～9):";:T$=INKEY$
52 M=(T$="n")*2+(T$="m")*3+(T$="g")*4+(T$="h")*5+(T$="j")*6+(T$="t")*7+(T$="y")*8+(T$="u")*9
53 IF T$=CHR$(27) THEN END
55 IF M=0 THEN BEEP:50
60 FOR I=1 TO M:CLS:A=0
70 FOR J=1 TO M+1RINT CHR$(J+96);I;"=";:INPUT "";A(0,0,I,J):NEXT
80 FOR K=1 TO M:A=A+A(0,0,I,K):NEXT :IF A=0 THEN BEEP:CLS:70
90 CLS:NEXT :GRAPH /60~90句键入数据/
97 FOR J=1 TO M: FOR I=1 TO M:IF A(0,0,I,I)=0 THEN 98 ELSE 99
98 FOR K=1 TO M+1:SWAP A(0,0,I,K),A(0,0,I-(I=M)*(M)+1,K):NEXT
99 NEXT :NEXT /97~99句适当调整方程的位置/
100 GRAPH: FOR I=1 TO M-1
110 FOR K=1 TO I: FOR J=1 TO M+1:A(0,I,K,J)=A(0,I-1,K,J):NEXT :NEXT
120 FOR J=I+1 TO M
130 FOR K=1 TO M+1
140 A(0,I,J,K)=A(0,I-1,I,K)*A(0,I-1,J,I)-A(0,I-1,I,I)*A(0,I-1,J,K)
150 NEXT :NEXT :NEXT
180 FOR I=1 TO M: FOR J=1 TO M+1:A(1,M-1,I,J)=A(0,M-1,I,J):NEXT :NEXT
190 FOR I=M-2 TO 0 STEP -1
200 FOR K=I+2 TO M: FOR J=1 TO M+1:A(1,I,K,J)=A(1,I+1,K,J):NEXT :NEXT
210 FOR J=I+1 TO 1 STEP -1
220 FOR K=1 TO M+1
230 A(1,I,J,K)=A(1,I+1,I+2,K)*A(1,I+1,J,I+2)-A(1,I+1,I+2,I+2)*A(1,I+1,J,K)
240 NEXT :NEXT :NEXT /100~240句解方程主体程序/
270 A=0:B=0:FOR I=1 TO M:A=A+(A(1,0,I,M+1)<>0):B=B+(A(1,0,I,I)<>0):NEXT /270句判断解的情况/
280 IF A=0 AND B=0 THEN PRINT "原方程有无数个解！":END
290 IF B=0 THEN PRINT "原方程无解！":END
300 FOR I=1 TO M
310 X$(I)=LEFT$((STR$(A(1,0,I,M+1)/A(1,0,I,I))+" "),7):NEXT
315 FOR I=1 TO M:IF VAL(X$(I))<1 AND VAL(X$(I))>0 THEN X$(I)=LEFT$(("0"+X$(I)),7)
316 IF VAL(X$(I))<0 AND VAL(X$(I))>-1 THEN X$(I)=LEFT$(("-0"+RIGHT$(X$(I),6)),7)
320 NEXT :IF M>5 THEN 340
330 FOR I=1 TO MOCATE I,1:PRINT A$(I);"=";X$(I);:NEXT :END
340 N=M-INT(M/2)
350 FOR I=1 TO M
360 LOCATE I-(I>N)*N,(I>N)*10+1
370 PRINT A$(I);"=";X$(I);:NEXT :END /280~370句显示解的三种结果,其是300~370句显示原方程有且只有一组时的解/
你懂了吗? |
|
发表于 2007-2-16 20:22:13
楼主