テイラー展開

!数値微分の誤差と計算量 - 再帰呼び出しによるｎ階微分

OPTION BASE 0

LET ep=1e-8 !誤差

LET A=0 !X=0　※マクローリン展開
LET N=9 !N次近似

DIM c(N) !係数
LET c(0)=DF(A,0) !f(a)
!IF ABS(c(0))>ep THEN PRINT c(0) !０以外なら表示する
PRINT c(0)

FOR i=1 TO N
LET c(i)=DF(A,i)/fact(i) !dfn(a)/n!
!IF ABS(c(i))>ep THEN PRINT c(i);"* X ^";i !０以外なら表示する
PRINT c(i);"* X ^";i
PRINT 1/fact(i) !debug debug
NEXT i

END

EXTERNAL FUNCTION F(X) !関数
!LET F=SIN(X)
!LET F=COS(X)
LET F=EXP(X)
END FUNCTION

EXTERNAL FUNCTION DF(x,n) !f(x)のn階微分
IF n>0 THEN !微分係数から求める

!LET h=1/256
!LET DF=(DF(x+h,n-1)-DF(x,n-1))/h

LET h=1/32
LET f1=DF(x+h,n-1)-DF(x-h,n-1) !中点法１次
LET DF=f1/(2*h)
!※LET DF=DF((x+h,n-1)-DF(x-h,n-1))/(2*h) !３点近似

!LET h=1/16
!LET f1=DF(x+h,n-1)-DF(x-h,n-1) !中点法２次
!LET f2=DF(x+2*h,n-1)-DF(x-2*h,n-1)
!LET DF=(2*f1/3-f2/12)/h
!※LET DF=(DF(x-2*h,n-1)-8*DF(x-h,n-1)+8*DF(x+h,n-1)-DF(x+2*h,n-1))/(12*h) !５点近似

!LET h=1/16
!LET f1=DF(x+h,n-1)-DF(x-h,n-1) !中点法３次
!LET f2=DF(x+2*h,n-1)-DF(x-2*h,n-1)
!LET f3=DF(x+3*h,n-1)-DF(x-3*h,n-1)
!LET DF=(f1*3/4-f2*3/20+f3/60)/h
!※LET DF=(DF(x+3*h,n-1)/60-DF(x+2*h,n-1)*3/20+DF(x+h,n-1)*3/4-DF(x-h,n-1)*3/4+DF(x-2*h,n-1)*3/20-DF(x-3*h,n-1)/60)/h !７点近似

!LET h=1/16
!LET f1=DF(x+h,n-1)-DF(x,n-1) !前進法２次　※f0=DF(x,n-1)として計算すれば速くなる
!LET f2=DF(x+2*h,n-1)-DF(x,n-1)
!LET DF=(2*f1-f2/2)/h

!LET h=1/8
!LET f1=DF(x+h,n-1)-DF(x,n-1) !前進法３次　※同上
!LET f2=DF(x+2*h,n-1)-DF(x,n-1)
!LET f3=DF(x+3*h,n-1)-DF(x,n-1)
!LET DF=(3*f1-f2*3/2+f3/3)/h

!LET h=1/16
!LET f1=DF(x-h,n-1)-DF(x,n-1) !後進法２次　※同上
!LET f2=DF(x-2*h,n-1)-DF(x,n-1)
!LET DF=-(2*f1-f2/2)/h

!LET h=1/8
!LET f1=DF(x-h,n-1)-DF(x,n-1) !後進法３次　※同上
!LET f2=DF(x-2*h,n-1)-DF(x,n-1)
!LET f3=DF(x-3*h,n-1)-DF(x,n-1)
!LET DF=-(3*f1-f2*3/2+f3/3)/h

ELSE
LET DF=F(x)
END IF
END FUNCTION

