|
|
無限区間積分
/∞
| f(x)dx
/-∞
ガウス・エルミート則の係数(分点、重み)を算出する。
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 HERMITEPARA(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"
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);")";
IF I < INTEGRAL THEN PRINT ",";
NEXT I
PRINT ")*EXP(";
FOR I=1 TO INTEGRAL
PRINT "X(I";STR$(I);")^2";
IF I < INTEGRAL THEN PRINT "+";
NEXT I
PRINT ")"
FOR I=1 TO INTEGRAL
PRINT "NEXT"
NEXT I
PRINT "PRINT S,PI^";STR$(INTEGRAL/2)
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);
IF I < INTEGRAL THEN PRINT ",";
NEXT I
PRINT ")"
PRINT "LET FUNC=EXP(";
FOR I=1 TO INTEGRAL
PRINT "-X";STR$(I);"^2";
NEXT I
PRINT ")"
PRINT "END FUNCTION"
ELSE
FOR N=2 TO MAXLEVEL
CALL HERMITEPARA(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 HERMITEPARA(N,A(),W())
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM H(MAXLEVEL),D(MAXLEVEL)
CALL HERMITEPOLY(N,H)
FOR I=0 TO N
LET H(I)=H(I)/H(N)
NEXT I
FOR I=1 TO N
CALL DERIVATIVE(H,D)
LET XX=-15
DO
LET X=XX
LET XX=X-HORNER(N,H,X)/HORNER(N,D,X)
LOOP UNTIL ABS(HORNER(N,H,XX)) < EPS AND ABS(X-XX) < EPS
LET A(I)=XX
LET W(I)=WEIGHT(N,XX)
CALL DIV(H,XX)
NEXT I
END SUB
EXTERNAL SUB HERMITEPOLY(N,NEWP()) !'エルミート多項式(係数)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM P(N+1),OLDP(N)
LET OLDP(0)=1
LET P(1)=2
FOR I=0 TO N
LET NEWP(I)=0
NEXT I
FOR K=2 TO N
FOR J=1 TO K
LET NEWP(J)=NEWP(J)+2*P(J-1)
LET NEWP(J-1)=NEWP(J-1)-2*(K-1)*OLDP(J-1)
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 HERMITE(NN,X) !'エルミート多項式(値)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM H(NN+1)
LET H(0)=1
LET H(1)=2*X
FOR N=1 TO NN-1
LET H(N+1)=2*X*H(N)-2*N*H(N-1)
NEXT N
LET HERMITE=H(NN)
END FUNCTION
EXTERNAL FUNCTION WEIGHT(N,X) !'重み
OPTION ARITHMETIC DECIMAL_HIGH
LET WEIGHT=2^(N+1)*FAC(N)*SQR(PI)/HERMITE(N+1,X)^2
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 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 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 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
|
|