自動微分でヤコビ行列を求める場合

　　連立方程式あるく 2007/11/21 18:11:37 (修正2回) 　　├ガウスの消去法は，通常，連立一次方程式を... 白石　和夫 2007/11/21 21:04:24 　　│└いえ、ガウス法ではなく、ガウスの消去法の... あるく 2007/11/22 23:17:59 　　│　└非線形の連立方程式を解くためのガウス消去... 白石　和夫 2007/11/23 09:06:48 　　│　　└!２元非線形連立方程式の解山中和義 2007/11/23 13:15:07 (修正1回) 　　│　　　├!オセッカイですが・・ SECOND 2007/11/24 00:21:32 　　│　　　│└質問なのですが、皆さんが使ってる"!"はどう... あるく 2007/11/24 19:05:44 　　│　　　│　└ヘルプファイルに、詳しい説明があります。 SECOND 2007/11/25 00:05:30 (修正1回) 　　│　　　├2.312376477871337.31237647787135 あるく 2007/11/27 17:43:16 (修正1回) 　　│　　　│└そうです。山中和義 2007/11/27 18:44:39 　　│　　　│　└ありがとうございます。疑問が解消しました... あるく 2007/11/27 19:17:46 　　│　　　├自動微分でヤコビ行列を求める場合山中和義 2007/11/28 21:30:22

自動微分でヤコビ行列を求める場合

精度の違いで多少収束が早い。

!２元非線形連立方程式の解

OPTION ARITHMETIC complex

!自動微分による微分係数の計算（ボトムアップ方式）
SUB DiffSQR(g(), f()) !合成 sqr(g)
LET f(1)=SQR(g(1)) !h(g)
FOR i=2 TO UBOUND(f)
LET f(i)=g(i)/(2*SQR(g(1))) !∂h(g)/∂g * ∂g/∂xi
NEXT i
END SUB

DIM x(3),y(3)
DATA 1,1,0 ![1, ∂f/∂x, ∂f/∂y]
MAT READ x
DATA 1,0,1 ![1, ∂g/∂x, ∂g/∂y]
MAT READ y

DIM c5(3) !定数5
DATA 5,0,0
MAT READ c5

FOR k=1 TO 15 !漸化式を収束させる

DIM xy(2) !⊿x、⊿y

!連立１次方程式Jx=bを得る
DIM v1(3),v2(3),v3(3),v4(3)
MAT v1=x+c5 !f(x,y)=x+5-y
MAT v2=v1-y
CALL DiffSQR(x, v3) !g(x,y)=SQR(x)-y
MAT v4=v3-y

DIM J(2,2)
LET J(1,1)=v2(2) !J=(∂f/∂x　∂f/∂y)　ヤコビ(Jacobi)行列
LET J(1,2)=v2(3) !　(∂g/∂x　∂g/∂y)
LET J(2,1)=v4(2)
LET J(2,2)=v4(3)

DIM b(2)
LET b(1)=-v2(1) !b=(-f)
LET b(2)=-v4(1) !　(-g)

!----------　※ガウスの消去法などで２元連立１次方程式を解く
DIM Ji(2,2) !逆行列を求める　Full BASIC版
MAT Ji=INV(J) !正則なら
MAT xy=Ji*b
!----------

LET x(1)=x(1)+xy(1) !Xi+1=Xi+⊿x
LET y(1)=y(1)+xy(2) !Yi+1=Yi+⊿y
PRINT x(1),y(1) !確認　∥Xi+1 - Xi∥≦ε(1+∥Xi+1∥)

NEXT k

DEF f(x,y)=x+5-y !f(x)=x+5
DEF g(x,y)=SQR(x)-y !g(x)=SQR(x)
PRINT f(x(1),y(1)),g(x(1),y(1)) !確認　∥f(xi)∥<δ

END

　　│　　　└たびたびすみません。また質問なのですが、あるく 2007/12/02 16:39:20 　　│　　　　└たとえば, 白石　和夫 2007/12/02 16:52:02 (修正1回) 　　│　　　　　└なるほど。あるく 2007/12/02 18:08:53 (修正1回) 　　│　　　　　　└わかりやすい記述に変更してください。山中和義 2007/12/03 10:04:15 　　└!f(x)=x+5 SECOND 2007/11/22 18:47:13 (修正1回) 　　　└ありがとうございます。あるく 2007/11/22 23:25:17 　　　　└!複素数平面（ガウス平面）で根の見当をつけ... 山中和義 2007/11/23 13:52:37 (修正1回)

Re: !２元非線形連立方程式の解	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/11/28 21:30:22
自動微分でヤコビ行列を求める場合精度の違いで多少収束が早い。 !２元非線形連立方程式の解 OPTION ARITHMETIC complex !自動微分による微分係数の計算（ボトムアップ方式） SUB DiffSQR(g(), f()) !合成 sqr(g) LET f(1)=SQR(g(1)) !h(g) FOR i=2 TO UBOUND(f) LET f(i)=g(i)/(2SQR(g(1))) !∂h(g)/∂g ∂g/∂xi NEXT i END SUB DIM x(3),y(3) DATA 1,1,0 ![1, ∂f/∂x, ∂f/∂y] MAT READ x DATA 1,0,1 ![1, ∂g/∂x, ∂g/∂y] MAT READ y DIM c5(3) !定数5 DATA 5,0,0 MAT READ c5 FOR k=1 TO 15 !漸化式を収束させる DIM xy(2) !⊿x、⊿y !連立１次方程式Jx=bを得る DIM v1(3),v2(3),v3(3),v4(3) MAT v1=x+c5 !f(x,y)=x+5-y MAT v2=v1-y CALL DiffSQR(x, v3) !g(x,y)=SQR(x)-y MAT v4=v3-y DIM J(2,2) LET J(1,1)=v2(2) !J=(∂f/∂x　∂f/∂y)　ヤコビ(Jacobi)行列 LET J(1,2)=v2(3) !　(∂g/∂x　∂g/∂y) LET J(2,1)=v4(2) LET J(2,2)=v4(3) DIM b(2) LET b(1)=-v2(1) !b=(-f) LET b(2)=-v4(1) !　(-g) !----------　※ガウスの消去法などで２元連立１次方程式を解く DIM Ji(2,2) !逆行列を求める　Full BASIC版 MAT Ji=INV(J) !正則なら MAT xy=Ji*b !---------- LET x(1)=x(1)+xy(1) !Xi+1=Xi+⊿x LET y(1)=y(1)+xy(2) !Yi+1=Yi+⊿y PRINT x(1),y(1) !確認　∥Xi+1 - Xi∥≦ε(1+∥Xi+1∥) NEXT k DEF f(x,y)=x+5-y !f(x)=x+5 DEF g(x,y)=SQR(x)-y !g(x)=SQR(x) PRINT f(x(1),y(1)),g(x(1),y(1)) !確認　∥f(xi)∥<δ END