数値積分

 投稿者:島村1243  投稿日:2010年 8月22日(日)11時57分32秒
  シンプソンの公式で数値定積分を行ったのですが

微小幅を作る分割数を増加すると、積分値が大幅に増加して行き

収束値らしき値が得られず、何処に誤りが有るのか見当が付きません。

お分かりになりましたらご教示お願い致します。



LET N=8

LET D=1.432

LET L= 5e-3

LET wd=0.3e-3

LET ws=L/(N-1)


DEF y1(x,h)=SQR(h^2+(D*SIN(x/2))^2)

DEF y2(x,h)=SQR((ws-h)^2+(D*SIN(x/2))^2)


LET ndiv=500 !計算分割数

LET dx=2*PI/ndiv


!     /2π

!p(h)=| {1/y1(x,h)-1/y2(x,h)}*dx

!     /0


LET V1=p(wd/2)

LET V2=p(ws-wd/2) !理屈上はV1=-V2になる。

LET Capa=1/(V1-V2)

PRINT "V1=";V1,"V2=";V2

PRINT "計算分割数=";ndiv,"Capa=";Capa



FUNCTION p(h)  !シンプソン公式で積分

   LET sump=0

   FOR n=0 TO ndiv STEP 2

      LET x0=n*dx

      LET x1=(n+1)*dx

      LET x2=(n+2)*dx

      LET fy1=1/y1(x0,h)+4*1/y1(x1,h)+1/y1(x2,h)

      LET fy2=1/y2(x0,h)+4*1/y2(x1,h)+1/y2(x2,h)

      LET sump=sump+fy1-fy2

   NEXT N

   LET p=sump*dx/3

END FUNCTION


END
 

Re: 数値積分

 投稿者:山中和義  投稿日:2010年 8月22日(日)13時11分14秒
  > No.1364[元記事へ]

島村1243さんへのお返事です。

> シンプソンの公式で数値定積分を行ったのですが

V1= 3.70090979805763    V2=-3.70090979805763
Capa= .135101914740645

でよろしいですか? 10万分割で精度は10桁程度です。
LET N=8

PUBLIC NUMERIC D,ws
LET D=1.432
LET L=5e-3
LET wd=0.3e-3
LET ws=L/(N-1)

!     /2π
!p(h)=| {1/y1(x,h)-1/y2(x,h)}*dx
!     /0

LET ndiv=100000 !計算分割数

PUBLIC NUMERIC h
LET h=wd/2
LET V1=Simpson(0,2*PI,ndiv)
LET h=ws-wd/2
LET V2=Simpson(0,2*PI,ndiv) !理屈上はV1=-V2になる。


LET Capa=1/(V1-V2)

PRINT "V1=";V1,"V2=";V2
PRINT "Capa=";Capa

END


EXTERNAL FUNCTION f(x)
LET y1=SQR(h^2+(D*SIN(x/2))^2)
LET y2=SQR((ws-h)^2+(D*SIN(x/2))^2)
LET f=1/y1-1/y2
END FUNCTION


EXTERNAL FUNCTION Simpson(a,b,N) !シンプソン法 1/3則
LET HH=(b-a)/(2*N)
LET S1=0
FOR i=1 TO N !h/3*{f(a) + 4*Σ[i=1,n]f(a+(2*i-1)*h) + 2*Σ[i=1,n-1]f(a+(2*i)*h) + f(b)}
   LET x=a+HH*(2*i-1)
   LET S1=S1+f(x)
NEXT i
LET S2=0
FOR i=1 TO N-1
   LET x=a+HH*(2*i)
   LET S2=S2+f(x)
NEXT i
LET Simpson=HH*(f(a)+f(b)+4*S1+2*S2)/3
END FUNCTION
 

Re: 数値積分

 投稿者:島村1243  投稿日:2010年 8月23日(月)07時22分11秒
  > No.1365[元記事へ]

山中和義さんへのお返事です。

> V1= 3.70090979805763    V2=-3.70090979805763
> Capa= .135101914740645
>
> でよろしいですか? 10万分割で精度は10桁程度です。

山中さん、早速のご教示有難う御座います。
私のプログラムの何処にロジック誤りが有るのかを、山中さんのプログラムと比較して
調べていたので返信が遅くなりました。
私のプログラムで分割数を約600万にしたら、山中さんの計算結果に近づいて来ました。
その結果「ロジックは正しいが、速度と計算機誤差の累積対策が無考慮なので使用に耐え
ない」と判断しました。

山中さんのコードを採用し、結果をSECONDさんの浮遊容量計算に適用したら、
 コイル端子間合成浮遊容量= 123.965055354134 pF
