|
|
不等間隔 数値積分
/x(n)
| f(x)dx
/x(1)
積分区間 x(1)～x(n)
PUBLIC NUMERIC MAXLEVEL
OPTION BASE 0
LET MAXLEVEL=10 !'最大次数 (MAXLEVEL > N)
LET N=5 !'データ数
DIM X(N),Y(N),A(MAXLEVEL),B(MAXLEVEL)
FOR I=1 TO N
READ X(I),Y(I) !'y=x^2
NEXT I
CALL LARGRANGE(N,X,Y,A) !'ラグランジュ多項式
CALL INTEGRAL(A,B) !'積分
PRINT HORNER(N,B,X(N))-HORNER(N,B,X(1));X(N)^3/3-X(1)^3/3
!'CALL DERIVATIVE(A,B) !'微分
!'INPUT PROMPT "f'(x) x=":XX
!'PRINT HORNER(N,B,XX);2*XX
DATA 0,0 !'離散値データ X値は等間隔でなくてもいい
DATA 1,1
DATA 3,9
DATA 4,16
DATA 6,36
END
EXTERNAL SUB LARGRANGE(N,X(),Y(),A())
OPTION BASE 0
DIM U(MAXLEVEL),V(MAXLEVEL)
CALL CLR(A)
LET U(1)=1
FOR I = 1 TO N
LET R = Y(I)
CALL CLR(V)
LET V(0)=1
FOR J = 1 TO N
IF I <> J THEN
LET U(0)=-X(J)
CALL MUL(V,U)
LET R = R / (X(I)-X(J))
END IF
NEXT J
CALL SHORTMUL(V,R)
CALL ADD(A,V)
NEXT I
END SUB
EXTERNAL FUNCTION HORNER(N,A(),XX)
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 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 SUB INTEGRAL(A(),B())
FOR I=MAXLEVEL-1 TO 0 STEP -1
LET B(I+1)=A(I)/(I+1)
NEXT I
LET B(0)=0
END SUB
EXTERNAL SUB ADD(A(),B())
FOR I=0 TO MAXLEVEL
LET A(I)=A(I)+B(I)
NEXT I
END SUB
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 SHORTMUL(A(),X)
FOR I=0 TO MAXLEVEL
LET A(I)=A(I)*X
NEXT I
END SUB
EXTERNAL SUB DERIVATIVE(A(),B())
FOR I=MAXLEVEL TO 1 STEP -1
LET B(I-1)=I*A(I)
NEXT I
LET B(MAXLEVEL)=0
END SUB
--------------------------------------------------------------
不等間隔 数値積分 その2
/x(m)
| f(x)dx
/x(1)
積分区間 x(1)～x(m)
LET N=5 !'データ数 条件...MOD(N-1,M-1)=0
LET M=3 !'小区間 間隔数
DIM X(N),Y(N),XX(M),YY(M)
FOR I=1 TO N
READ X(I),Y(I) !'y=x^2
NEXT I
LET S=0
FOR I=1 TO N-M+1 STEP M-1 !'小区間に分割する X(1)～X(3),X(3)～X(5),X(5)～X(7),...
FOR K=0 TO M-1
LET XX(K+1)=X(I+K)
LET YY(K+1)=Y(I+K)
NEXT K
LET S=S+INTEGRAL(XX,YY,M) !'全区間分を足し合わせる
!'LET SS=SS+INTEGRAL2(XX,YY,M) !'下記参照
NEXT I
PRINT S;(X(N)^3-X(1)^3)/3 !' ;SS
!'INPUT PROMPT "f'(x) x=":XX
!'PRINT DERIVATIVE(N,X,Y,XX);2*XX !'下記参照
DATA 0,0 !'離散値データ X値は等間隔でなくてもいい
DATA 1,1
DATA 3,9
DATA 4,16
DATA 6,36
END
EXTERNAL FUNCTION INTEGRAL(X(),Y(),M)
DECLARE EXTERNAL FUNCTION COMB
DIM H(M),XX(M),TEMP(M)
FOR I=1 TO M
LET H(I)=1
LET K=0
FOR J=1 TO M
IF I<>J THEN
LET H(I)=H(I)/(X(I)-X(J))
LET K=K+1
LET XX(K)=X(J)
END IF
NEXT J
LET SIGN=1
FOR J=M TO 1 STEP -1
LET SM=SM+SIGN*X(M)^J/J*COMB(XX,M,M-J,TEMP,1)*H(I)*Y(I)
LET S1=S1+SIGN*X(1)^J/J*COMB(XX,M,M-J,TEMP,1)*H(I)*Y(I)
LET SIGN=-SIGN
NEXT J
NEXT I
LET INTEGRAL=SM-S1
END FUNCTION
EXTERNAL FUNCTION COMB(X(),N,R,A(),K)
IF R=0 THEN
LET S=1
FOR I=1 TO N
IF A(I)=1 THEN LET S=S*X(I)
NEXT I
LET COMB=S
ELSE
FOR I=K TO N-R+1
LET A(I)=1
LET SS=SS+COMB(X,N,R-1,A,I+1)
LET A(I)=0
NEXT I
LET COMB=SS
END IF
END FUNCTION
--------------------------------------------------------------
不等間隔 数値積分 その3(上記参照)
未定係数法
EXTERNAL FUNCTION INTEGRAL2(X(),Y(),N)
DIM A(N,N),B(N)
FOR I=1 TO N
FOR J=1 TO N
LET A(I,J)=X(J)^(I-1)
NEXT J
LET B(I)=(X(N)^I-X(1)^I)/I
NEXT I
MAT A=INV(A)
MAT B=A*B
FOR I=1 TO N
LET S=S+B(I)*Y(I)
NEXT I
LET INTEGRAL2=S
END FUNCTION
--------------------------------------------------------------
不等間隔 数値微分(上記参照)
EXTERNAL FUNCTION DERIVATIVE(N,X(),Y(),XX)
DIM A(N,N),B(N)
FOR I=1 TO N
FOR J=1 TO N
LET A(I,J)=X(J)^I
NEXT J
LET B(I)=I*XX^(I-1)
!'LET B(I)=I*(I-1)*XX^(I-2) !'2階微分
NEXT I
MAT A=INV(A)
MAT B=A*B
FOR I=1 TO N
LET S=S+B(I)*Y(I)
NEXT I
LET DERIVATIVE=S
END FUNCTION
--------------------------------------------------------------
不等間隔 数値微分(定義のみ)
EXTERNAL FUNCTION DIFF(N,X(),Y(),XX)
DIM A(N)
FOR I=1 TO N
LET L=1
LET KK=0
FOR J=1 TO N
IF J<>I THEN
LET KK=KK+1
LET L=L*(X(I)-X(J))
LET A(KK)=X(J)
END IF
NEXT J
LET S1=0
FOR J=1 TO N-1
LET S=1
FOR K=1 TO N-1
IF K<>J THEN LET S=S*(XX-A(K))
NEXT K
LET S1=S1+S
NEXT J
LET SS=SS+S1*Y(I)/L
NEXT I
LET DIFF=SS
END FUNCTION
|
|