|
|
しばっちさんへのお返事です。
> !' dy^m/dx^m=F(x,dy(0),dy(1),dy(2),...,dy(m-1))
しばっちさんが作られたM階常微分方程式のRungeKutta法プログラムを、RLC直列電気回路の過渡現象解析(M=2)に応用し、厳密解との差がどの程度あるかを確認する下記プログラムを書いてRUNしたら、凄い精度（分割数50程度でも一致）であることが分かりました。
しばっちさんの書かれたプログラムは汎用として使えるのでとても良いですね！！
-----以下は確認したプログラムコード----
!しばっちさんの投稿プログラム(数値解法)にグラフ表示を追加し
!RLC直列回路に交流電圧vs=V*cos(wt+fai)を投入した場合の過渡電流
!計算に適用し、その数値解とラプラス変換に依る厳密解とを比較したもの。
OPTION BASE 0
OPTION ARITHMETIC COMPLEX
!--- RungeKutta法を適用する２階常微分方程式の設定 -----
!RLC直列回路の連立微分方程式は、vcをCの端子電圧、icを電流とすると、
! L*dic/dt+R*ic+vc=vs と ic=C*dvc/dt
!となるので、上記連立方程式をvcの２階微分方程式に纏めたものが下記FUNCTION
! di^m/dt^m=F(t,di(0),di(1),di(2),...,di(m-1))の設定
FUNCTION F(t,DY()) !' DY(0)=i , DY(1)=i' , DY(2)=i''...
LET F=V*COS(w*t+fai)/CL-R/L*DY(1)-DY(0)/CL
END FUNCTION
!--- 共通計算定数・条件の設定 ----
LET M=2 !常微分方程式の階数
LET Ts=0 ![msec] Ts～Teまで計算
LET Te=20 ![msec]
LET N=100 !Ts～Teまでの分割数
LET V=1 !投入電圧[p.u.]
LET fai=45 !Vの投入位相角[度]
LET Hz=50 !Vの周波数[Hz]
LET R=10 ![Ω]
LET L=0.5 ![H]
LET C=4.946E-4 ![F]
!---定数の内部処理---
LET Ts=Ts/1000 ![sec]に換算
LET Te=Te/1000
LET CL=C*L
LET w=2*PI*Hz
LET fai=RAD(fai)
LET Im=V/SQR(R^2+(w*L-1/w/C)^2) !定常電流の最大値
SET WINDOW 0,N,-1.8,1.8
PLOT LINES:0,0;N,0
DIM DY(M-1),J1(M),J2(M),J3(M),J4(M)
DIM Y1(M-1),Y2(M-1),Y3(M-1)
LET h=(Te-Ts)/N
LET t=Ts
LET DY(0)=0 !'iの初期値
pause "最初に、RungeKutta法による数値解を赤色線で表示します。"
SET POINT STYLE 7
SET POINT COLOR "black"
PLOT POINTS:0,V*COS(fai)
SET POINT COLOR "red"
PLOT POINTS:0,DY(0)
FOR I=1 TO N
FOR K=1 TO M-1
LET J1(K)=h*DY(K)
NEXT K
LET J1(M)=h*F(t,DY)
FOR K=1 TO M-1
LET J2(K)=h*(DY(K)+J1(K+1)/2)
NEXT K
FOR K=0 TO M-1
LET Y1(K)=DY(K)+J1(K+1)/2
NEXT K
LET J2(M)=h*F(t+h/2,Y1)
FOR K=1 TO M-1
LET J3(K)=h*(DY(K)+J2(K+1)/2)
NEXT K
FOR K=0 TO M-1
LET Y2(K)=DY(K)+J2(K+1)/2
NEXT K
LET J3(M)=h*F(t+h/2,Y2)
FOR K=1 TO M-1
LET J4(K)=h*(DY(K)+J3(K+1))
NEXT K
FOR K=0 TO M-1
LET Y3(K)=DY(K)+J3(K+1)
NEXT K
LET J4(M)=h*F(t+h,Y3)
LET t=t+h
FOR K=0 TO M-1
LET DY(K)=DY(K)+(J1(K+1)+2*J2(K+1)+2*J3(K+1)+J4(K+1))/6
NEXT K
SET POINT COLOR "black"
PLOT POINTS:i,V*COS(w*t+fai)
LET ic=C*DY(1)
SET POINT COLOR "red"
PLOT POINTS:i,ic/Im
NEXT I
!--- ラプラス変換法に依る厳密解計算式 ---
pause "続けて、厳密解を緑色線で表示しますか？"
LET j=SQR(-1)
LET dt =(Te-Ts)/N !計算微分時間[秒]を設定
LET Zst=SQR(R^2+(w*L-1/w/C)^2)
LET Ist=V/Zst
LET alfa = -R / 2 /L
LET beta = SQR(1/C /L - alfa^2)
LET delta1=alfa+j*beta
LET delta2=alfa-j*beta
LET K1 =delta1/(delta1-delta2)/(delta1-j*w)
LET K2 =delta2/(delta2-delta1)/(delta2-j*w)
LET K3 =j*w/(j*w-delta1)/(j*w-delta2)
PRINT "【１】ラプラス変換法により得られた定数を表示します。"
PRINT "α=";alfa;" , β=";beta
PRINT "K1=";K1
PRINT "K2=";K2
PRINT "K3=";K3
PRINT
FOR nn=0 TO N
LET t = nn * dt
LET wt=w*t
LET i1 =K1* EXP(delta1*t)
LET i2 =K2* EXP(delta2*t)
LET i3 =K3*EXP(j*wt)
LET i = EXP(j*fai)*V/L*(i1 +i2 + i3)
LET ipu=Re(i)/Ist
LET vpu=V*COS(wt+fai)
!*** 描画指示 ***
SET POINT COLOR "black"
PLOT POINTS: nn,vpu !電源波形
SET POINT COLOR "green"
PLOT POINTS: nn,ipu !電流波形
NEXT nn
PRINT "●黒色曲線は、印加電源電圧の波形です。"
PRINT "●数値解と厳密解が一致した場合は、赤色線の上に緑色線が上書き"
PRINT "され、緑色線のみが見える事に留意してください。"
END
|
|