となりました。実機データに近づいたのか否かは分かりません。

PUBLIC NUMERIC D,ws
LET N=8
LET D=1.432
LET L=5e-3
LET wd=0.3e-3
LET ws=L/(N-1)
LET ndiv=100000 !計算分割数

PUBLIC NUMERIC h
LET c0=2.99792458e8                     !m/s      光速度
LET e0=1e7/(4*PI*c0^2)                  !F/m 真空の誘電率
LET h=wd/2
LET V1=Simpson(0,2*PI,ndiv)
LET h=ws-wd/2
LET V2=Simpson(0,2*PI,ndiv) !理屈上はV1=-V2になる。
LET capa=8*PI^2*e0/(V1-V2)

PRINT "コイル端子間合成浮遊容量=";css1(N)*1e12;"pF"
!PRINT "V1=";V1,"V2=";V2

FUNCTION css1(N)
   LET sumc=0
   FOR k=1 TO N-1
      LET c=(N-k)*capa/k/(N/k)^2
      LET sumc=sumc+c
   NEXT k
   LET css1=sumc
END FUNCTION

END

EXTERNAL FUNCTION f(x)
  LET y1=SQR(h^2+(D*SIN(x/2))^2)
  LET y2=SQR((ws-h)^2+(D*SIN(x/2))^2)
  LET f=1/y1-1/y2
END FUNCTION

EXTERNAL FUNCTION Simpson(a,b,N) !シンプソン法 1/3則
  LET HH=(b-a)/(2*N)
  LET S1=0
  FOR i=1 TO N
     LET x=a+HH*(2*i-1)
     LET S1=S1+f(x)
  NEXT i
  LET S2=0
  FOR i=1 TO N-1
     LET x=a+HH*(2*i)
     LET S2=S2+f(x)
  NEXT i
  LET Simpson=HH*(f(a)+f(b)+4*S1+2*S2)/3
END FUNCTION
 

Re: 数値積分

 投稿者:SECOND  投稿日:2010年 8月23日(月)18時20分42秒
  > No.1366[元記事へ]

島村1243さんへのお返事です。

プログラムで検算、して頂いたようで、ありがとうございます。

元々は、ゲルマ・ラジオの、ループ・アンテナ・コイルで、
その両端子に換算される 浮遊容量が、分ると
バリコン無しで、特定ラジオ局に同調するように巻ける・・程度の実験と、
方法ですので、誤差の10数%くらいは、区別できないです。

この方面の、論文を、見つけましたが、私には少々難解です。
長岡係数完全楕円積分総集編 24~25ページ付近の項、4.分布容量
 

Re: 数値積分

 投稿者:島村1243  投稿日:2010年 8月23日(月)20時22分3秒
  > No.1367[元記事へ]

SECONDさんへのお返事です。

> この方面の、論文を、見つけました。
> 長岡係数完全楕円積分総集編 24~25ページ付近の項、4.分布容量

(実行プログラム有りませんがご容赦願います。)
さて、貴重な論文「4.分布容量」見させて頂きました。
私が今回試算した式は、上の論文に出ている考え方と全く同じ考え方で導出
している(ことが確認できた)ので、論文の内容が弱電の世界で一般的に
通用しているものならば、実用に耐えると思います。
 

Re: 数値積分

 投稿者:島村1243  投稿日:2010年 8月23日(月)21時05分6秒
  > No.1365[元記事へ]

山中和義さんへのお返事です。

> > シンプソンの公式で数値定積分を行ったのですが
>
> V1= 3.70090979805763    V2=-3.70090979805763
> Capa= .135101914740645
>
> でよろしいですか? 10万分割で精度は10桁程度です。

度々の投稿で恐縮です。私の最初のプログラムコードは内部関数の中からDEF関数を
呼んでいることで変な挙動になる事が分かりました。その経緯を書きます。

(1)私の最初にお尋ねしたプログラムコードが何故ダメなのかを調べていたら
内部関数「FUNCTION p(h)」の中のコードにロジックミスが見つかりました。

誤り:「FOR n=0 TO ndiv STEP 2」
正 :「FOR n=0 TO ndiv-2 STEP 2」

(2)上記を修正してRUNしましたが、結論はやはり山中さんのプログラム計算結果に
なりません。

(3)そこで試しに、FUNCTION p(h)を外部関数に変更し、DEF関数を外部関数内で定義
する方法に変更して見ました。末尾に変更後のコードを示します。すると、山中さんの
プログラムと同じ(最後の桁で少し異なるが)値が出ました。

