|
|
ニュートン・コーツ則 数値積分式
有理数モードでお試しください
OPTION BASE 0
PUBLIC NUMERIC MAXLEVEL
LET MAXLEVEL=10 !'次数
DIM X(1),Y(MAXLEVEL),L(MAXLEVEL)
LET DISPMODE=0 !' 0 or else
LET INTEGRAL=1 !' INTEGRAL >= 1
LET SWITCH=0 !' 0...閉じた公式 else...開いた公式
IF DISPMODE<>0 THEN
IF INTEGRAL > 1 THEN
PRINT "DIM A(";STR$(INTEGRAL);"),B(";STR$(INTEGRAL);"),N(";STR$(INTEGRAL);")"
PRINT "FOR I=1 TO";INTEGRAL
LET A$="(" & CHR$(34) & " & STR$(I) & " & CHR$(34) & ")=" & CHR$(34) & ":"
LET B$="(I)"
ELSE
LET A$="=" & CHR$(34) & ":"
LET B$=""
END IF
LET C$="! INPUT PROMPT " & CHR$(34)
PRINT C$;"下限 ";A$;"A";B$
PRINT C$;"上限 ";A$;"B";B$
PRINT C$;"分割数";A$;"N";B$
PRINT "READ A";B$;",B";B$;",N";B$
IF INTEGRAL > 1 THEN PRINT "NEXT"
FOR I=1 TO INTEGRAL
PRINT "DATA 0,1,10"
NEXT I
FOR I=2 TO MAXLEVEL
PRINT "PRINT INTEGRAL";STR$(I);"(";
IF INTEGRAL=1 THEN
PRINT "A,B,N)"
ELSE
FOR J=1 TO INTEGRAL
PRINT "A(";STR$(J);"),B(";STR$(J);"),";
NEXT J
FOR J=1 TO INTEGRAL
PRINT "N(";STR$(J);")";
IF J < INTEGRAL THEN PRINT ",";
NEXT J
PRINT ")"
END IF
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"
PRINT
END IF
LET X(1)=1
FOR N=1 TO MAXLEVEL-1
FOR I=0 TO N
CALL CLR(Y)
LET P=1
LET Y(0)=1
FOR J=0 TO N
IF I<>J THEN
LET X(0)=-J
CALL MUL(Y,X)
LET P=P*(I-J)
END IF
NEXT J
CALL INTEGRAL(Y)
IF SWITCH=0 THEN
LET L(I)=HORNER(Y,N)/P
ELSE
LET L(I)=(HORNER(Y,N+1)-HORNER(Y,-1))/P
END IF
NEXT I
IF DISPMODE<>0 THEN
PRINT "EXTERNAL FUNCTION INTEGRAL";STR$(N+1);"(";
IF INTEGRAL > 1 THEN
FOR J=1 TO INTEGRAL
PRINT "A";STR$(J);",B";STR$(J);",";
NEXT J
FOR J=1 TO INTEGRAL
PRINT "N";STR$(J);
IF J < INTEGRAL THEN PRINT ",";
NEXT J
ELSE
PRINT "A,B,N";
END IF
PRINT ")"
IF SWITCH=0 THEN LET A$=STR$(N) ELSE LET A$=STR$(N+2)
IF INTEGRAL=1 THEN
PRINT "LET H=(B-A)/N/";A$
PRINT "LET S=0"
PRINT "FOR K=0 TO N-1"
PRINT "LET S=S";
FOR I=0 TO N
IF SWITCH=0 THEN LET B$=STR$(I) ELSE LET B$=STR$(I+1)
IF L(I) < 0 THEN PRINT "-"; ELSE PRINT "+";
PRINT STR$(ABS(L(I)));"*H*FUNC(A+H*(";A$;"*K+";B$;"))";
NEXT I
PRINT
PRINT "NEXT"
ELSE
IF SWITCH=0 THEN
PRINT "DIM R(0 TO ";STR$(N);")"
ELSE
PRINT "DIM R(";STR$(N+1);")"
END IF
FOR I=0 TO N
IF SWITCH=0 THEN LET B$=STR$(I) ELSE LET B$=STR$(I+1)
PRINT "R(";B$;")=";
IF L(I) < 0 THEN PRINT "-";
PRINT STR$(ABS(L(I)))
NEXT I
FOR J=1 TO INTEGRAL
PRINT "LET H";STR$(J);"=(B";STR$(J);"-A";STR$(J);")/N";STR$(J);"/";A$
NEXT J
PRINT "LET S=0"
FOR J=INTEGRAL TO 1 STEP -1
PRINT "FOR K";STR$(J);"=0 TO N";STR$(J);"-1"
NEXT J
FOR J=1 TO INTEGRAL
IF SWITCH=0 THEN
PRINT "FOR J";STR$(J);"=0 TO";N
ELSE
PRINT "FOR J";STR$(J);"=1 TO";N+1
END IF
NEXT J
PRINT "LET S=S+";
FOR J=1 TO INTEGRAL
PRINT "R(J";STR$(J);")*";
NEXT J
FOR I=1 TO INTEGRAL
PRINT "H";STR$(I);"*";
NEXT I
PRINT "FUNC(";
FOR I=1 TO INTEGRAL
PRINT "A";STR$(I);"+H";STR$(I);"*(";A$;"*K";STR$(I);"+J";STR$(I);")";
IF I < INTEGRAL THEN PRINT ",";
NEXT I
PRINT ")"
FOR I=1 TO INTEGRAL*2
PRINT "NEXT"
NEXT I
END IF
PRINT "LET INTEGRAL";STR$(N+1);"=S"
PRINT "END FUNCTION"
ELSE
PRINT "∫(x";STR$(N);",x0)f(x)dx=";
FOR I=0 TO N
IF L(I) < 0 THEN
PRINT "-";
ELSE
IF I > 0 THEN PRINT "+";
END IF
PRINT STR$(ABS(L(I)));"*h*f(x";STR$(I);")";
NEXT I
PRINT
END IF
PRINT
NEXT N
END
EXTERNAL SUB MUL(A(),B())
OPTION BASE 0
DIM C(MAXLEVEL)
FOR I=0 TO MAXLEVEL-1
FOR J=0 TO 1
LET C(I+J)=C(I+J)+A(I)*B(J)
NEXT J
NEXT I
CALL COPY(A,C)
END SUB
EXTERNAL FUNCTION HORNER(A(),XX)
FOR N=MAXLEVEL TO 0 STEP -1
IF A(N)<>0 THEN EXIT FOR
NEXT N
LET Y=A(N)
FOR I=N-1 TO 0 STEP -1
LET Y=Y*XX+A(I)
NEXT I
LET HORNER=Y
END FUNCTION
EXTERNAL SUB COPY(X(),Y())
FOR I=0 TO MAXLEVEL
LET X(I)=Y(I)
NEXT I
END SUB
EXTERNAL SUB CLR(X())
FOR I=0 TO MAXLEVEL
LET X(I)=0
NEXT I
END SUB
EXTERNAL SUB INTEGRAL(A())
OPTION BASE 0
DIM B(MAXLEVEL)
FOR I=MAXLEVEL-1 TO 0 STEP -1
LET B(I+1)=A(I)/(I+1)
NEXT I
CALL COPY(A,B)
END SUB
|
|