高階数値微分

 投稿者:しばっち  投稿日:2019年 8月30日(金)20時01分53秒
  PUBLIC NUMERIC H
LET H=1/128
LET X=.5
FOR I=1 TO 6
   PRINT DIFF(X,I);DIFFN(X,I)
NEXT I
END

EXTERNAL  FUNCTION F(X)
LET F=X^7
END FUNCTION

EXTERNAL FUNCTION DIFF(X,N) !'再帰呼び出しによる高階数値微分
IF N=0 THEN
   LET DIFF=F(X)
ELSE
!' LET DIFF=(DIFF(X+H,N-1)-DIFF(X-H,N-1))/(2*H) !'3点微分
!' LET DIFF=(-DIFF(X+2*H,N-1)+8*DIFF(X+H,N-1)-8*DIFF(X-H,N-1)+DIFF(X-2*H,N-1))/(12*H) !'5点微分
   LET DIFF=(DIFF(X+3*H,N-1)-9*DIFF(X+2*H,N-1)+45*DIFF(X+H,N-1)-45*DIFF(X-H,N-1)+9*DIFF(X-2*H,N-1)-DIFF(X-3*H,N-1))/(60*H) !'7点微分
   !' LET DIFF=(-3*DIFF(X+4*H,N-1)+32*DIFF(X+3*H,N-1)-168*DIFF(X+2*H,N-1)+672*DIFF(X+H,N-1)-672*DIFF(X-H,N-1)+168*DIFF(X-2*H,N-1)-32*DIFF(X-3*H,N-1)+3*DIFF(X-4*H,N-1))/(840*H) !'9点微分
END IF
END FUNCTION

EXTERNAL FUNCTION DIFFN(X,K)
FOR J=0 TO K
   LET S=S+(-1)^J*COMB(K,J)*F(X+(K/2-J)*H)
NEXT J
LET DIFFN=S/(H^K)
END FUNCTION
 

戻る