(4)この結果から、内部関数の中で、外に有るDEF関数を使用したことにトラブルの
原因があると分かりました。

お尋ねしたいことは、
1)内部関数中で、メインに有るDEF関数を使用すると一般的に挙動が変になるのか?
2)一般的に使用できるが、私のコードがダメなので挙動が変になったのか?
のいずれなのだろうか、ということです。

なお、内部関数の中でDEF関数を定義したらエラーになって受け付けてくれませんでした。

----<変更後のコード>----
PUBLIC NUMERIC D,ws,wd,ndiv
LET N=8
LET D=1.432
LET L= 5e-3
LET wd=0.3e-3
LET ws=L/(N-1)
LET ndiv=100000 !計算分割数

LET V1=p(wd/2)
LET V2=p(ws-wd/2)
LET Capa=1/(V1-V2)
PRINT "V1=";V1,"V2=";V2
PRINT "計算分割数=";ndiv,"Capa=";Capa
END

EXTERNAL FUNCTION p(h)  !シンプソン公式で積分
DEF y1(x)=SQR(h^2+(D*SIN(x/2))^2)
DEF y2(x)=SQR((ws-h)^2+(D*SIN(x/2))^2)
LET dx=2*PI/ndiv
LET sump=0
FOR n=0 TO ndiv-2 STEP 2
   LET x0=n*dx
   LET x1=(n+1)*dx
   LET x2=(n+2)*dx
   LET fy1=1/y1(x0)+4*1/y1(x1)+1/y1(x2)
   LET fy2=1/y2(x0)+4*1/y2(x1)+1/y2(x2)
   LET sump=sump+fy1-fy2
NEXT N
LET p=sump*dx/3
END FUNCTION
 

Re: 数値積分

 投稿者:SECOND  投稿日:2010年 8月24日(火)00時01分30秒
  > No.1366[元記事へ]

島村1243さんへのお返事です。

面白いことが分りました。

ws=L/(N-1)で計算されているようですが、私の計算では、ws=L/N で計算しています。
この ws を、島村さんの式へ適用すると、なんと同じ値で、出てきます。
微小区間の積分な為か、あまり曲率の影響が無いという結果が出たのでしょう。

142.489293249106 pF と コイル端子間合成浮遊容量= 142.489267952594 pF
巻径が大き過ぎるせいかとも思いますが、14mm でも、1/1000 もないようです。
D=1.432e-1   14.2489293249106 pF   コイル端子間合成浮遊容量= 14.2487466579388 pF
D=1.432e-2   1.42489293249106 pF   コイル端子間合成浮遊容量= 1.42377272689129 pF
 

Re: 数値積分

 投稿者:山中和義  投稿日:2010年 8月24日(火)08時21分15秒
  > No.1369[元記事へ]

島村1243さんへのお返事です。

> 1)内部関数中で、メインに有るDEF関数を使用すると一般的に挙動が変になるのか?

前プログラムと改修プログラムの計算結果は同じになります。
DEF文や内部関数との関連性は、特に問題はないと思います。


> 2)一般的に使用できるが、私のコードがダメなので挙動が変になったのか?

小生のプログラムとの違いは、以下が理由だと思います。

精度の保証

・分割数は実際は、2倍で計算している。
 分割数に奇数を指定された場合を考慮しています。したがって、10万分割ではなく、20万分割です。

 LET dx=2*PI/ndiv

 LET HH=(b-a)/(2*N)


・桁落ちが発生する
 ほぼ等しい数どうしの差では、有効桁が半減する。 Σ(a-b)とΣa-Σbとの違い。
 EXTERNAL FUNCTION p(h)  !シンプソン公式で積分
 DEF y1(x)=SQR(h^2+(D*SIN(x/2))^2)
 DEF y2(x)=SQR((ws-h)^2+(D*SIN(x/2))^2)
 LET dx=2*PI/ndiv
 LET sump=0
 FOR n=0 TO ndiv-2 STEP 2
    LET x0=n*dx
    LET x1=(n+1)*dx
    LET x2=(n+2)*dx
    LET fy1=1/y1(x0)+4*1/y1(x1)+1/y1(x2)
    LET fy2=1/y2(x0)+4*1/y2(x1)+1/y2(x2)
    LET sump=sump+fy1-fy2 !※←←←←←
 NEXT N
 LET p=sump*dx/3
 END FUNCTION
 EXTERNAL FUNCTION p(h)  !シンプソン公式で積分
 DEF y1(x)=SQR(h^2+(D*SIN(x/2))^2)
 DEF y2(x)=SQR((ws-h)^2+(D*SIN(x/2))^2)
 LET dx=2*PI/ndiv
 LET s1=0
 LET s2=0
 FOR n=0 TO ndiv-2 STEP 2
    LET x0=n*dx
    LET x1=(n+1)*dx
    LET x2=(n+2)*dx
    LET s1=s1 + 1/y1(x0)+4*1/y1(x1)+1/y1(x2)
    LET s2=s2 + 1/y2(x0)+4*1/y2(x1)+1/y2(x2)
 NEXT N
 LET p=(s1-s2)*dx/3
 END FUNCTION
 

