|
|
半無限区間積分
/∞
| f(x)dx
/0
ガウス・ラゲール則の係数(分点、重み)を算出する。
1000桁モードを使用し、ラゲール多項式をニュートン法 + 組立除法で解く
OPTION BASE 0
OPTION ARITHMETIC DECIMAL_HIGH
PUBLIC NUMERIC MAXLEVEL,EPS
LET MAXLEVEL=10
DIM X(MAXLEVEL),W(MAXLEVEL)
LET KETA=16
LET EPS=10^(-KETA)
LET DISPMODE=0 !' 0 or else
LET INTEGRAL=1 !' INTEGRAL >= 1
IF DISPMODE<>0 THEN
CALL LAGUERREPARA(MAXLEVEL,X,W)
PRINT "DIM X(";STR$(MAXLEVEL);"),W(";STR$(MAXLEVEL);")"
PRINT "FOR I=1 TO";MAXLEVEL
PRINT "READ X(I),W(I)"
PRINT "NEXT"
FOR I=1 TO INTEGRAL
PRINT "INPUT PROMPT ";CHR$(34);"GAMMA(";
FOR J=1 TO INTEGRAL
PRINT "X";STR$(J);
IF J < INTEGRAL THEN PRINT ",";
NEXT J
PRINT ") X";STR$(I);"=";CHR$(34);":U";STR$(I)
NEXT I
PRINT "LET S=0"
FOR I=1 TO INTEGRAL
PRINT "FOR I";STR$(I);"=1 TO";MAXLEVEL
NEXT I
PRINT "LET S=S+";
FOR I=1 TO INTEGRAL
PRINT "W(I";STR$(I);")*";
NEXT I
PRINT "FUNC(";
FOR I=1 TO INTEGRAL
PRINT "X(I";STR$(I);"),";
NEXT I
FOR I=1 TO INTEGRAL
PRINT "U";STR$(I);
IF I < INTEGRAL THEN PRINT ",";
NEXT I
PRINT ")*EXP(";
FOR I=1 TO INTEGRAL
PRINT "X(I";STR$(I);")";
IF I < INTEGRAL THEN PRINT "+";
NEXT I
PRINT ")"
FOR I=1 TO INTEGRAL
PRINT "NEXT"
NEXT I
PRINT "PRINT S"
FOR I=1 TO MAXLEVEL
PRINT "DATA ";
PRINT USING "##." & REPEAT$("#",KETA):X(I);
PRINT ",";
PRINT USING "#." & REPEAT$("#",KETA) & "^^^^":W(I)
NEXT I
PRINT "END"
PRINT
PRINT "EXTERNAL FUNCTION FUNC(";
FOR I=1 TO INTEGRAL
PRINT "X";STR$(I);",";
NEXT I
FOR I=1 TO INTEGRAL
PRINT "U";STR$(I);
IF I < INTEGRAL THEN PRINT ",";
NEXT I
PRINT ")"
PRINT "FUNC=";
FOR I=1 TO INTEGRAL
PRINT "EXP(-X";STR$(I);")*X";STR$(I);"^(U";STR$(I);"-1)";
IF I < INTEGRAL THEN PRINT "*";
NEXT I
PRINT
PRINT "END FUNCTION"
ELSE
FOR N=2 TO MAXLEVEL
CALL LAGUERREPARA(N,X,W)
PRINT TAB(8+KETA/2);"分点";TAB(8+KETA*1.5);" 重み"
FOR I=1 TO N
PRINT "No.";I;":";
PRINT USING "##." & REPEAT$("#",KETA):X(I);
PRINT " ";
PRINT USING "#." & REPEAT$("#",KETA) & "^^^^":W(I)
NEXT I
NEXT N
END IF
END
EXTERNAL SUB LAGUERREPARA(N,A(),W())
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM LA(MAXLEVEL+1),D(MAXLEVEL)
CALL LAGUERREPOLY(N,LA)
FOR I=0 TO N
LET LA(I)=LA(I)/LA(N)
NEXT I
FOR I=1 TO N
CALL DERIVATIVE(LA,D)
LET XX=0
DO
LET X=XX
LET XX=X-HORNER(N,LA,X)/HORNER(N,D,X)
LOOP UNTIL ABS(HORNER(N,LA,XX)) < EPS AND ABS(X-XX) < EPS
LET A(I)=XX
LET W(I)=WEIGHT(N,XX)
CALL DIV(LA,XX)
NEXT I
END SUB
EXTERNAL SUB LAGUERREPOLY(N,NEWP()) !'ラゲール多項式(係数)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM P(N+1),OLDP(N)
LET OLDP(0)=1
LET P(1)=-1
LET P(0)=1
FOR I=0 TO N
LET NEWP(I)=0
NEXT I
FOR K=2 TO N
FOR J=0 TO K
LET NEWP(J)=NEWP(J)+(2*K-1)*P(J)-(K-1)^2*OLDP(J)
LET NEWP(J+1)=NEWP(J+1)-P(J)
NEXT J
IF K < N THEN
FOR I=0 TO K
LET OLDP(I)=P(I)
LET P(I)=NEWP(I)
LET NEWP(I)=0
NEXT I
END IF
NEXT K
END SUB
EXTERNAL FUNCTION LAGUERRE(NN,X) !'ラゲール多項式(値)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM L(NN+1)
LET L(0)=1
LET L(1)=1-X
FOR N=1 TO NN-1
LET L(N+1)=(2*N+1-X)*L(N)-N*N*L(N-1)
NEXT N
LET LAGUERRE=L(NN)
END FUNCTION
EXTERNAL FUNCTION WEIGHT(N,X)
OPTION ARITHMETIC DECIMAL_HIGH
LET WEIGHT=FAC(N)^2/(X*LAGUERREDIFF(N,X)^2)
END FUNCTION
EXTERNAL FUNCTION LAGUERREDIFF(N,X)
OPTION ARITHMETIC DECIMAL_HIGH
LET LAGUERREDIFF=(LAGUERRE(N+1,X)-(N+1-X)*LAGUERRE(N,X))/X
END FUNCTION
EXTERNAL FUNCTION FAC(X)
OPTION ARITHMETIC DECIMAL_HIGH
LET S=1
FOR I=2 TO X
LET S=S*I
NEXT I
LET FAC=S
END FUNCTION
EXTERNAL SUB DERIVATIVE(A(),B())
OPTION ARITHMETIC DECIMAL_HIGH
FOR I=MAXLEVEL TO 1 STEP -1
LET B(I-1)=I*A(I)
NEXT I
LET B(MAXLEVEL)=0
END SUB
EXTERNAL FUNCTION HORNER(N,A(),X)
OPTION ARITHMETIC DECIMAL_HIGH
LET Y=A(N)
FOR I=N-1 TO 0 STEP -1
LET Y=Y*X+A(I)
NEXT I
LET HORNER=Y
END FUNCTION
EXTERNAL SUB DIV(A(),P)
OPTION BASE 0
OPTION ARITHMETIC DECIMAL_HIGH
DIM C(MAXLEVEL)
FOR I=MAXLEVEL TO 1 STEP -1
LET C(I-1)=A(I)+C(I)*P
NEXT I
CALL COPY(A,C)
END SUB
EXTERNAL SUB COPY(X(),Y())
OPTION ARITHMETIC DECIMAL_HIGH
FOR I=0 TO MAXLEVEL
LET X(I)=Y(I)
NEXT I
END SUB
|
|