簡単な常微分方程式の解を求めるシミュレーション

　　ここでは、

１）　ラプラス変換を使って解析解を求めることが大変な、高次の線形常微分方程式系の解
２）　むだ時間系のステップ応答
３）　特性方程式の根が実数であるシステムと複素数であるシステムのステップ応答
４）　インパルス入力信号に対する二次遅れ系の応答

を数値シミュレーションで求めることにしょう。

（注）　むだ時間系とは、その系（システム）へ入力しても、応答が出てくるまでに有限時間の遅れが存在するシステムのことである。イメージとしては、衛星回線を通じて放送されるテレビ映像の場合、伝送の時間遅れのため会話が相手の話にかぶったりして、会話の応答が不自然になる場合がる。これは、まさしくむだ時間系であり、相手側の発言入力に対して、こちら側が聞き取るまでに遅れ時間が発生している。
（注）　線形常微分方程式系　 a_n*dⁿx(t)/dtⁿ + a_n-1*d^n-1x(t)/dt^n-1 +　・・・　+ a₁*dx(t)/dt + a₀*x(t) = u(t),　
の特性方程式とは、微分を　s　で置き換え,uを省いた次のｓの多項式のことである。
　　　　 a_n*sⁿ + a_n-1*s^n-1 +　・・・　+ a₁*s + a₀ = 0,　

１）　高次のシステムとして、以下を考える。
d^３x(t)/dt^３ + 4*d^２x(t)/dt^２ + 3*dx(t)/dt + x(t) = u(t),

入力： u(t) = 1, for　t≧0
　　　　 u(t) = 0, for t＜0
初期状態： d^２x(0)/dt^２=0, dx(0)/dt=0,   x(0) = 0,

このシステムのｘを求める計算、およびそのグラフを求めるソフトは以下である。

Sub sim2()
Dim x(1001) As Double, u(1001) As Double, dxdt(1001) As Double
Dim ddxddt(1001) As Double, dddxdddt(1001) As Double
Dim time(1001) As Double, dt As Double
Dim i As Integer

   Set WS = Worksheets("Sheet2")
   WS.Activate
   Range("A1:I1300").Select
   Selection.Clear

Cells(1, 2) = "時刻t(sec)"
Cells(2, 3) = "x(t)"

dt = 0.02
time(0) = 0
u(0) = 0
For i = 1 To 1000
time(i) = dt * i
u(i) = 1
Next

ddxddt(0) = 0
dxdt(0) = 0
x(0) = 0

For i = 1 To 1000
x(i) = x(i - 1) + dxdt(i - 1) * dt
dxdt(i) = dxdt(i - 1) + ddxddt(i - 1) * dt
ddxddt(i) = ddxddt(i - 1) + dddxdddt(i - 1) * dt
dddxdddt(i) = -4 * ddxddt(i) - 3 * dxdt(i) - x(i) + u(i)
Next
For i = 0 To 1000
Cells(3 + i, 2) = time(i)
Cells(3 + i, 3) = x(i)
Next

    Range("B2:C1003").Select
    Charts.Add
    ActiveChart.ChartType = xlLineMarkers
    ActiveChart.SetSourceData Source:=Sheets("Sheet2").Range("B2:C1003"), PlotBy _
        :=xlColumns
    ActiveChart.Location Where:=xlLocationAsObject, Name:="Sheet2"
    With ActiveChart
        .HasTitle = True
        .ChartTitle.Characters.Text = "三次遅れ系のステップ応答波形"
        .Axes(xlCategory, xlPrimary).HasTitle = True
        .Axes(xlCategory, xlPrimary).AxisTitle.Characters.Text = "time(sec)"
        .Axes(xlValue, xlPrimary).HasTitle = False
    End With

End Sub

このときのグラフは、次である。

２）　むだ時間系として、以下を考える。
y(t) = x(t-0.2),
dx(t)/dt + x(t) = u(t),
入力： u(t) = 1, for　t≧0
　　　　 u(t) = 0, for t＜0
初期状態：   x(0) = 0,

このシステムのｘを求める計算、およびそのグラフを求めるソフトは以下である。

Sub sim3()

Dim x(1001) As Double, u(1001) As Double
Dim time(1001) As Double, dt As Double
Dim y(1001) As Double, dydt(1001) As Double
Dim z(1001) As Double, dzdt(1001) As Double
Dim i As Integer

   Set WS = Worksheets("Sheet3")
   WS.Activate
   Range("A1:I1300").Select
   Selection.Clear

Cells(1, 2) = "時刻t(sec)"
Cells(2, 3) = "y(t)"
Cells(2, 4) = "x(t)"

dt = 0.005
time(0) = 0
u(0) = 0
x(0) = 0
For i = 1 To 1000
u(i) = 1
time(i) = dt * i
Next

For i = 1 To 1000
x(i) = x(i - 1) + (u(i - 1) - x(i - 1)) * dt
Next
For i = 0 To 40
y(i) = 0
Next
For i = 41 To 1000
y(i) = x(i - 40)
Next
For i = 0 To 1000
Cells(3 + i, 2) = time(i)
Cells(3 + i, 3) = y(i)
Cells(3 + i, 4) = x(i)
Next

    Range("B2:D1003").Select
    Charts.Add
    ActiveChart.ChartType = xlLineMarkers
    ActiveChart.SetSourceData Source:=Sheets("Sheet3").Range("B2:D1003"), PlotBy _
        :=xlColumns
    ActiveChart.Location Where:=xlLocationAsObject, Name:="Sheet3"
    With ActiveChart
        .HasTitle = True
        .ChartTitle.Characters.Text = "一次遅れ＋むだ時間系のステップ応答波形"
        .Axes(xlCategory, xlPrimary).HasTitle = True
        .Axes(xlCategory, xlPrimary).AxisTitle.Characters.Text = "time(sec)"
        .Axes(xlValue, xlPrimary).HasTitle = False
    End With