Re: 数値積分

 投稿者:島村1243  投稿日:2010年 8月24日(火)09時16分28秒
  > No.1372[元記事へ]

山中和義さんへのお返事です。

> 前プログラムと改修プログラムの計算結果は同じになります。
> DEF文や内部関数との関連性は、特に問題はないと思います。
> 小生のプログラムとの違いは、以下が理由だと思います。
>
> 精度の保証
>
> ・分割数は実際は、2倍で計算している。
>  分割数に奇数を指定された場合を考慮しています。したがって、10万分割ではなく、20万分割です。
> ・桁落ちが発生する
>  ほぼ等しい数どうしの差では、有効桁が半減する。 Σ(a-b)とΣa-Σbとの違い。

内部関数とDEF文との関係も問題ない(当方の操作ミスだった様です)こと確認しました。
精度についても詳細な解説を頂き良く分かりました。
どうも有難う御座いました。
 

Re: 数値積分

 投稿者:島村1243  投稿日:2010年 8月24日(火)21時55分44秒
  > No.1371[元記事へ]

SECONDさんへのお返事です。

> 面白いことが分りました。
>
> ws=L/(N-1)で計算されているようですが、私の計算では、ws=L/N で計算しています。
> この ws を、島村さんの式へ適用すると、なんと同じ値で、出てきます。
> 微小区間の積分な為か、あまり曲率の影響が無いという結果が出たのでしょう。
pF


8巻のコイルは、計算上リング状導体が8個と想定するので、リング間隔は(8-1)個になると考えws=L/(N-1)としました。

さて今までの合成計算では、最初のリングと直ぐ隣(1番目)のリング間静電容量をCoとすると、
 最初のリングと2番目隣リング間静電容量は、Coが2個直列なのでCo/2
 最初のリングとk番目隣リング間静電容量は、Coがk個直列なのでCo/k
としています。

しかし、よくよく考えると
 最初のリングと2番目リング間は距離が2倍になるが静電容量は単純に1/2にならない。
と思いました。
そこでリング間の距離を変えて各グループの静電容量を全て計算し、合成したらどうなるかを試算して見ました。(プログラムは末尾)

結果は
k= 1 ,capk= 94.4500797302147 pF <--- Co
k= 2 ,capk= 58.3980276605647 pF <---単純に1/2となっていない。
k= 3 ,capk= 48.3779952764142 pF <---単純に1/3となっていない。
k= 4 ,capk= 43.255485010341 pF
k= 5 ,capk= 40.0164539074001 pF
k= 6 ,capk= 37.7288933988779 pF
k= 7 ,capk= 35.9996371056512 pF
226.40250993776 pF
と大幅に増えましたが、もしかしたら、こちらの値の方が実機に近いのではないでしょうか。

PUBLIC NUMERIC D,ws,wd,e0,ndiv
LET ndiv=100000 !計算分割数
LET N=8       !回
LET D=1.432   !m
LET L= 5e-3   !m
LET wd=0.3e-3 !m
LET ws=L/(N-1) !m

LET c0=2.99792458e8                     !m/s      光速度
LET e0=1e7/(4*PI*c0^2)                  !F/m 真空の誘電率

PRINT css1(N)*1e12;"pF"

FUNCTION css1(N)
   LET sumc=0
   FOR k=1 TO N-1
      LET h=wd/2       !最初の導体表面位置指定
      LET V1=p(h,k)    !最初の導体表面電位
      LET h=k*ws-wd/2  !k番目隣の導体表面位置指定
      LET V2=p(h,k)    !k番目隣の導体表面電位
      LET capk=1/(V1-V2)  !最初の導体とk番目隣の導体間の静電容量
      PRINT "k=";k;",capk=";capk*1e12;"pF"
      LET c=(N-k)*capk/(N/k)^2  !端子間に換算
      LET sumc=sumc+c
   NEXT k
   LET css1=sumc
END FUNCTION

