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

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

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

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