微分方程式の数値解法

 投稿者：しまむら1243  投稿日：2009年 4月10日(金)21時05分30秒
  次の微分方程式を、数値計算で解く場合の方法を教えてください。

回路の微分方程式
（一次側の電圧平衡式）L1*Δi1/Δt+R1*(1+sin(w*t))*i1-M*Δi2/Δt=E
（二次側の電圧平衡式）M*Δi1/Δt=L2*Δi2/Δt+R2*i2

変圧器で結合された電気回路の式で記号の意味は次のとおりです。
　t :時間[s]
　E :一次側直流定電圧電源[V] 例えば10[V]
　i1:一次側電流[A]
　L1:変圧器一次側自己インダクタンス[H] 例えば0.1[H]
  r :一次側直列変動抵抗で、r=R1*(1+sin(ω*t)
      R1:一定抵抗[Ω] 例えば1000[Ω]
　　　 ω :可変抵抗rの変動角周波数[rad/s]　例えば2*PI*50[rad/s]
　i2:二次側電流[A]
　L2:変圧器二次側自己インダクタンス[H] 例えば0.1[H]
  R2:変圧器二次負荷抵抗[Ω]  例えば10[Ω]
　M:変圧器一次二次相互インダクタンス[H]=SQR(L1*L2)

ルンゲ・クッタ法が使えそうも無く、単純なオイラー法（修正無）で解こうとしたが発散してうまくいきません。よろしくお願いします。 		 

Re: 微分方程式の数値解法

 投稿者：SECOND  投稿日：2009年 4月13日(月)08時56分12秒
  > No.314[元記事へ]

!こんにちは、おひさしぶりです。ご参考程度に、自信がなくて・・

!以下は、ルンゲクッタ法で、時系列の数値解と、その動きを描いてみましたが、
!クリチカルな点など、いくつか不明な面もあります。


!回路の微分方程式
!（一次側の電圧平衡式）L1*Δi1/Δt+R1*(1+sin(w*t))*i1-M*Δi2/Δt=E
!（二次側の電圧平衡式）M*Δi1/Δt=L2*Δi2/Δt+R2*i2
!---------------------------
!計算用の微分方程式
! f1(t, i1)= (di1/dt)= ( E-R1*(1+sin(w*t))*i1+M*(di2/dt) )/L1
! f2(t, i2)= (di2/dt)= ( M*(di1/dt)-R2*i2 )/L2
!
!※上のままでは、正帰還 暴走するので、f1(),f2() に、バッファ f1_,f2_ を付けた。
DEF f1(t, i1)=( E-R1*(1+SIN(w9*t))*i1+M*f2_ )/L1 ! f2_=f2(t, i2)
DEF f2(t, i2)=( M*f1_-R2*i2 )/L2                 ! f1_=f1(t, i1)

LET E=10
LET R1=1000
LET i1=E/R1/2
LET i2=0
LET L1=50 !0.1 小さいと、オーバーフロー？
LET L2=2 !0.1 小さいと、オーバーフロー？
LET R2=10
LET M=SQR(L1*L2)
LET w9=2*PI*1 !50 小さいと、描画ピッチ演算ピッチ共に無理
!
LET dt=0.05 !sec. 演算ピッチ。描画ピッチが遅れない程度に。

SUB RungeKutta4_1(t, i1)
   LET k1=f1(t, i1)
   LET k2=f1(t+dt/2, i1+k1*dt/2)
   LET k3=f1(t+dt/2, i1+k2*dt/2)
   LET k4=f1(t+dt,   i1+k3*dt )
   LET i1=i1+(k1+2*k2+2*k3+k4)*dt/6
   LET f1_=f1(t, i1)
END SUB

SUB RungeKutta4_2(t, i2)
   LET k1=f2(t, i2)
   LET k2=f2(t+dt/2, i2+k1*dt/2)
   LET k3=f2(t+dt/2, i2+k2*dt/2)
   LET k4=f2(t+dt,   i2+k3*dt )
   LET i2=i2+(k1+2*k2+2*k3+k4)*dt/6
   LET f2_=f2(t, i2)
END SUB

!----run
LET t=0
LET w=.5 !13
SET WINDOW -w,w,-w,w
SET LINE width 3
SET COLOR MIX(15) .5,.5,.5
SET TEXT background "OPAQUE"
LET t0=TIME
DO
   LET t1=TIME
   IF dt=< ABS(t1-t0) THEN
      SET DRAW mode hidden
      CLEAR
      DRAW grid(5,5)
      PLOT TEXT,AT w*0.25,w*0.9:"マウス 右ボタンで、終了。"
      PLOT TEXT,AT w*0.4,w*0.76,USING"演算ピッチ=#.### 秒":dt
      PLOT TEXT,AT w*0.4,w*0.69,USING"描画ピッチ=#.### 秒":t1-t0
      LET t0=t1
      PRINT t;i1;i2;f1_;f2_
      !---
      PLOT TEXT,AT-w/2.5,i1:"i1"
      PLOT TEXT,AT w/2.9,i2:"i2"
      PLOT LINES :-w/3,i1;0,i1
      PLOT LINES : w/3,i2;0,i2
      !---
      SET DRAW mode explicit
      CALL RungeKutta4_1(t,i1) !　次のi1 へ更新
      CALL RungeKutta4_2(t,i2) !　次のi2 へ更新
      LET t=t+dt
   END IF
   WAIT DELAY 0 !　省電力効果
   MOUSE POLL mx,my,mlb,mrb
LOOP UNTIL mrb=1

END 		 

Re: 微分方程式の数値解法

 投稿者：しまむら1243  投稿日：2009年 4月13日(月)19時13分9秒
  > No.315[元記事へ]

SECONDさんへのお返事です。

SECONDさん、こんにちわ。ご無沙汰しております。変な題材にも関わらず、早速のご教示ありがとうございます。

未だ数値計算の内容を読み終えていないのですが、早速、作成していただいたプログラムを走らせて見ました。
電流が振動している様子は、横軸に時間、縦軸に大きさを採って表すのが一般的ですが、SECONDさんの発想は斬新的だと思いました。

さて質問の意図を記載しなくて申し訳ありませんでしたが、トランジスタの教科書では、「変圧器結合のエミッタ接地A級増幅回路のコレクタ～エミッタ間電圧は、電圧E（一般にはVccと書く）を中心として±Eの範囲で振動する」と書かれています。
つまり、変圧器の一次側電圧のピーク値は2Eになることを意味しているのですが、何故E以上になりうるのか？を数値計算で確認したかったものです。

提示した微分方程式の抵抗rは、トランジスタのベース電流を正弦波に振動させたときの電流増幅機能を、変動する抵抗rで置き換えて模擬したつもりです。
（教科書では電流増幅機能で説明しているが、トランジスタが電流を増幅するのではなく、ベース電流が変化することによってトランジスタのC～E間抵抗が変化し、その結果、外部電源Eから流れる電流i1が大きく変動して電流が増幅された様に見える、と思っているので。）

そして最終的に知りたいことは、トランスの一次側電圧 「Vt=E-r*i1 」は、直流電源電圧がEしか無いにも関わらず、Eを中心として０～2Eの範囲で正弦波状に振動する、つまり確かに「 2E>=Vt>=0 」で正弦波杖に振動する事を見たいのです。

この様な意図であるため、波形図は横軸に時間t[s]、縦軸にi1[A]（最終的にはVt[V]）として頂くと有難いです。勝って言って申し訳ありません。 		 

Re: 微分方程式の数値解法

 投稿者：SECOND  投稿日：2009年 4月13日(月)20時30分23秒
  > No.316[元記事へ]

しまむら1243さんへのお返事です。

大変、申し訳有りませんが、ご期待には答えられません。それよりこのテーマにて
ルンゲクッタでの描画とその動きに、いくつか不明な部分が見つかりまして、それに
気をとられております。できれば、その追跡にお力添えを、頂きたかったのですが。
トランス結合Ａ級増幅器にスピーカがつながった状態にしては、安定性を欠き、
正帰還暴走したり、クリチカルに過ぎるのでは、ないだろうか？ 		 

Re: 微分方程式の数値解法

 投稿者：しまむら1243  投稿日：2009年 4月14日(火)04時18分58秒
  > No.317[元記事へ]

SECONDさんへのお返事です。

> しまむら1243さんへのお返事です。
>
> 大変、申し訳有りませんが、ご期待には答えられません。それよりこのテーマにて
> ルンゲクッタでの描画とその動きに、いくつか不明な部分が見つかりまして、それに
> 気をとられております。できれば、その追跡にお力添えを、頂きたかったのですが。
> トランス結合Ａ級増幅器にスピーカがつながった状態にしては、安定性を欠き、
> 正帰還暴走したり、クリチカルに過ぎるのでは、ないだろうか？

波形を描くように手を入れさせて頂きました。すると
１）変圧器一次電圧波形は正弦波Vt=Ｅsin(wt)[V]が出ました。 感激!!です。
でも私はVt=E+Esin(wt)を想定していたので正しい結果なのかが判断できず。

２）電流波形は正弦波振動しません。１）と併せて何処かにロジック上の問題があるようですが、未だ分かりません。

!------- 以下変更したプログラム------
!※上のままでは、正帰還 暴走するので、f1(),f2() に、バッファ f1_,f2_ を付けた。
DEF f1(t, i1)=( E-R1*(1+SIN(w9*t))*i1+M*f2_ )/L1 ! f2_=f2(t, i2)
DEF f2(t, i2)=( M*f1_-R2*i2 )/L2                 ! f1_=f1(t, i1)

LET E=10
LET R1=1000
LET L1=0.5 !0.2以下にするとオーバーフロー？
LET L2=0.1
LET R2=10
LET M=SQR(L1*L2)
LET frq=50 !信号周波数[Hz]
LET dt=0.001 !微分時間[s]

LET drawTime=0.1 ![s] !描画時間で、frq×drawTimeサイクル分だけ描画する。
LET iBairitu=500 !電流波形描画拡大倍率

LET i1_=0 ! i1の初期値
LET i2_=0 ! i2の初期値
LET Vt_=E !変圧器一次電圧Vtの初期値
LET t_=0  ! tの初期値

SUB RungeKutta4_1(t_, i1_,i1)
   LET k1=f1(t_, i1_)
   LET k2=f1(t_+dt/2, i1_+k1*dt/2)
   LET k3=f1(t_+dt/2, i1_+k2*dt/2)
   LET k4=f1(t_+dt,   i1_+k3*dt )
   LET i1=i1_+(k1+2*k2+2*k3+k4)*dt/6
   LET f1_=f1(t_, i1)
END SUB

SUB RungeKutta4_2(t_, i2_,i2)
   LET k1=f2(t_, i2_)
   LET k2=f2(t_+dt/2, i2_+k1*dt/2)
   LET k3=f2(t_+dt/2, i2_+k2*dt/2)
   LET k4=f2(t_+dt,   i2_+k3*dt )
   LET i2=i2_+(k1+2*k2+2*k3+k4)*dt/6
   LET f2_=f2(t_, i2)
END SUB

!----run
LET w=2*PI*frq !交流信号の角周波数[rad/s]
LET nmax=drawTime/dt !計算点数

SET WINDOW -0,drawTime,-15,15
DRAW grid

FOR n=1 TO nmax
   LET t=n*dt
   CALL RungeKutta4_1(t_,i1_,i1)
   CALL RungeKutta4_2(t_,i2_,i2)
   LET Vt=E-R1*(1+SIN(w*t))*i1

   SET LINE COLOR "red" !電流i1波形描画色
   PLOT LINES:t_,i1_*iBairitu;t,i1*iBairitu

   SET LINE COLOR "black" !電圧波形描画色
   PLOT LINES:t_,Vt_;t,Vt

   PLOT LINES: t_,0;t,0 !時間軸描画

   LET t_=t   !　次のt へ更新
   LET i1_=i1 !　次のi1 へ更新
   LET i2_=i2 !　次のi2 へ更新
   LET Vt_=Vt !　次のi2 へ更新
Next N

END 		 

Re: 微分方程式の数値解法

 投稿者：山中和義  投稿日：2009年 4月14日(火)08時28分2秒
  > No.318[元記事へ]

しまむら1243さんへのお返事です。

●このプログラムについて

　DEF f1(t, i1)=( E-R1*(1+SIN(w9*t))*i1+M*f2_ )/L1 ! f2_=f2(t, i2)
と
　LET Vt=E-R1*(1+SIN(w*t))*i1

で、ωの変数名が違う。

改修案
　DEF文のw9をwとして、
　刻み幅を、LET dt=0.000001 !微分時間[s]　とする。



●オイラー法、ルンゲ・クッタ法の適用できる式の形式
　２変数の連立１次常微分方程式
　　Δi1/Δt=f1(t,i1,i2)
　　Δi2/Δt=f2(t,i1,i2)
　が一般ですが、右辺に微分が含まれていても可能なのか？

改修案
　SECONDさんの工夫でＯＫ!?



私も正しい結果かどうか判断できませんが、波形を表示するプログラムを掲載します。

●その１
!回路図
!　　　→i1　　　 N:1　 →i2
!　　┌ｒ───・┐ ┌────┬
!　　E　　　 Vt↑L1 L2↑V2　　R2
!　　└─────┘ └・───┴

!抵抗rから変圧器までのFパラメータ（二端子対回路としてみると）
!　┌ 1　r ┐┌ N　0 　┐=┌ N　r/N ┐
!　└ 0　1 ┘└ 0　1/N ┘ └ 0　1/N ┘

!ただし、巻き数比N=SQR(L1/L2)とする。

!二端子対の両端の電圧、電流は
!　┌ E　┐=┌ N　r/N ┐┌ V2 ┐
!　└ i1 ┘ └ 0　1/N ┘└ i2 ┘
!となり、整理すると
!　┌ E　┐=┌ N*V2+r/N*i2 ┐ … ①
!　└ i1 ┘ └ 1/N*i2　　　┘ … ②

!二次側電圧V2=R2*i2より、i2=V2/R2 … ③

!③を①に代入して、V2=(N*R2)/(r+N^2*R2)*E … ④

!これを③に代入して、i2=N/(r+N^2*R2)*E … ⑤

!②に代入して、i1=E/(r+N^2*R2) … ⑥

!また、変圧器の一次側電圧Vt=E-r*i1 … ⑦ である。


LET E=10
LET L1=0.5
LET R1=1000
DEF r=R1*(1+SIN(w*t))
LET w=2*PI*50
LET R2=10
LET L2=0.1

LET N=SQR(L1/L2)
SUB routine
   LET V2=N*R2*E/(r+N^2*R2)
   LET i2=N*E/(r+N^2*R2)
   LET i1=E/(r+N^2*R2)
   LET Vt=E-r*i1
END SUB


LET yy=12 !縦軸の範囲

LET TT=50e-3 !時間区間 [0,T]


!!!SET bitmap SIZE 600,600 !画面を大きくする
SET WINDOW -TT/8,TT,-yy,yy !表示領域
DRAW grid(TT/5,yy/6)

PLOT TEXT ,AT TT*7/8,0: "[秒]"


SET LINE COLOR 1
PLOT LINES: t,Vt; t,Vt
SET LINE COLOR 4
PLOT LINES: t,i1; t,i1

FOR t=0 TO TT STEP TT/10000
   CALL routine
   !!!PRINT t;i1,V2;i2

   SET LINE COLOR 1
   PLOT LINES: t,Vt; t,Vt
   SET LINE COLOR 4
   PLOT LINES: t,i1; t,i1 !※
NEXT t
PLOT LINES


END


●その２
!回路図
!　　　　　　→i1　　 →i2
!　　┌ｒ───・┐ ┌────┬
!　　E　　　 V1↑L1 L2↑V2　　R2
!　　└─────┘ └・───┴

!　E-r*i1=V1=L1*Δi1/Δt-M*Δi2/Δt
!　　　　　　M*Δi1/Δt-L2*Δi2/Δt=V2=R2*i2

!２変数i1,i2の連立１次常微分方程式で表される。

!行列で表現すると
!　┌ E-r*i1 ┐=┌ L1　-M　┐┌ Δi1/Δt ┐
!　└ R2*i2　┘ └ M 　-L2 ┘└ Δi2/Δt ┘

!　┌ Δi1/Δt ┐=1/(-L1*L2+M*M)┌ -L2　M ┐┌ E-r*i1 ┐
!　└ Δi2/Δt ┘　　　　　　　 └ -M　L1 ┘└ R2*i2　┘

!ところで、M=SQR(L1*L2)より、上記の変形はできない。

!したがって、M=k*SQR(L1*L2)として考える。


!精度が悪いため、計算の後半では誤差が累積されやすくなる。



!２変数の連立１次常微分方程式の解

!----- ↓↓↓↓↓ -----

LET yy=12 !縦軸の範囲

LET TT=50e-3 !時間区間 [0,T]

DEF r(t)=R1*(1+SIN(w*t))
!DEF f1(t,i1,i2)=1/(L1*L2-M*M)*(L2*(E-r(t)*i1) - M*R2*i2)
!DEF f2(t,i1,i2)=1/(L1*L2-M*M)*(-M*(E-r(t)*i1) +L1*R2*i2)
DEF f1(t,i1,i2)=1/(-L1*L2+M*M)*(-L2*(E-r(t)*i1) + M*R2*i2)
DEF f2(t,i1,i2)=1/(-L1*L2+M*M)*( -M*(E-r(t)*i1) +L1*R2*i2)


LET E=10 ![V]
LET L1=0.5 ![H]
LET R1=1000 ![Ω]
LET w=2*PI*50 ![rad/s]
LET L2=0.1 ![H]
LET R2=10 ![Ω]

LET K=0.999 !理想変圧器 K=1　※
LET M=K*SQR(L1*L2) ![H]

!----- ↑↑↑↑↑ -----


!!!SET bitmap SIZE 600,600 !画面を大きくする
SET WINDOW -TT/8,TT,-yy,yy !表示領域
DRAW grid(TT/5,yy/6)

PLOT TEXT ,AT TT*7/8,0: "[秒]"


!●オイラー法
LET t=0 !初期値 i1(0)=0、i2(0)=0
LET i1=0
LET i2=0
SET LINE COLOR 1
PLOT LINES: t,E-r(t)*i1; t,E-r(t)*i1
SET LINE COLOR 4
PLOT LINES: t,i1; t,i1

LET h=TT/500000 !刻み幅　※
FOR i=1 TO 500000
   LET k1=h*f1(t,i1,i2)
   LET k2=h*f2(t,i1,i2)
   LET i1=i1+k1
   LET i2=i2+k2
   LET t=t+h

   !!!PRINT t,i1,i2

   SET LINE COLOR 1
   PLOT LINES: t,E-r(t)*i1; t,E-r(t)*i1

   SET LINE COLOR 4
   PLOT LINES: t,i1; t,i1
NEXT i
PLOT LINES


END
 

Re: 微分方程式の数値解法

 投稿者：山中和義  投稿日：2009年 4月14日(火)14時12分40秒
  > No.318[元記事へ]

しまむら1243さんへのお返事です。

> ２）電流波形は正弦波振動しません。

i=E/Rだから、分母のRを0～2*R（正弦波）としても、iは正弦波になりません。



SET WINDOW -1,5,-1,15
DRAW grid

LET E=1
LET R1=5
FOR t=0 TO 5 STEP 0.001
   LET r=R1*(1+SIN(5*t))
   LET i=E/r
   PLOT LINES: t,r; t,r !正弦波
   !PLOT LINES: t,i; t,i !正弦波？
NEXT t

END 		 

Re: 微分方程式の数値解法

 投稿者：しまむら1243  投稿日：2009年 4月14日(火)15時24分22秒
  > No.320[元記事へ]

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

山中さん、ご教示ありがとうございます。
確かにｒが正弦波状に変化しても、その逆数は正弦波状にはならないですね。根本的な誤りをしていました。
と言うことは、トランジスタの電流増幅機能（ベースに正弦波信号を入力すると、ほぼ正弦波状の電流出力が得られる）を模擬する可変抵抗は、正弦波の逆数で設定したほうが良いということになりそうですね。
SECONDさん、山中さんからご教示いただいたコードを、その観点で見てみます。

取り敢えずSECONDさんの原版の作画部分を再修正したものを下記に示します。
周波数を高くし、結合係数を小さくすると正弦波状に近づきます。
変圧器一次側電圧が±Eで振動するのも確認出来る様になりつつあります。

!※上のままでは、正帰還 暴走するので、f1(),f2() に、バッファ f1_,f2_ を付けた。
DEF f1(t, i1)=( E-R1*(1+SIN(w*t))*i1+M*f2_ )/L1 ! f2_=f2(t, i2)
DEF f2(t, i2)=( M*f1_-R2*i2 )/L2                 ! f1_=f1(t, i1)

LET E=10
LET R1=10000
LET ketugo=0.8 !結合係数で、1に近づけると波形が鋭利になり発振する
LET L1=0.2
LET L2=0.1
LET R2=2
LET frq=10000     !信号周波数[Hz]
LET ndiv=5000     !１サイクルの計算分割数
LET drawHz=10     ! 1画面で表示するサイクル数[Hz]
LET pass_gamen=1  !初期過渡時の描画面をパスする回数
LET iBairitu=1000 !電流波形描画拡大倍率。1000で1目盛り1[mA]になる

LET nmax=drawHz*ndiv !1画面の描画点数
LET dt=1/frq/ndiv    !計算微分時間[s]

LET M=ketugo*SQR(L1*L2)
LET w=2*PI*frq       !交流正弦波信号の角周波数[rad/s]
LET i1_=0            ! i1の初期値
LET i2_=0            ! i2の初期値
LET Vt_=0            !変圧器一次電圧Vtの初期値
LET t_=0             ! tの初期値

!----run

SET WINDOW -0,nmax,-2*E,2*E
DRAW grid

FOR nn=0 TO pass_gamen
   FOR n=0 TO nmax
      IF n=0 THEN LET xn_=0
      LET t=n*dt
      LET xn=n
      CALL RungeKutta4_1(t_,i1_,i1)
      CALL RungeKutta4_2(t_,i2_,i2)
      LET Vt=E-R1*(1+SIN(w*t))*i1

      IF nn=pass_gamen THEN
         SET LINE COLOR "red" !電流i1波形描画色
         PLOT LINES:xn_,i1_*iBairitu;xn,i1*iBairitu
         !SET LINE COLOR "blue" !電流i2波形描画色
         !PLOT LINES:xn_,i2_*iBairitu;xn,i2*iBairitu
         SET LINE COLOR "black" !電圧波形描画色
         PLOT LINES:xn_,Vt_;xn,Vt
         PLOT LINES: xn_,0;xn,0 !時間軸描画
         SET LINE COLOR "green" !電圧±E[V]ライン
         PLOT LINES:xn_,-E;xn,-E
         PLOT LINES:xn_,E;xn,E
      END IF

      LET t_=t   !　次のt へ更新
      LET i1_=i1 !　次のi1 へ更新
      LET i2_=i2 !　次のi2 へ更新
      LET Vt_=Vt !　次のVt へ更新
      LET xn_=xn
   next N
NEXT nn

SUB RungeKutta4_1(t_, i1_,i1)
   LET k1=f1(t_, i1_)
   LET k2=f1(t_+dt/2, i1_+k1*dt/2)
   LET k3=f1(t_+dt/2, i1_+k2*dt/2)
   LET k4=f1(t_+dt,   i1_+k3*dt )
   LET i1=i1_+(k1+2*k2+2*k3+k4)*dt/6
   LET f1_=f1(t_, i1)
END SUB

SUB RungeKutta4_2(t_, i2_,i2)
   LET k1=f2(t_, i2_)
   LET k2=f2(t_+dt/2, i2_+k1*dt/2)
   LET k3=f2(t_+dt/2, i2_+k2*dt/2)
   LET k4=f2(t_+dt,   i2_+k3*dt )
   LET i2=i2_+(k1+2*k2+2*k3+k4)*dt/6
   LET f2_=f2(t_, i2)
END SUB

END



> しまむら1243さんへのお返事です。
> > ２）電流波形は正弦波振動しません。
> i=E/Rだから、分母のRを0～2*R（正弦波）としても、iは正弦波になりません。
>
> SET WINDOW -1,5,-1,15
> DRAW grid
>
> LET E=1
> LET R1=5
> FOR t=0 TO 5 STEP 0.001
>    LET r=R1*(1+SIN(5*t))
>    LET i=E/r
>    PLOT LINES: t,r; t,r !正弦波
>    !PLOT LINES: t,i; t,i !正弦波？
> NEXT t
>
> END 		 

Re: 微分方程式の数値解法

 投稿者：SECOND  投稿日：2009年 5月 2日(土)17時18分22秒
  > No.315[元記事へ]

!どなたか御知恵 拝借 お願いします。

!先のルンゲクッタ法で、他の微分要素を含み、暴走を生じ、それを押えるために
!安定化バッファ、f1_ f2_ を用いたが、その是非を、別な等価回路と網回路で計算
!し直し、比較して確かめた。　一致する様で、正しい波形だった。
!計算間隔 0.0105秒 では互いの差( ルンゲ・クッタ１、ルンゲ・クッタ２ )が目立ち
!ますが、　～0.0005秒くらいで重なります。(全文コピー、RUN)
!
!問題は、ラプラス変換とも比較しようとした際に生じ、R1(b+sinωt)i(t)部分の変換で、
!時間関数の、抵抗と電流の積が、

! R1(b+sinωt)i(t) → R1*b*i(t)+R1*j/2*( e^(j*w*t)-e^(-j*w*t) )*i(t)
!      ラプラス変換 → R1*b*i(s)+R1*j/2*(  i(s-jw) - i(s+jw)   )

!のようになって、i(s) にまとまらず、i(s)を決定できずに、強引に jw=0 での i(s)に
! s-jw, s+jw を代入して求めた i(s) を逆変換して i(t) を出した。(添付画像の式参照）
!しかし・・

!R1=30 付近の低い値では、比較的一致するも、R1=100以上あたりから著しく外れる。
!原因は、強引な i(s)の前提ですが、参考例題も見当たらず難渋している。

!下の laplace_ON を１にすると、ルンゲ・クッタ２ との比較テストモードに入ります。
!ルンゲ・クッタ１は、オーバーフローしやすい為、除いています。

!問題を簡単にするため、２次回路開放同然 R2→1e9 で、
! E, L1, R1(1+sinωt) の３素子のみの直列回路として、その電流i(t)と
!                                   R1(1+sinωt)の、両端電圧v(t)を、描きます。
! ラプラス変換　　　v(t)灰, i(t)緑
! ルンゲ・クッタ２　v(t)赤, i(t)青  … ２次側は単独、L2開放端v2(t)赤, i2(t)青≒０
!---------------------------
!
OPTION ARITHMETIC COMPLEX
LET laplace_ON=0

!---------------------------ルンゲ・クッタ１で計算。(バッファが無いと暴走する)
!回路の微分方程式 　　r(t)=R1*(1+sin(w*t))
!
!  i1┌→┐M┌→┐i2    （一次側の電圧平衡式）L1*(di1/dt)+r(t)*i1-M*(di2/dt)=E
!    ┌─┐ ┌─┐      （二次側の電圧平衡式）M*(di1/dt)=L2*(di2/dt)+R2*i2
!    │  L1 L2  R2
!    E   │ └─┘
!    │  r(t)=R1*(1+sin(w*t))
!    └─┘
!  バッファ無しでは    f1(t, i1)= (di1/dt)= ( E-R1*(1+sin(w*t))*i1+M*(di2/dt) )/L1
!  暴走する 微分方程式 f2(t, i2)= (di2/dt)= ( M*(di1/dt)-R2*i2 )/L2

!---------------------------ルンゲ・クッタ２で計算。(バッファ不要)
!相互インダクタンス M= k*(√L1*L2)の結合係数 k=1 として、漏れ磁束の無い場合、
!回路のi1側に揃えるＴ型等価回路で、網路を変更。
!
!    ┌ i31→┐         上の回路と等価な i2 =(i31+i21)*SQR(L1/L2)
!        ┌→┐i21
!    ┌─┬─┐                 (i31の網路) R2*L1/L2*(i31+i21) +r(t)*i31= E
!    │  L1  R2*L1/L2           (i21の網路) R2*L1/L2*(i31+i21) +L1*(di21/dt)= 0
!    E   ├─┘
!    │  r(t)=R1*(1+sin(w*t))
!    └─┘
!                                              i31= (E-R2*L1/L2*i21) /(r(t)+R2*L1/L2)
!                              E+L1*(di21/dt)= r(t)*i31
!                          　　E+L1*(di21/dt)= r(t)*(E-R2*L1/L2*i21) /(r(t)+R2*L1/L2)
!                                    (
!                                     )
!問題の無い 微分方程式 f21(t, i21)= (di21/dt)= -(E+r(t)*i21) /(L2/R2*r(t)+L1)

!-------------------------------
LET E=10
LET R1=1000
LET L1=20
LET L2=20
LET R2=1000
LET M=SQR(L1*L2) !結合係数 k=1
!
LET b=1
LET w=2*PI*1 !1Hz
DEF r(t)=R1*(b+SIN(w*t))
!
IF laplace_ON=1 THEN
   LET R1=30
   LET R2=1e9
END IF
!
!-----ルンゲ・クッタ１
DEF f1(t,i1)=( E-r(t)*i1+M*f2_ )/L1 ! f2_…直接のf2()は、不可
DEF f2(t,i2)=( M*f1_-R2*i2 )/L2     ! f1_…直接のf1()は、不可

SUB RungeKutta_1(t,i1,i2)
   LET k1=f1(t,i1)
   LET k2=f1(t+dt/2, i1+k1*dt/2)
   LET k3=f1(t+dt/2, i1+k2*dt/2)
   LET k4=f1(t+dt,   i1+k3*dt )
   LET i1=i1+(k1+2*k2+2*k3+k4)*dt/6
   LET f1_=f1(t,i1) ! f1_…f1()のバッファ
   !
   LET k1=f2(t,i2)
   LET k2=f2(t+dt/2, i2+k1*dt/2)
   LET k3=f2(t+dt/2, i2+k2*dt/2)
   LET k4=f2(t+dt,   i2+k3*dt )
   LET i2=i2+(k1+2*k2+2*k3+k4)*dt/6
   LET f2_=f2(t,i2) ! f2_…f2()のバッファ
END SUB

!-----ルンゲ・クッタ２
DEF i31=(E-R2*L1/L2*i21)/(r(t)+R2*L1/L2)
DEF f21(t,i21)=-(E+r(t)*i21)/(L2/R2*r(t)+L1)
DEF i20=(i31+i21)*SQR(L1/L2)

SUB RungeKutta_2(t,i21)
   LET k1=f21(t,i21)
   LET k2=f21(t+dt/2, i21+k1*dt/2)
   LET k3=f21(t+dt/2, i21+k2*dt/2)
   LET k4=f21(t+dt,   i21+k3*dt )
   LET i21=i21+(k1+2*k2+2*k3+k4)*dt/6
END SUB

!-----ラプラス逆変換
DIM Xr(5)
!   xr(0)=0 …根は６個、０は自明で因数分解から除いた。
LET Xr(1)=-R1*b/L1
LET Xr(2)=COMPLEX(0,w)
LET Xr(3)=COMPLEX(-R1*b/L1,w)
LET Xr(4)=COMPLEX(0,-w)
LET Xr(5)=COMPLEX(-R1*b/L1,-w)
FOR j=1 TO 5
   PRINT Xr(j)
NEXT j
DEF Gs(s)=E/L1*( (s^2+w^2)*((s+R1*b/L1)^2+w^2) -R1/L1*s*w*(2*s+R1*b/L1) )*EXP(s*t)
DEF k_0=Gs(0    )/ (       (0    -Xr(1))*(0    -Xr(2))*(0    -Xr(3))*(0    -Xr(4))*(0    -Xr(5)) )
DEF k_1=Gs(Xr(1))/ ( Xr(1)              *(Xr(1)-Xr(2))*(Xr(1)-Xr(3))*(Xr(1)-Xr(4))*(Xr(1)-Xr(5)) )
DEF k_2=Gs(Xr(2))/ ( Xr(2)*(Xr(2)-Xr(1))              *(Xr(2)-Xr(3))*(Xr(2)-Xr(4))*(Xr(2)-Xr(5)) )
DEF k_3=Gs(Xr(3))/ ( Xr(3)*(Xr(3)-Xr(1))*(Xr(3)-Xr(2))              *(Xr(3)-Xr(4))*(Xr(3)-Xr(5)) )
DEF k_4=Gs(Xr(4))/ ( Xr(4)*(Xr(4)-Xr(1))*(Xr(4)-Xr(2))*(Xr(4)-Xr(3))              *(Xr(4)-Xr(5)) )
DEF k_5=Gs(Xr(5))/ ( Xr(5)*(Xr(5)-Xr(1))*(Xr(5)-Xr(2))*(Xr(5)-Xr(3))*(Xr(5)-Xr(4))               )
DEF k0_5=re(k_0+k_1+k_2+k_3+k_4+k_5) !留数の和

!-----run
SET TEXT background "OPAQUE"
DIM baki(10),bakv(10) ! drawing channels
!
FOR dt=0.0105 TO 0.000499 STEP -0.001 !演算ピッチ sec. pitch time
   IF laplace_ON=1 THEN LET dt=.0005
   SET DRAW mode hidden
   CLEAR
   SET DRAW mode explicit
   LET t=0 !計算開始 sec. calculation start
   LET ts=0 !描画開始 sec. drawing start
   LET tw=3 !描画時間 sec. drawing time
   LET Vw=28 !＋目盛上限(上側グラフ)１の倍数、Volt.Ampere. maximum scale
   LET ofs=15 !＋目盛上限(下側グラフ)５の倍数、Under_Graph offset !CEIL(Vw/10)*5
   !
   LET i21= -E/R1   ! ルンゲ・クッタ２ 初期値( i1=i31←i21, i2=i20←i21+i31) at t=0
   LET i1=i31       !┐
   LET i2=i20       !┼ルンゲ・クッタ１ 初期値! at t=0
   LET f2_=f2(t,i2) !┘
   !
   IF laplace_ON=1 THEN LET i21=0
   !
   SET WINDOW -.06*tw, tw, -Vw,Vw
   SET COLOR MIX(15) .5,.5,.5
   DRAW grid(1,5)
   PLOT TEXT,AT .1,Vw*0.92,USING"計算間隔=#.###### 秒。　計算開始後###.### 秒からの描画":dt,ts
   PLOT TEXT,AT .1,Vw*0.86,USING"ω=#.##rad/s E=##V L1=##H L2=##H k=1 R1=####Ω R2=####Ω":w,E,L1,L2,R1,R2
   PLOT TEXT,AT .1,-0.1*Vw:"r(t)="& STR$(R1)& "*("& STR$(b)& "+sinωt)"
   DO
      IF laplace_ON=1 THEN
      !-----ラプラス
         CALL Grph( "green","i1(x20mA)", K0_5*50, "gray","v1(x1V)", K0_5*r(t) , 8, 0)
      ELSE
      !-----ルンゲ・クッタ１
         CALL Grph( "blue","i1(x20mA)", i1*50 , "red","v1(x1V)",r(t)*i1          , 1, 0)
         CALL Grph( "blue","i1(x20mA)", i1*50 , "red","v1(x1V)", E-(L1*f1_-M*f2_), 2, 0)
         CALL Grph( "blue","i2(x2mA)" ,-i2*500, "red","v2(x1V)", L2*f2_-M*f1_    , 3, ofs )
         CALL Grph( "blue","i2(x2mA)" ,-i2*500, "red","v2(x1V)",-R2*i2           , 4, ofs )
         CALL RungeKutta_1(t,i1,i2)
      END IF
      !-----ルンゲ・クッタ２
      CALL Grph( "blue","i1(x20mA)", i31*50 , "red","v1(x1V)", r(t)*i31       , 5, 0)
      CALL Grph( "blue","i1(x20mA)", i31*50 , "red","v1(x1V)", E+L1*f21(t,i21), 6, 0)
      CALL Grph( "blue","i2(x2mA)" ,-i20*500, "red","v2(x1V)",-R2*i20         , 7, ofs)
      CALL RungeKutta_2(t, i21)
      !-----
      LET t=t+dt
   LOOP UNTIL ts+tw< t
NEXT dt

!-----draw
SUB Grph( icol$,i$,i, vcol$,v$,v, n, ofs)
   SET WINDOW -.06*tw, tw, -Vw+ofs, Vw+ofs
   IF ts< t THEN
      SET LINE COLOR vcol$
      PLOT LINES :t-ts-dt, bakv(n); t-ts, v
      SET LINE COLOR icol$
      PLOT LINES :t-ts-dt, baki(n); t-ts, i
   ELSEIF t=ts THEN
      SET TEXT COLOR vcol$
      PLOT TEXT,AT .1, .17*Vw :v$
      SET TEXT COLOR icol$
      PLOT TEXT,AT .1, .11*Vw :i$
      IF ofs<>0 THEN
         SET LINE COLOR 15
         PLOT LINES :-.06*tw, 0; tw,0
         SET TEXT COLOR 15
         ASK PIXEL SIZE (0,0 ;tw,Vw) px,py
         FOR j=5*INT((-Vw+ofs)/5+1) TO ofs-1 STEP 5
            PLOT TEXT,AT -22*tw/px, j-15*Vw/py, USING"###":j
         NEXT j
      END IF
      SET TEXT COLOR 1
   END IF
   LET baki(n)=i
   LET bakv(n)=v
END SUB

END
 

Re: 微分方程式の数値解法

 投稿者：しまむら1243  投稿日：2009年 5月 5日(火)21時04分16秒
  > No.343[元記事へ]

SECONDさんへのお返事です。

SECONDさん、いろいろご検討頂いて有り難うございます。

!---------------------------ルンゲ・クッタ２で計算。(バッファ不要)
!相互インダクタンス M= k*(√L1*L2)の結合係数 k=1 として、漏れ磁束の無い場合、
!回路のi1側に揃えるＴ型等価回路で、網路を変更。
!問題の無い 微分方程式 f21(t, i21)= (di21/dt)= -(E+r(t)*i21) /(L2/R2*r(t)+L1)

変圧器の一次側換算等価回路で考えるこの案、いいですね。
原式にしたがって忠実に計算するのが一番正確だと拘っていて気付きませんでした。

ただ一寸電気的に理解できないのは、等価回路方式ではk=1の場合を扱えるのに、原式から行列を使ってルンゲ・クッタ法を適用できる連立微分式を導かれた山中さんの式だとk=1の場合が扱えない、という点です。等価回路も山中さんの連立微分式も、同じ原式から組み立てられたものなのですが。。でもこれはプログラミングとは別の問題ですね。

さてラプラス変換を利用してルンゲ・クッタの修正法を検証しようと

!のようになって、i(s) にまとまらず、i(s)を決定できずに、強引に jw=0 での i(s)に
! s-jw, s+jw を代入して求めた i(s) を逆変換して i(t) を出した。(添付画像の式参照）

とされていますが、t=＋0で（は未だ交流になっていないから）jw=0とする、という手法はsin(wt)のラプラス変換の原型が崩されてしまうので、理屈的に不合理になって無理ではないかと思います。
では、何か巧い方法は？と言われると、ラプラス変換は定係数の微分式に有用であって、変係数の微分式に適用出来ないのではないでしょうか。残念ながら手元の工学用数学書とネット検索した限りではその様な記載しか無かったです。

私は山中さんが導かれた連立微分式にルンゲ・クッタ法を適用したものと、SECONDさんが作成されたルンゲ・クッタ修正法とを同じ条件で計算してみましたが、t=＋0の立ち上がり付近でSECONDさんの波形には高周波振動が含まれていた点を除けば、ぴったり合っていました。 		 

Re: 微分方程式の数値解法

 投稿者：SECOND  投稿日：2009年 5月 5日(火)22時50分21秒
  > No.345[元記事へ]

しまむら1243さんへのお返事です。

こんにちは、調査して頂いたようで、大変感謝します。でも今の私には、

> では、何か巧い方法は？と言われると、ラプラス変換は定係数の微分式に有用であって、変係数の微分式に適用出来ないのではないでしょうか。残念ながら手元の工学用数学書とネット検索した限りではその様な記載しか無かったです。

この文だけしか、目に入りません。本当に不可能なのか否か、にわかには受け入れがたく
もうしばらくは、探索にこだわってみます。ありがとうございました。 		 

戻る