自動微分による微分値の計算

!自動微分による微分係数の計算（ボトムアップ方式）

!関数値f、偏微分値∂f/∂xi（微分値f'）の組を、ベクトル[f,∂f/∂x1,∂f/∂x2,…,∂f/∂xn]で表す。
!四則演算などは次のようなベクトル演算で表現できる。
SUB DiffAdd(g(),h(), f()) !加算 g+h
FOR i=1 TO UBOUND(f) !MAT f=g+hでも可
LET f(i)=g(i)+h(i)
NEXT i
END SUB
SUB DiffSub(g(),h(), f()) !減算 g-h
FOR i=1 TO UBOUND(f) !MAT f=g-hでも可
LET f(i)=g(i)-h(i)
NEXT i
END SUB
SUB DiffMul(g(),h(), f()) !乗算 g*h
LET f(1)=g(1)*h(1)
FOR i=2 TO UBOUND(f)
LET f(i)=g(i)*h(1)+g(1)*h(i)
NEXT i
END SUB
SUB DiffDiv(g(),h(), f()) !除算 g/h
LET f(1)=g(1)/h(1)
FOR i=2 TO UBOUND(f)
LET f(i)=(g(i)*h(1)-g(1)*h(i))/h(1)^2
NEXT i
END SUB
SUB DiffCompo(h$,g(), f()) !合成 h(g)
DIM t(2)
CALL FxDFx(h$,g(1), t) !h,h'
LET f(1)=t(1) !h(g)
FOR i=2 TO UBOUND(f)
LET f(i)=t(2)*g(i) !∂h(g)/∂g * ∂g/∂xi
NEXT i
END SUB

SUB DiffSetConst(c, f()) !定数を設定する
MAT f=ZER
LET f(1)=c
END SUB

SUB FxDFx(f$,x, fx()) !関数値f(x),微分値f'(x)を得る
IF UCASE$(left$(f$,2))="X^" THEN !べき乗 x^n
LET n=VAL(right$(f$,LEN(f$)-2))
LET fx(1)=x^n !微分型データ[f,f']
LET fx(2)=n*x^(n-1)
ELSEIF UCASE$(right$(f$,2))="^X" THEN !a^x
LET a=VAL(left$(f$,LEN(f$)-2))
LET fx(1)=a^x
LET fx(2)=a^x*LOG(a)
ELSE
SELECT CASE UCASE$(f$) !関数
CASE "SQR"
LET fx(1)=SQR(x)
LET fx(2)=1/(2*SQR(x))
CASE "SIN" !三角関数
LET fx(1)=SIN(x)
LET fx(2)=COS(x)
CASE "COS"
LET fx(1)=COS(x)
LET fx(2)=-SIN(x)
CASE "TAN"
LET fx(1)=TAN(x)
LET fx(2)=1/COS(x)^2
CASE "ATN"
LET fx(1)=ATN(x)
LET fx(2)=1/(1+x^2)
CASE "EXP" !指数関数
LET fx(1)=EXP(x)
LET fx(2)=EXP(x)
CASE "LOG" !対数関数
LET fx(1)=LOG(x)
LET fx(2)=1/x
CASE ELSE
PRINT f$;"は未定義です。"
STOP
END SELECT
END IF
END SUB

DIM c1(N+1) !定数１
MAT c1=ZER
LET c1(1)=1
!-------------------- ここまでサブルーチン

自動微分による微分値の計算	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/11/25 19:43:53 ** この記事は1回修正されてます
!自動微分による微分係数の計算（ボトムアップ方式） !関数値f、偏微分値∂f/∂xi（微分値f'）の組を、ベクトル[f,∂f/∂x1,∂f/∂x2,…,∂f/∂xn]で表す。 !四則演算などは次のようなベクトル演算で表現できる。 SUB DiffAdd(g(),h(), f()) !加算 g+h FOR i=1 TO UBOUND(f) !MAT f=g+hでも可 LET f(i)=g(i)+h(i) NEXT i END SUB SUB DiffSub(g(),h(), f()) !減算 g-h FOR i=1 TO UBOUND(f) !MAT f=g-hでも可 LET f(i)=g(i)-h(i) NEXT i END SUB SUB DiffMul(g(),h(), f()) !乗算 gh LET f(1)=g(1)h(1) FOR i=2 TO UBOUND(f) LET f(i)=g(i)h(1)+g(1)h(i) NEXT i END SUB SUB DiffDiv(g(),h(), f()) !除算 g/h LET f(1)=g(1)/h(1) FOR i=2 TO UBOUND(f) LET f(i)=(g(i)h(1)-g(1)h(i))/h(1)^2 NEXT i END SUB SUB DiffCompo(h$,g(), f()) !合成 h(g) DIM t(2) CALL FxDFx(h$,g(1), t) !h,h' LET f(1)=t(1) !h(g) FOR i=2 TO UBOUND(f) LET f(i)=t(2)g(i) !∂h(g)/∂g ∂g/∂xi NEXT i END SUB SUB DiffSetConst(c, f()) !定数を設定する MAT f=ZER LET f(1)=c END SUB SUB FxDFx(f$,x, fx()) !関数値f(x),微分値f'(x)を得る IF UCASE$(left$(f$,2))="X^" THEN !べき乗 x^n LET n=VAL(right$(f$,LEN(f$)-2)) LET fx(1)=x^n !微分型データ[f,f'] LET fx(2)=nx^(n-1) ELSEIF UCASE$(right$(f$,2))="^X" THEN !a^x LET a=VAL(left$(f$,LEN(f$)-2)) LET fx(1)=a^x LET fx(2)=a^xLOG(a) ELSE SELECT CASE UCASE$(f$) !関数 CASE "SQR" LET fx(1)=SQR(x) LET fx(2)=1/(2*SQR(x)) CASE "SIN" !三角関数 LET fx(1)=SIN(x) LET fx(2)=COS(x) CASE "COS" LET fx(1)=COS(x) LET fx(2)=-SIN(x) CASE "TAN" LET fx(1)=TAN(x) LET fx(2)=1/COS(x)^2 CASE "ATN" LET fx(1)=ATN(x) LET fx(2)=1/(1+x^2) CASE "EXP" !指数関数 LET fx(1)=EXP(x) LET fx(2)=EXP(x) CASE "LOG" !対数関数 LET fx(1)=LOG(x) LET fx(2)=1/x CASE ELSE PRINT f$;"は未定義です。" STOP END SELECT END IF END SUB DIM c1(N+1) !定数１ MAT c1=ZER LET c1(1)=1 !-------------------- ここまでサブルーチン