End Sub

このときのグラフは、次である。

３）　特性方程式の根が実数の場合と、複素数の場合の応答の違い。自由応答で初期値を非ゼロの場合、

d^２x(t)/dt^２ + 1.1*dx(t)/dt + 0.3*x(t) = 0,　　（特性方程式の根は、-0.5, -0.6）

初期状態： dx(0)/dt=0,   x(0) = 1,

d^２y(t)/dt^２ + 1.1*dy(t)/dt + 6*y(t) = 0,     （特性方程式の根は、-0.55±4.77*j）

初期状態： dy(0)/dt=0,   y(0) = 1,

このシステムのｘを求める計算、およびそのグラフを求めるソフトは以下である。

Sub sim4()
Dim x(1001) As Double, dxdt(1001) As Double
Dim time(1001) As Double, dt As Double
Dim y(1001) As Double, dydt(1001) As Double
Dim i As Integer

   Set WS = Worksheets("Sheet4")
   WS.Activate
   Range("A1:I1300").Select
   Selection.Clear

Cells(1, 2) = "時刻t(sec)"
Cells(2, 3) = "x(t)"
Cells(2, 4) = "y(t)"

time(0) = 0
dt = 0.01
dxdt(0) = 0
x(0) = 1
dydt(0) = 0
y(0) = 1

For i = 1 To 1000
time(i) = dt * i
x(i) = x(i - 1) + dxdt(i - 1) * dt
dxdt(i) = dxdt(i - 1) + (-1.1 * dxdt(i - 1) - 0.3 * x(i - 1)) * dt
y(i) = y(i - 1) + dydt(i - 1) * dt
dydt(i) = dydt(i - 1) + (-1.1 * dydt(i - 1) - 6 * y(i - 1)) * dt
Next
For i = 0 To 1000
Cells(3 + i, 2) = time(i)
Cells(3 + i, 3) = x(i)
Cells(3 + i, 4) = y(i)
Next

    Range("B2:D1003").Select
    Charts.Add
    ActiveChart.ChartType = xlLineMarkers
    ActiveChart.SetSourceData Source:=Sheets("Sheet4").Range("B2:D1003"), PlotBy _
        :=xlColumns
    ActiveChart.Location Where:=xlLocationAsObject, Name:="Sheet4"
    With ActiveChart
        .HasTitle = True
        .ChartTitle.Characters.Text = "実根と複素根を持つ各系の自由応答波形の相違"
        .Axes(xlCategory, xlPrimary).HasTitle = True
        .Axes(xlCategory, xlPrimary).AxisTitle.Characters.Text = "time(sec)"
        .Axes(xlValue, xlPrimary).HasTitle = False
    End With

End Sub

このときのグラフは、次である。

上記のグラフから、特性方程式の根が複素数のとき、出力が振動成分を含むことが分かる。

４）　二次遅れ系　　　dx/dt + x(t)=y(t), d(y)/dt + 2*y(t) = u(t),　初期条件： y(0)=0,   x(0) = 0,　に対して
　　　インパルス信号　u(t) = δ(t), for　t=0　 u(t) = 0, 　　 for t≠0　　　（単位インパルス関数）
が入力された場合の　x(t)　の応答出力は次のソフトを走らせて求められる。

Sub sim5()
Dim x(10001) As Double, u(10001) As Double, dxdt(10001) As Double
Dim y(10001) As Double, dydt(10001) As Double
Dim time(10001) As Double, dt As Double
Dim i As Integer

   Set WS = Worksheets("Sheet5")
   WS.Activate
   Range("A1:I1300").Select
   Selection.Clear

Cells(1, 2) = "時刻t(sec)"
Cells(2, 3) = "x(t)"

dt = 0.001
time(0) = 0
For i = 0 To 9
time(i) = dt * i
u(i) = 0.1 / dt
Next
For i = 10 To 10000
time(i) = dt * i
u(i) = 0
Next

y(0) = 0
x(0) = 0

For i = 1 To 10000
x(i) = x(i - 1) + (y(i - 1) - x(i - 1)) * dt
y(i) = y(i - 1) + (u(i - 1) - 2 * y(i - 1)) * dt
Next
For i = 0 To 1000
Cells(3 + i, 2) = time(i * 10)
Cells(3 + i, 3) = x(i * 10)
Cells(3 + i, 6) = Abs(Exp(-dt * i * 10) - Exp(-2 * dt * i * 10) - x(i * 10))
Next

    Range("B2:C1003").Select
    Charts.Add
    ActiveChart.ChartType = xlLineMarkers
    ActiveChart.SetSourceData Source:=Sheets("Sheet5").Range("B2:C1003"), PlotBy _
        :=xlColumns
    ActiveChart.Location Where:=xlLocationAsObject, Name:="Sheet5"
    With ActiveChart
        .HasTitle = True
        .ChartTitle.Characters.Text = "二次遅れ系のインパルス応答波形"
        .Axes(xlCategory, xlPrimary).HasTitle = True
        .Axes(xlCategory, xlPrimary).AxisTitle.Characters.Text = "time(sec)"
        .Axes(xlValue, xlPrimary).HasTitle = True
        .Axes(xlValue, xlPrimary).AxisTitle.Characters.Text = "x(t)"
    ActiveChart.HasLegend = False
End With

End Sub

（注）　m-adachi@02.246.ne.jp へメールで依頼があれば、上記エクセルのソフトをメールの添付ファイルで送ります。

数値シミュレーションの最初のページへ

PID制御のシミュレーションへ

ホームページのトップへ