上記システムは、u(t)
という一入力、y(t) と φ(t) の二出力を持
つ非線形システムである。
ここで、システムのパラメー タが以下である場合を考える。
m=4.735×10-3[kg]、M=0.628
[kg]、λ=1.003×10-4[kg・m2/s]、
υ=4.01[kg/s]、L=0.2465[m]、I=1.155×10-4[kg・m2]、
G=3.18[N/V]、
ここで、初期値として、y(0)=0,dy(0)/dt=0, φ(0)=1, dφ(0)/dt=0, としたときの、このシステムの自由運動をシミュレーションした結果が次のグラフである。この結果を得るためのソフトはその下に示すものである。
Sub touritsu1()
Dim A(1 To 2, 1 To 2) As Double, Ainv(1 To 2, 1 To 2) As Double
Dim b(1 To 2, 1 To 2) As Double, dummy(1 To 2) As Double
Dim time(1001) As Double, dt As Double
Dim y(1001) As Double, dydt(1001) As Double, ddyddt(1001)
Dim phi(1001) As Double, dphidt(1001) As Double, ddphiddt(1001)
Dim u(1001) As Double
Dim mp As Double, mc As Double, lambda As Double, ny As Double
Dim leng As Double, inam As Double, ga As Double, g As Double
Dim i As Integer
Set WS = Worksheets("Sheet1")
WS.Activate
Range("A1:I1300").Select
Selection.Clear
mp = 0.004735
mc = 0.628
lambda = 0.0001003
ny = 4.01
leng = 0.2465
inam = 0.0001155
ga = 3.18
g = 9.8
dt = 0.005
For i = 0 To 1000
time(i) = i * dt
u(i) = 0
Next
phi(0) = 1
dphidt(0) = 0
y(0) = 0
dydt(0) = 0
ddphiddt(0) = 0
ddyddt(0) = 0
For i = 1 To 1000
phi(i) = phi(i - 1) + dphidt(i - 1) * dt
y(i) = y(i - 1) + dydt(i - 1) * dt
dphidt(i) = dphidt(i - 1) + ddphiddt(i - 1) * dt
dydt(i) = dydt(i - 1) + ddyddt(i - 1) * dt
A(1, 1) = mp + mc
A(1, 2) = mp * leng * Cos(phi(i - 1))
A(2, 1) = mp * leng * Cos(phi(i - 1))
A(2, 2) = inam + mp * leng * leng
Ainv(1, 1) = A(2, 2) / (A(1, 1) * A(2, 2) - A(1, 2) *
A(2, 1))
Ainv(2, 2) = A(1, 1) / (A(1, 1) * A(2, 2) - A(1, 2) *
A(2, 1))
Ainv(1, 2) = -A(1, 2) / (A(1, 1) * A(2, 2) - A(1, 2) *
A(2, 1))
Ainv(2, 1) = -A(2, 1) / (A(1, 1) * A(2, 2) - A(1, 2) *
A(2, 1))
dummy(1) = Ainv(1, 1) * (mp * leng * Sin(phi(i - 1)) * dphidt(i - 1)
* dphidt(i - 1) - ny * dydt(i - 1)) + Ainv(1, 2) * (-lambda * dphidt(i
- 1) + mp * leng * g * Sin(phi(i - 1))) + Ainv(1, 1) * ga * u(i - 1)
dummy(2) = Ainv(2, 1) * (mp * leng * Sin(phi(i - 1)) * dphidt(i - 1)
* dphidt(i - 1) - ny * dydt(i - 1)) + Ainv(2, 2) * (-lambda * dphidt(i
- 1) + mp * leng * g * Sin(phi(i - 1))) + Ainv(2, 1) * ga * u(i - 1)
ddyddt(i) = dummy(1)
ddphiddt(i) = dummy(2)
Cells(3 + i, 12) = Ainv(1, 1)
Cells(3 + i, 13) = Ainv(1, 2)
Cells(3 + i, 14) = Ainv(2, 1)
Cells(3 + i, 15) = Ainv(2, 2)
Next
Cells(1, 2) = "時刻t(sec)"
Cells(2, 3) = "y(t)×1000"
Cells(2, 4) = "Φ(t)"
For i = 0 To 1000
Cells(3 + i, 2) = time(i)
Cells(3 + i, 3) = 1000 * y(i)
Cells(3 + i, 4) = phi(i)
Next
Range("B2:D1003").Select
Charts.Add
ActiveChart.ChartType = xlLineMarkers
ActiveChart.SetSourceData Source:=Sheets("Sheet1").Range("B2:D1003"),
PlotBy _
:=xlColumns
ActiveChart.Location Where:=xlLocationAsObject,
Name:="Sheet1"
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
= "y(t)×1000[m],Φ(t)[rad]"
End With
End Sub
（注） m-adachi@02.246.ne.jp へメールで依頼があれば、上記エクセルのソフトをメールの添付ファイルで送ります。
簡単な常微分方程式の解を求めるシミュレーションへ