!-----------------------------------------
!    ループの数 1ループ当りの容量  両端子に写る倍率
! c1 =  (N-1)      *capk[k=1の値]      /(N/1)^2
! c2 =  (N-2)      *capk[k=2の値]      /(N/2)^2
! c3 =  (N-3)      *capk[k=3の値]      /(N/3)^2
!  ・
!  ・
! cn-1=(N-(N-1))   *capk[k=N-1の値]   /(N/(N-1))^2

END

EXTERNAL FUNCTION p(h,k)
DEF y1(x)=SQR(h^2+(D*SIN(x/2))^2)
DEF y2(x)=SQR((k*ws-h)^2+(D*SIN(x/2))^2)
LET dx=2*PI/ndiv
LET sump1=0
LET sump2=0
FOR n=0 TO ndiv-2 STEP 2
   LET x0=n*dx
   LET x1=(n+1)*dx
   LET x2=(n+2)*dx
   LET fy1=1/y1(x0)+4*1/y1(x1)+1/y1(x2)
   LET fy2=1/y2(x0)+4*1/y2(x1)+1/y2(x2)
   LET sump1=sump1+fy1
   LET sump2=sump2+fy2
NEXT N
LET p=(sump1-sump2)*dx/3/(8*PI^2*e0)
END FUNCTION
 
 

Re: 数値積分

 投稿者:SECOND  投稿日:2010年 8月24日(火)23時03分13秒
  > No.1375[元記事へ]

島村1243さんへのお返事です。

あとは、正確なコイルを、実際に作成され、実測で確かめる事をしないと進めませんので、
ご自身で追求される事を期待します。結果がでましたら、教えてください。
この件については、これで打ち切らせてください。ありがとうございました。


    L    コイルは、ヘリカル構造なため、巻き始め終りの端子側と、その反対側で
  ||    Lの長さが 異なります。どちら側をLとするかで 以下。 (左図は N=2 )
  ○○
               (端子の無い側のL)(端子の有る側のL)  同じコイルのWSの
 ○○○   WS=  ─────────=─────────    値は、同値であること
 | |           (N-1)        N
  L
 

Re: 数値積分

 投稿者:みなと  投稿日:2010年11月14日(日)21時28分57秒
  > No.1367[元記事へ]

SECONDさんへのお返事です。

私の本の記事を読んでくださり、ありがとうございます。
執筆しました、みなとと申します。
円形対の静電容量に関する私の記事が話題に上った、
一連の流れを拝見させていただきました。

私もコイルの自己共振周波数について関心があり、
コイルの巻線と巻線の間に現れる静電容量や、
巻線と地面との静電容量をどのように把握し、
自己共振周波数を求めるためにどのように用いるか、
という点でまだまだ悩んでおります。

現在つまづいている点は、巻線の位置によって
電流が異なることはどう考慮するのかということです。

http://www.geocities.jp/kouyoubako/index.htm

 

Re: 数値積分

 投稿者:SECOND  投稿日:2010年11月15日(月)21時36分2秒
  > No.1444[元記事へ]

みなと 様
 思いがけないご返事、おそれいります。あの当時、探す事のできる唯一のもの
 でした。ありがとうございました。・・私には少し難しいです。
 

Re: 数値積分

 投稿者:島村1243  投稿日:2010年11月20日(土)08時34分37秒
  > No.1444[元記事へ]

みなとさんへのお返事です。

(1.1)式の基になる方程式はd/dtを含んでいたが、d/dtをjωに置き換えて計算した結果が(1.1)式になる、と解釈(意味は理解できていません)しました。

みなとさんの書かれた

>現在つまづいている点は、巻線の位置によって
>電流が異なることはどう考慮するのか

の主旨は、「1本の導体上を通過する電流は同一値である筈なのに、(1.1)式を使って求めると位置によって異なる結果になるのは何故か?」ということでしょうか?
もしもその主旨であるなら、下記送電線での現象(進行波と定在波の違い)が参考になるかも知れません。
「電流が異なる(境界条件)は、どう(式に)考慮(反映)すれば良いのか」という主旨でしたら下記は的外れですので無視してください。

無損失長距離送電線の一端に正弦波電圧を加えると、送電線上の電圧・電流は、送電線に分布する自己インダクタンスL[H/m]と選間静電容量C[F/m]の影響によって時間tと位置zを変数とする波動方程式で表されます。
その解をt,zの関数である進行波(或いは後進波)として求めると、進行波の最大値は電線上を進行するとき変わりません。しかし、この波動方程式中のd/dtをjωに置き換えて計算すると、進行波としてではなく定在波を得ることになり、1本の電線上であるにも関わらず電流の最大値はA点からの位置に応じて異なります。
 

戻る