テイラー展開	返事を書くノートメニュー
しばっち <dihjvcfsyu> 2008/02/05 23:17:12 ** この記事は1回修正されてます
PUBLIC NUMERIC H LET H=1/128 LET X=0 LET N=5 PRINT DIFF(X,A) FOR A=1 TO N PRINT DIFF(X,A)/FAC(A);"* (X";SIGN$(X);") ^";A NEXT A END EXTERNAL FUNCTION F(X) LET F=1/(1+X) !'LET F=COS(X)^SIN(X) !'LET F=1/COS(X) !'LET F=TAN(X) !'LET F=EXP(SIN(X)) END FUNCTION EXTERNAL FUNCTION DIFF(X,M) IF M>0 THEN LET DIFF=(DIFF(X+H,M-1)-DIFF(X-H,M-1))/(2H) EXIT FUNCTION END IF IF M=0 THEN LET DIFF=F(X) END FUNCTION EXTERNAL FUNCTION FAC(X) LET S=1 FOR I=2 TO X LET S=SI NEXT I LET FAC=S END FUNCTION EXTERNAL FUNCTION SIGN$(X) IF X<0 THEN LET SIGN$=" + "& STR$(ABS(X)) ELSE LET SIGN$=" - "& STR$(ABS(X)) END FUNCTION 2変数版 PUBLIC NUMERIC H LET H=1/128 LET X=0 LET Y=0 LET M=5 PRINT F(X,Y) FOR N=1 TO M FOR B=0 TO N LET A=N-B PRINT DIFF(X,A,Y,B)/FAC(A)/FAC(B);"* (X";SIGN$(X);") ^";A;"* (Y";SIGN$(Y);") ^";B NEXT B NEXT N END EXTERNAL FUNCTION F(X,Y) LET F=EXP(X)SIN(Y) END FUNCTION EXTERNAL FUNCTION DIFF(X,M,Y,N) IF M>0 THEN LET DIFF=(DIFF(X+H,M-1,Y,N)-DIFF(X-H,M-1,Y,N))/(2H) EXIT FUNCTION END IF IF N>0 THEN LET DIFF=(DIFF(X,M,Y+H,N-1)-DIFF(X,M,Y-H,N-1))/(2H) EXIT FUNCTION END IF IF M=0 OR N=0 THEN LET DIFF=F(X,Y) END FUNCTION EXTERNAL FUNCTION FAC(X) LET S=1 FOR I=2 TO X LET S=SI NEXT I LET FAC=S END FUNCTION EXTERNAL FUNCTION SIGN$(X) IF X<0 THEN LET SIGN$=" + "& STR$(ABS(X)) ELSE LET SIGN$=" - "& STR$(ABS(X)) END FUNCTION

Re: テイラー展開	返事を書くノートメニュー
山中和義 <drdlxujciw> 2008/02/06 22:34:27 ** この記事は1回修正されてます
!数値微分の誤差と計算量 - 再帰呼び出しによるｎ階微分 OPTION BASE 0 LET ep=1e-8 !誤差 LET A=0 !X=0　※マクローリン展開 LET N=9 !N次近似 DIM c(N) !係数 LET c(0)=DF(A,0) !f(a) !IF ABS(c(0))>ep THEN PRINT c(0) !０以外なら表示する PRINT c(0) FOR i=1 TO N LET c(i)=DF(A,i)/fact(i) !dfn(a)/n! !IF ABS(c(i))>ep THEN PRINT c(i);"* X ^";i !０以外なら表示する PRINT c(i);"* X ^";i PRINT 1/fact(i) !debug debug NEXT i END EXTERNAL FUNCTION F(X) !関数 !LET F=SIN(X) !LET F=COS(X) LET F=EXP(X) END FUNCTION EXTERNAL FUNCTION DF(x,n) !f(x)のn階微分 IF n>0 THEN !微分係数から求める !LET h=1/256 !LET DF=(DF(x+h,n-1)-DF(x,n-1))/h LET h=1/32 LET f1=DF(x+h,n-1)-DF(x-h,n-1) !中点法１次 LET DF=f1/(2h) !※LET DF=DF((x+h,n-1)-DF(x-h,n-1))/(2h) !３点近似 !LET h=1/16 !LET f1=DF(x+h,n-1)-DF(x-h,n-1) !中点法２次 !LET f2=DF(x+2h,n-1)-DF(x-2h,n-1) !LET DF=(2f1/3-f2/12)/h !※LET DF=(DF(x-2h,n-1)-8DF(x-h,n-1)+8DF(x+h,n-1)-DF(x+2h,n-1))/(12h) !５点近似 !LET h=1/16 !LET f1=DF(x+h,n-1)-DF(x-h,n-1) !中点法３次 !LET f2=DF(x+2h,n-1)-DF(x-2h,n-1) !LET f3=DF(x+3h,n-1)-DF(x-3h,n-1) !LET DF=(f13/4-f23/20+f3/60)/h !※LET DF=(DF(x+3h,n-1)/60-DF(x+2h,n-1)3/20+DF(x+h,n-1)3/4-DF(x-h,n-1)3/4+DF(x-2h,n-1)3/20-DF(x-3h,n-1)/60)/h !７点近似 !LET h=1/16 !LET f1=DF(x+h,n-1)-DF(x,n-1) !前進法２次　※f0=DF(x,n-1)として計算すれば速くなる !LET f2=DF(x+2h,n-1)-DF(x,n-1) !LET DF=(2f1-f2/2)/h !LET h=1/8 !LET f1=DF(x+h,n-1)-DF(x,n-1) !前進法３次　※同上 !LET f2=DF(x+2h,n-1)-DF(x,n-1) !LET f3=DF(x+3h,n-1)-DF(x,n-1) !LET DF=(3f1-f23/2+f3/3)/h !LET h=1/16 !LET f1=DF(x-h,n-1)-DF(x,n-1) !後進法２次　※同上 !LET f2=DF(x-2h,n-1)-DF(x,n-1) !LET DF=-(2f1-f2/2)/h !LET h=1/8 !LET f1=DF(x-h,n-1)-DF(x,n-1) !後進法３次　※同上 !LET f2=DF(x-2h,n-1)-DF(x,n-1) !LET f3=DF(x-3h,n-1)-DF(x,n-1) !LET DF=-(3f1-f23/2+f3/3)/h ELSE LET DF=F(x) END IF END FUNCTION