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

続き

!●問題
!　関数 f(x1,x2,x3)=sin(1.5*(x1+x2))/(x2*x3) のとき、x1=x2=x3=1における関数値、偏微分係数は？

!●数値微分の場合（前進差分）
DEF f(x1,x2,x3)=SIN(1.5*(x1+x2))/(x2*x3) !関数を定義する
LET N=3 !変数の数

LET xx1=1
LET xx2=1
LET xx3=1
LET d=1E-8 !⊿
PRINT f(xx1,xx2,xx3), !f
PRINT (f(xx1+d,xx2,xx3)-f(xx1,xx2,xx3))/d, !∂f/∂x1
PRINT (f(xx1,xx2+d,xx3)-f(xx1,xx2,xx3))/d, !∂f/∂x2
PRINT (f(xx1,xx2,xx3+d)-f(xx1,xx2,xx3))/d !∂f/∂x3
PRINT

!●自動微分の場合

!x1=[1,1,0,0]　※[x1, ∂x1/∂x1, ∂x1/∂x2, ∂x1/∂x3]
!x2=[1,0,1,0]
!x3=[1,0,0,1]
DIM x1(N+1),x2(N+1),x3(N+1) !変数の定義
DATA 1,1,0,0
MAT READ x1
DATA 1,0,1,0
MAT READ x2
DATA 1,0,0,1
MAT READ x3

!関数値、偏微分値がすぐにわかるような関数同士の四則演算、合成の形に分解する。
!　v1=x1+x2
!　v2=1.5*v1
!　v3=sin(v2)
!　v4=x2*x3
!　v5=v3/v4
!
!※式を手コンパイルで求める。
!多言語（C++など）では、オペランドのオーバーロード（演算子の多重定義）機能で実現する。

!v1=x1+x2=[1,1,0,0]+[1,0,1,0]=[2,1,1,0]
DIM v1(N+1)
CALL DiffAdd(x1,x2, v1)
MAT PRINT v1

!v2=1.5*v1=[1.5,0,0,0]*[2,1,1,0]=[3,1.5,1.5,0]
DIM c15(N+1) !定数1.5
CALL DiffSetConst(1.5, c15)
DIM v2(N+1)
CALL DiffMul(c15,v1, v2)
!!!MAT v2=1.5*v1でも可
MAT PRINT v2

!v3=sin(v2)=[sin(v2),(∂sin(v2)/∂v2)*(∂v2/∂x1),(∂sin(v2)/∂v2)*(∂v2/∂x2),(∂sin(v2)/∂v2)*(∂v2/∂x3)]=[sin(3),1.5*cos(3),1.5*cos(3),0]
DIM v3(N+1)
!LET v3(1)=SIN(FxVal(v2))
!LET dfx=DFxVal("SIN",v2) !※∂sin(v2)/∂v2
!LET v3(2)=dfx*1.5 !※∂v2/∂x1=∂{1.5*(x1+x2)}/∂x1=1.5
!LET v3(3)=dfx*1.5
!LET v3(4)=dfx*0
CALL DiffCompo("SIN",v2, v3)

!v4=x2*x3=[1,0,1,0]*[1,0,0,1]=[1,0,1,1]
DIM v4(N+1)
CALL DiffMul(x2,x3, v4)
MAT PRINT v4

!v5=v3/v4
DIM v5(N+1)
CALL DiffDiv(v3,v4, v5)
MAT PRINT v5

!確認
PRINT SIN(3),1.5*COS(3),1.5*COS(3)-SIN(3),-SIN(3)

END

自動微分による微分値の計算	返事を書くノートメニュー
山中和義 <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 !-------------------- ここまでサブルーチン

Re: 自動微分による微分値の計算	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/11/25 19:44:51 ** この記事は1回修正されてます
続き !●問題 !　関数 f(x1,x2,x3)=sin(1.5(x1+x2))/(x2x3) のとき、x1=x2=x3=1における関数値、偏微分係数は？ !●数値微分の場合（前進差分） DEF f(x1,x2,x3)=SIN(1.5(x1+x2))/(x2x3) !関数を定義する LET N=3 !変数の数 LET xx1=1 LET xx2=1 LET xx3=1 LET d=1E-8 !⊿ PRINT f(xx1,xx2,xx3), !f PRINT (f(xx1+d,xx2,xx3)-f(xx1,xx2,xx3))/d, !∂f/∂x1 PRINT (f(xx1,xx2+d,xx3)-f(xx1,xx2,xx3))/d, !∂f/∂x2 PRINT (f(xx1,xx2,xx3+d)-f(xx1,xx2,xx3))/d !∂f/∂x3 PRINT !●自動微分の場合 !x1=[1,1,0,0]　※[x1, ∂x1/∂x1, ∂x1/∂x2, ∂x1/∂x3] !x2=[1,0,1,0] !x3=[1,0,0,1] DIM x1(N+1),x2(N+1),x3(N+1) !変数の定義 DATA 1,1,0,0 MAT READ x1 DATA 1,0,1,0 MAT READ x2 DATA 1,0,0,1 MAT READ x3 !関数値、偏微分値がすぐにわかるような関数同士の四則演算、合成の形に分解する。 !　v1=x1+x2 !　v2=1.5v1 !　v3=sin(v2) !　v4=x2x3 !　v5=v3/v4 ! !※式を手コンパイルで求める。 !多言語（C++など）では、オペランドのオーバーロード（演算子の多重定義）機能で実現する。 !v1=x1+x2=[1,1,0,0]+[1,0,1,0]=[2,1,1,0] DIM v1(N+1) CALL DiffAdd(x1,x2, v1) MAT PRINT v1 !v2=1.5v1=[1.5,0,0,0][2,1,1,0]=[3,1.5,1.5,0] DIM c15(N+1) !定数1.5 CALL DiffSetConst(1.5, c15) DIM v2(N+1) CALL DiffMul(c15,v1, v2) !!!MAT v2=1.5v1でも可 MAT PRINT v2 !v3=sin(v2)=[sin(v2),(∂sin(v2)/∂v2)(∂v2/∂x1),(∂sin(v2)/∂v2)(∂v2/∂x2),(∂sin(v2)/∂v2)(∂v2/∂x3)]=[sin(3),1.5cos(3),1.5cos(3),0] DIM v3(N+1) !LET v3(1)=SIN(FxVal(v2)) !LET dfx=DFxVal("SIN",v2) !※∂sin(v2)/∂v2 !LET v3(2)=dfx1.5 !※∂v2/∂x1=∂{1.5(x1+x2)}/∂x1=1.5 !LET v3(3)=dfx1.5 !LET v3(4)=dfx0 CALL DiffCompo("SIN",v2, v3) !v4=x2x3=[1,0,1,0][1,0,0,1]=[1,0,1,1] DIM v4(N+1) CALL DiffMul(x2,x3, v4) MAT PRINT v4 !v5=v3/v4 DIM v5(N+1) CALL DiffDiv(v3,v4, v5) MAT PRINT v5 !確認 PRINT SIN(3),1.5COS(3),1.5COS(3)-SIN(3),-SIN(3) END