|
|
有限区間積分
/1
| f(x)dx
/-1
ガウス・ルジャンドル則の係数(分点、重み)を算出する。
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 LEGENDREPARA(MAXLEVEL,X,W)
PRINT "DIM X(";STR$(MAXLEVEL);"),W(";STR$(MAXLEVEL);")"
FOR I=1 TO INTEGRAL
PRINT "! INPUT PROMPT ";CHR$(34);"下限 ";STR$(I);"=";CHR$(34);":A";STR$(I)
PRINT "! INPUT PROMPT ";CHR$(34);"上限 ";STR$(I);"=";CHR$(34);":B";STR$(I)
NEXT I
FOR I=1 TO INTEGRAL
PRINT "LET A";STR$(I);"=0"
PRINT "LET B";STR$(I);"=1"
NEXT I
FOR J=1 TO INTEGRAL
PRINT "LET U";STR$(J);"=(B";STR$(J);"+A";STR$(J);")/2"
PRINT "LET V";STR$(J);"=(B";STR$(J);"-A";STR$(J);")/2"
NEXT J
PRINT "FOR I=1 TO";MAXLEVEL
PRINT "READ X(I),W(I)"
PRINT "NEXT"
PRINT "LET S=0"
FOR J=1 TO INTEGRAL
PRINT "FOR K";STR$(J);"=1 TO";MAXLEVEL
NEXT J
PRINT "LET S=S+";
FOR J=1 TO INTEGRAL
PRINT "W(K";STR$(J);")*";
NEXT J
PRINT "FUNC(";
FOR I=1 TO INTEGRAL
PRINT "U";STR$(I);"+V";STR$(I);"*X(K";STR$(I);")";
IF I < INTEGRAL THEN PRINT ",";
NEXT I
PRINT ")";
FOR J=1 TO INTEGRAL
PRINT "*V";STR$(J);
NEXT J
PRINT
FOR J=1 TO INTEGRAL
PRINT "NEXT"
NEXT J
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);
IF I < INTEGRAL THEN PRINT ",";
NEXT I
PRINT ")"
PRINT "LET S=1";
FOR I=1 TO INTEGRAL
PRINT "-X";STR$(I);"*X";STR$(I);
NEXT I
PRINT
PRINT "IF S > 0 THEN"
PRINT "LET FUNC=SQR(S)"
PRINT "ELSE"
PRINT "LET FUNC=0"
PRINT "END IF"
PRINT "END FUNCTION"
ELSE
FOR N=2 TO MAXLEVEL
PRINT TAB(8+KETA/2);"分点";TAB(8+KETA*1.5);" 重み"
CALL LEGENDREPARA(N,X,W)
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 LEGENDREPARA(N,A(),W())
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM P(MAXLEVEL),D(MAXLEVEL)
CALL LEGENDREPOLY(N,P)
FOR I=0 TO N
LET P(I)=P(I)/P(N)
NEXT I
FOR I=1 TO N
CALL DERIVATIVE(P,D) !'微分
LET XX=-1 !'初期値
DO
LET X=XX
LET XX=X-HORNER(N,P,X)/HORNER(N,D,X) !'ニュートン法
LOOP UNTIL ABS(HORNER(N,P,XX)) < EPS AND ABS(X-XX) < EPS
LET A(I)=XX !'分点
LET W(I)=WEIGHT(N,XX) !'重み
CALL DIV(P,XX) !'組立除法
NEXT I
END SUB
EXTERNAL SUB LEGENDREPOLY(KK,NEWP()) !'ルジャンドル多項式(係数)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM P(KK+1),OLDP(KK),PP(KK)
LET OLDP(0)=1
LET P(1)=1
FOR I=0 TO KK
LET NEWP(I)=0
NEXT I
FOR K=2 TO KK
FOR J=1 TO K
LET NEWP(J)=NEWP(J)+(2*K-1)/K*P(J-1)
LET NEWP(J-1)=NEWP(J-1)-(K-1)/K*OLDP(J-1)
NEXT J
IF K < KK 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 LEGENDRE(K,X) !'ルジャンドル多項式(値)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
DIM PP(K+1)
LET PP(0)=1
LET PP(1)=X
FOR N=1 TO K-1
LET PP(N+1)=((2*N+1)*X*PP(N)-N*PP(N-1))/(N+1)
NEXT N
LET LEGENDRE=PP(K)
END FUNCTION
EXTERNAL FUNCTION WEIGHT(N,X) !'重み
OPTION ARITHMETIC DECIMAL_HIGH
LET WEIGHT=2*(1-X^2)/(N*LEGENDRE(N-1,X))^2
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) !'組立除法
!'A(N)*X^N+A(N-1)*X^(N-1)+...+A(2)*X^2+A(1)*X+A(0)=(X-P)(C(N-1)*X^(N-1)+...+C(2)*X^2+C(1)*X+C(0))
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
---------------------------------------------------
ルジャンドル多項式表示
有理数モードでお試しください
OPTION BASE 0
PUBLIC NUMERIC MAXLEVEL
LET MAXLEVEL=20 !'次数
DIM P(MAXLEVEL)
PRINT "P(0)=1"
PRINT "P(1)=X"
FOR K=2 TO MAXLEVEL
CALL LEGENDREPOLY(K,P) !'上記参照(OPTION ARITHMETIC DECIMAL_HIGHを外す)
PRINT "P(";STR$(K);")=";
CALL DISPLAY(P)
NEXT K
END
EXTERNAL SUB DISPLAY(A())
FOR N=MAXLEVEL TO 0 STEP -1
IF A(N)<>0 THEN EXIT FOR
NEXT N
IF N > 1 THEN
IF A(N) < 0 THEN PRINT "-";
IF ABS(A(N))<>1 THEN
PRINT STR$(ABS(A(N)));"*X^";STR$(N);
ELSE
PRINT "X^";STR$(N);
END IF
END IF
FOR I=N-1 TO 2 STEP -1
IF A(I)<>0 THEN
IF A(I) < 0 THEN PRINT "-"; ELSE PRINT "+";
IF ABS(A(I))<>1 THEN
PRINT STR$(ABS(A(I)));"*X^";STR$(I);
ELSEIF ABS(A(I))=1 THEN
PRINT "X^";STR$(I);
END IF
END IF
NEXT I
IF A(1)<>0 THEN
IF N > 1 THEN
IF A(1) < 0 THEN PRINT "-"; ELSE PRINT "+";
END IF
IF ABS(A(1))<>1 THEN
PRINT STR$(ABS(A(1)));"*X";
ELSEIF ABS(A(1))=1 THEN
PRINT "X";
END IF
END IF
IF A(0)<>0 THEN
IF A(0) < 0 THEN PRINT "-"; ELSE PRINT "+";
PRINT STR$(ABS(A(0)));
END IF
PRINT
END SUB
|
|