数値微分式
数値微分式 しばっち 2008/03/19 21:46:52
├どのように使うのですか 島村1243 2008/03/20 10:01:11
│└返事 しばっち 2008/03/20 13:22:28
│ └有難うございます。 島村1243 2008/03/21 22:18:20
│ └!Richardsonの補外法(extrapolation) 山中和義 2008/03/22 13:40:09
│ └こんなに短いコードで 島村1243 2008/03/22 17:14:23
├可変精度数値微分 しばっち 2008/03/20 21:07:22 (修正1回)
│├F=X^2を試したら 島村1243 2008/03/22 10:23:28 (修正1回)
││└ご指摘ありがとうございます。こちらのミス... しばっち 2008/03/22 14:59:11
││ └完璧でした。 島村1243 2008/03/22 16:36:45
│└可変精度高階数値微分 しばっち 2008/03/23 09:54:12
└高階数値微分式 しばっち 2008/03/22 17:21:13
数値微分式 しばっち 2008/03/19 21:46:52 
ツリーへ
|
数値微分式
|
 返事を書く  ノートメニュー
|
|
しばっち <dihjvcfsyu> 2008/03/19 21:46:52
|
LET MAXLEVEL=19
DIM L(MAXLEVEL),M(MAXLEVEL)
FOR N=3 TO MAXLEVEL STEP 2
FOR I=1 TO N
LET H=1
FOR J=1 TO N
IF I<>J THEN LET H=H*(I-J)
NEXT J
LET L(I)=H
NEXT I
!'FOR P=1 TO N
LET P=INT(N/2+1) !'中央差分式
FOR K=1 TO N
LET S=0
FOR I=1 TO N
LET H=1
FOR J=1 TO N
IF K<>J AND I<>J AND K<>I THEN
LET H=H*(P-J)
LET FL=1
END IF
NEXT J
IF FL=1 THEN
LET S=S+H
LET FL=0
END IF
NEXT I
LET M(K)=S
NEXT K
PRINT "F'(X)=";
FOR I=1 TO N
LET GM=GCD(M(I),L(I))
LET M(I)=M(I)/GM
LET L(I)=L(I)/GM
NEXT I
LET LM=L(1)
FOR I=2 TO N
LET LM=LCM(LM,L(I))
NEXT I
FOR I=1 TO N
LET B=LM/L(I)
LET M(I)=M(I)*B
LET L(I)=L(I)*B
NEXT I
PRINT "(";
FOR I=1 TO N
IF M(I)*L(I)<0 THEN PRINT "-"; ELSE PRINT "+";
PRINT STR$(ABS(M(I)));"*F(X";
IF (P-I)<0 THEN PRINT "+"; ELSE PRINT "-";
PRINT STR$(ABS(P-I));"*H)";
NEXT I
PRINT ")/(";STR$(ABS(L(1)));"*H)"
!'NEXT P
NEXT N
END
EXTERNAL FUNCTION LCM(M,N)
LET LCM=M*N/GCD(M,N)
END FUNCTION
EXTERNAL FUNCTION GCD(M,N)
DO WHILE N<>0
LET T=MOD(M,N)
LET M=N
LET N=T
LOOP
LET GCD=M
END FUNCTION
|
├どのように使うのですか 島村1243 2008/03/20 10:01:11 ツリーへ
|
Re: 数値微分式
|
返事を書く ノートメニュー
|
|
島村1243 <bjllmpcujp> 2008/03/20 10:01:11
|
どのように使うのですか
いつも興味を持って見させて頂いておりますが、どの様に
使うのかが分からず残念です。最初に目的や使い方の簡単
なコメントを記載して頂けると有難いです。
|
│└返事 しばっち 2008/03/20 13:22:28 ツリーへ
|
Re: どのように使うのですか
|
返事を書く ノートメニュー
|
|
しばっち <dihjvcfsyu> 2008/03/20 13:22:28
|
返事
単に数値微分式を求めさせているだけで、特に使い道は考えておりませんが
計算された式の1つをテキストウィンドウからコピペし、FUNCTIONで定義します。
使う式によって精度が違います。
FUNCTION DIFF(X)
LET DIFF=(-1*F(X-3*H)+9*F(X-2*H)-45*F(X-1*H)+0*F(X-0*H)+45*F(X+1*H)-9*F(X+2*H)+1*F(X+3*H))/(60*H)
END FUNCTION
微分値を求めたい関数も定義します。
FUNCTION F(X)
LET F=X^X-N
END FUNCTION
変数Hには1/100〜1/10000程度の数値を代入しておきます。
LET H=1/100
この計算された微分式はニュートン法等で使うことができます。
変数EPSは求める精度です。
LET EPS=1E-12
INPUT PROMPT "N=": N
LET X=N
DO
LET X=X-F(X)/DIFF(X)
LOOP UNTIL ABS(X^X-N)<EPS
PRINT X,X^X
END
|
│ └有難うございます。 島村1243 2008/03/21 22:18:20 ツリーへ
|
Re: 返事
|
返事を書く ノートメニュー
|
|
島村1243 <bjllmpcujp> 2008/03/21 22:18:20
|
有難うございます。
しばっちさん、大変良く分かりました。
微係数の求め方は中心進差分近似だけしか知らなかった
ので、しばっちさんのプログラムを動かした時にプリン
トされる沢山の式の意味が分かりませんでした。
説明を書いて頂いて、ようやくその深い意味が分かりま
した。有難うございました。
|
│ └!Richardsonの補外法(extrapolation) 山中和義 2008/03/22 13:40:09 ツリーへ
|
Re: 有難うございます。
|
返事を書く ノートメニュー
|
|
山中和義 <drdlxujciw> 2008/03/22 13:40:09
|
!Richardsonの補外法(extrapolation)
DEF f(x)=x*EXP(x) !関数の定義
FUNCTION N(j,h) !O(h^(2*j))の近似
IF j=1 THEN
LET N=(f(x0+h)-f(x0-h))/(2*h)
ELSE
LET N=N(j-1,h/2) + (N(j-1,h/2)-N(j-1,h))/(4^(j-1)-1) !j=2,3,…
END IF
END FUNCTION
LET x0=2
LET h=0.2
PRINT N(1,h) !N1(h)近似
PRINT N(2,h) !N2(h)近似
PRINT N(3,h) !N3(h)近似
PRINT N(4,h) !N4(h)近似
PRINT N(5,h) !N5(h)近似
PRINT
PRINT x0*EXP(x0)+EXP(x0) !検算
END
|
│ └こんなに短いコードで 島村1243 2008/03/22 17:14:23 ツリーへ
├可変精度数値微分 しばっち 2008/03/20 21:07:22 (修正1回) ツリーへ
|
Re: 数値微分式
|
返事を書く ノートメニュー
|
|
しばっち <dihjvcfsyu> 2008/03/20 21:07:22 ** この記事は1回修正されてます
|
可変精度数値微分
INPUT PROMPT "X=":X
LET H=1/1024
LET N=5 !'次数
LET M=(N+1)/2 !'中央差分
PRINT "微分値";DIFF(X,H,M,N)
PRINT "関数値";F(X)
END
EXTERNAL FUNCTION DIFF(XX,H,M,N)
DIM X(N),Y(N),A(N)
FOR I=1-M TO N-M
LET X(I+M)=XX+I*H
LET Y(I+M)=F(X(I+M))
NEXT I
FOR I=1 TO N
LET L=1
LET KK=1
FOR J=1 TO N
IF J<>I THEN
LET L=L*(X(I)-X(J))
LET A(KK)=X(J)
LET KK=KK+1
END IF
NEXT J
LET S1=0
FOR J=1 TO N-1
LET S=1
FOR K=1 TO N-1
IF K<>J THEN LET S=S*(XX-A(K))
NEXT K
LET S1=S1+S
NEXT J
LET SS=SS+S1*Y(I)/L
NEXT I
LET DIFF=SS
END FUNCTION
EXTERNAL FUNCTION F(X)
LET F=X^2
END FUNCTION
|
│├F=X^2を試したら 島村1243 2008/03/22 10:23:28 (修正1回) ツリーへ
|
Re: 可変精度数値微分
|
返事を書く ノートメニュー
|
|
島村1243 <bjllmpcujp> 2008/03/22 10:23:28 ** この記事は1回修正されてます
|
F=X^2を試したら
早速、可変精度数値微分のプログラムコードをBASICにコピーし
f(x)=x^2のx=3における微係数を求めてみようとして
EXTERNAL FUNCTION F(X)
LET F=X^2
END FUNCTION
と直してXの入力ダイアログに3を入力したら
X=3
微分値=9
関数値=9
と出力されてしまいました。
プログラムの使い方を誤解しているのだと思うのですが、何処が誤っているのか見当が付きません。
ご指摘頂けると有り難いです。
|
││└ご指摘ありがとうございます。こちらのミス... しばっち 2008/03/22 14:59:11 ツリーへ
|
Re: F=X^2を試したら
|
返事を書く ノートメニュー
|
|
しばっち <dihjvcfsyu> 2008/03/22 14:59:11
|
ご指摘ありがとうございます。こちらのミスでした。すみません。
早速検証し、修正しましたので再度コピペして使ってみて下さい。
多分これで大丈夫かと思います。何かあればご報告お願いします。
|
││ └完璧でした。 島村1243 2008/03/22 16:36:45 ツリーへ
│└可変精度高階数値微分 しばっち 2008/03/23 09:54:12 ツリーへ
|
Re: 可変精度数値微分
|
返事を書く ノートメニュー
|
|
しばっち <dihjvcfsyu> 2008/03/23 09:54:12
|
可変精度高階数値微分
LET H=1/512
INPUT PROMPT "X=":X
INPUT PROMPT "微分階数=":M
LET N=7 !'次数
LET NN=(N+1)/2 !'中央差分
!'LET NN=1 !'前進差分
!'LET NN=N !'後進差分
PRINT DIFF(X,H,NN,N,M);F(X)
END
EXTERNAL FUNCTION DIFF(XX,H,MM,NN,P)
DIM X(NN),Y(NN),A(NN),TMP(NN)
FOR I=1-MM TO NN-MM
LET X(I+MM)=XX+I*H
LET Y(I+MM)=F(X(I+MM))
NEXT I
FOR I=1 TO NN
LET KK=1
LET LL=1
FOR J=1 TO NN
IF J<>I THEN
LET LL=LL*(X(I)-X(J))
LET A(KK)=X(J)
LET KK=KK+1
END IF
NEXT J
LET SS=SS+FAC(P)*DIFFCALC(NN-1,XX,A,TMP,1,NN-1-P)*Y(I)/LL
NEXT I
LET DIFF=SS
END FUNCTION
EXTERNAL FUNCTION F(X)
LET F=X^10
END FUNCTION
EXTERNAL FUNCTION FAC(X)
LET S=1
FOR I=2 TO X
LET S=S*I
NEXT I
LET FAC=S
END FUNCTION
EXTERNAL FUNCTION DIFFCALC(N,XX,X(),T(),K,R)
IF R=0 THEN
LET S=1
FOR I=1 TO N
IF T(I)=1 THEN LET S=S*(XX-X(I))
NEXT I
LET DIFFCALC=S
ELSE
FOR I=K TO N-R+1
LET T(I)=1
LET SS=SS+DIFFCALC(N,XX,X,T,I+1,R-1)
LET T(I)=0
NEXT I
LET DIFFCALC=SS
END IF
END FUNCTION
|
└高階数値微分式 しばっち 2008/03/22 17:21:13 ツリーへ
|
Re: 数値微分式
|
返事を書く ノートメニュー
|
|
しばっち <dihjvcfsyu> 2008/03/22 17:21:13
|
高階数値微分式
OPTION BASE 0
PUBLIC NUMERIC N
LET MAXLEVEL=15
DIM X(1),Y(MAXLEVEL),L(MAXLEVEL),M(MAXLEVEL)
INPUT PROMPT "微分階数=":NN
IF MOD(NN,2)=0 THEN LET N1=1 ELSE LET N1=2
FOR N=NN+N1 TO MAXLEVEL STEP 2
LET P=INT(N/2+1) !'中央差分式
!' FOR P=1 TO N
FOR K=1 TO N
FOR I=1 TO N
LET Y(I)=0
NEXT I
LET Y(0)=1
FOR I=1 TO N
IF K<>I THEN
LET X(0)=-I
LET X(1)=1
CALL MUL(Y,X)
END IF
NEXT I
LET L(K)=HORNER(Y,K)
FOR I=1 TO NN
CALL DERIVATIVE(Y)
NEXT I
LET M(K)=HORNER(Y,P)
NEXT K
PRINT "d^";STR$(NN);"/dx^";STR$(NN);"F(X)=";
FOR I=1 TO N
LET GM=GCD(M(I),L(I))
LET M(I)=M(I)/GM
LET L(I)=L(I)/GM
NEXT I
LET LM=L(1)
FOR I=2 TO N
LET LM=LCM(LM,L(I))
NEXT I
FOR I=1 TO N
LET B=LM/L(I)
LET M(I)=M(I)*B
LET L(I)=L(I)*B
NEXT I
PRINT "(";
FOR I=1 TO N
IF M(I)*L(I)<0 THEN PRINT "-"; ELSE PRINT "+";
PRINT STR$(ABS(M(I)));"*F(X";
IF (P-I)<0 THEN PRINT "+"; ELSE PRINT "-";
PRINT STR$(ABS(P-I));"*H)";
NEXT I
PRINT ")/(";STR$(ABS(L(1)));"*H^";STR$(NN);")"
!'NEXT P
!'PRINT
NEXT N
END
EXTERNAL SUB MUL(A(),B())
OPTION BASE 0
DIM C(N)
FOR I=0 TO 1
FOR J=0 TO N-I
LET C(I+J)=C(I+J)+A(J)*B(I)
NEXT J
NEXT I
FOR I=0 TO N
LET A(I)=C(I)
NEXT I
END SUB
EXTERNAL SUB DERIVATIVE(A())
OPTION BASE 0
DIM B(N)
FOR I=N TO 1 STEP -1
LET B(I-1)=I*A(I)
NEXT I
FOR I=0 TO N
LET A(I)=B(I)
NEXT I
END SUB
EXTERNAL FUNCTION HORNER(A(),XX)
LET Y=A(N)
FOR I=N-1 TO 0 STEP -1
LET Y=Y*XX+A(I)
NEXT I
LET HORNER=Y
END FUNCTION
EXTERNAL FUNCTION GCD(M,N)
DO WHILE N<>0
LET T=MOD(M,N)
LET M=N
LET N=T
LOOP
LET GCD=M
END FUNCTION
EXTERNAL FUNCTION LCM (A, B)
LET LCM=A*B/GCD(A,B)
END FUNCTION
|
インデックスへ
EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。
新しく発言する