MaTXで微分方程式の数値解法を行う。

作成：2000-04-20

更新：2000-10-05

数値解法っていうのは、与えられた微分方程式を解析的に解く（つまりf(x)=sin(x)等の関数を求める）のではなく、グラフが描ける様に値を求めて行くという意味である。一般に微分方程式の解法にはオイラー法、Runge-Kutta法等があり、MaTXではRKF法（多分Runge-Kutta-Fehlberg法）が使われている。（これらの解法は微分方程式の数値解法としては最も一般的なものなので、数値解析の本を見れば必ず載っている）

自由度１の振動

サブテーマ「行列の固有値はなぜ「固有値」と言うのか？」予告編(笑)

バネ定数がｋ[N/m]であるバネにm[kg]の重りをつけると という微分方程式になる。これをさらに と変形する。実は、これは単振動の式なので、 という解が判っているが、ここでは、それが判っていないとして、この微分方程式をMaTXで数値解析してみる。

MaTXでは１階の（連立）微分方程式を解くコマンド（Ode)が用意されているので、 ここは状態変数を使って２階微分方程式を１階の連立微分方程式にする。

において

　とおけば、　なので、

　となって

を解く事になる。

周期が１になるように w=ω^2=(2π)^2とすると

で

Func Matrix diff_eqs(t,x,u)
  Real t;
  Matrix x,u
{
      Matrix M;
      Matrix dx;
      M=[[    0      , 1]
         [ -(2*PI)^2 , 0]];
      
      dx=M*x;
      return dx;
}

という微分方程式のMatrix関数を定義しておいて、

{T,X}=Ode45Auto(0.0,5.0,([1,0]'),diff_eqs); //微分方程式を数値的に解く関数

Y=Z(1,length(T));

Y(1,:)=X(1,:); // x(1)だけを取り出す。x(2)はx'なので、あんまり意味が無い。

mgplot(1,T,Y); // y=x(1)だけをプロット。

mgplot_cmd(1,"set grid");

mgplot_replot(1);

と流すと、、周期１の正弦波が出来ました。

そんな面倒な関数持ち出さなくても、微分の導関数が定義されているのだから、

で近似すれば、いいぢゃんという考えが湧くかと思う。
その様な近似で最も簡単なのは、 というオイラー前進差分近似法であるが、 一般的に言って（つまり僕の独断偏見では）１階微分ならともかく、２階微分とか偏微分の場合は止めたほうが無難である。
オイラー法を使うなら、後退近似とか中央近似とか修正オイラー法とかがあるのでそれらを試しながら使うべきであるが、どうせ複雑になるのなら、最初からRK法とかRKF法とかのライブラリを使った方が手っ取り早い。

それを示す為に、RKF法とオイラー法の比較を行ってみる。オイラー法では、

によりｘを求めていくので、 　の様に2階の場合には　とおいておいて

まず、１次微係数をにより求め、これからさらに、

によりｘを求めて行く。

ddxim1 = W * xim1;
dxi=dxim1 + ddxim1 * (T(i)-T(i-1));
xi1 = xi + dxi * ( T(i+1) - T(i) );
xx = xi1;
dxim1 = dxi;
xim1=xi;
xi=xi1;

実際のスクリプトとしては、こんな感じになる ( d0.mmで保存してmatxでloadすると走るはず。 )

結果は下記の通りであり、一般に（繰り返すが、独断偏見とも言う(笑)）前進近似のオイラー法を２階微分方程式の解法に使うと、解が振動する傾向にある様に思う。

まぁ、実はこれはレトリックというか、なんというか、RKF法の場合は、同じ刻み幅（Δｔ）でも計算量がオイラー法の何倍にもなっているので、それだけの計算量をこなすならば、オイラー法を使う場合は、ΔｔをRKFに比べて小さく刻まないと比較では不公平なんだが、あまり刻みすぎるとケタ落ちの危険性もあり、結局の所アルゴリズムの簡単さという長点はあるものの、例えばアセンブラとかDSPとかとにかく精度は５の次で、高速に計算するのが命ならばオイラー法が使える場面もあるとは思うが、 実際的では無いと思う。

自由度２の振動

サブテーマ「行列の固有値はなぜ「固有値」と言うのか？」(笑)

0000000

バネが１個ならば、 だったわけだが、バネを２個つけてみると、 上の重りと下の重りはそれぞれ m1,m2[kg]で、バネ定数は k1,k2[N/m]として、 初期位置からの変位をX1,X2とすると

上側(X1)の重りにおいて、
バネK1の伸びにより K1*X1の力が加わるが、さらにK2のバネからも力を受ける。
もしX1とX2が同じだけ移動すれば、K2のバネは伸び縮みしないので、力は受けないが違った場合は(X1-X2)の差だけ力を受ける。

従ってX1において

X2においては、

の微分方程式となり、

となり、結局

という２階連立微分方程式となる。

ここでは、超簡単にするために、上の重りと下の重りは同じ重さ(m=1[kg])で、 バネは２つとも同じバネ定数(K=1[N/m])とすると、

となります。詳しくは、その辺の本屋で振動関連の本の最初の方でも見てください。

この行列の固有値をMaTXで求めてみると

**************************************************************************
  MaTX Interpreter (matx)
  Windows95/NT(Visual C++) version 5.0.9
  last modified Sat Mar 18 21:34:53 JST 2000
  Copyright (C) 1989-2000, Masanobu Koga
**************************************************************************

>>K=[[-2, 1]
     [ 1,-1]]

  === [K] : ( 2, 2) ===
             ( 1)         ( 2)
  ( 1)-2.00000000E+000  1.00000000E+000
  ( 2) 1.00000000E+000 -1.00000000E+000


>>eigval(K) //固有値
  === [ans] : ( 2, 1) ===
     [   ( 1)-Real       ( 1)-Imag ]
  ( 1)-3.81966011E-001 0.00000000E+000
  ( 2)-2.61803399E+000 0.00000000E+000

>>eigvec(K) //固有ベクトル
  === [ans] : ( 2, 2) ===
  [ ( 1)-Real               ( 1)-Imag ]  [ ( 2)-Real        ( 2)-Imag ]
  ( 1)-5.25731112E-001 0.00000000E+000 -8.50650808E-001 0.00000000E+000
  ( 2)-8.50650808E-001 0.00000000E+000  5.25731112E-001 0.00000000E+000

この固有値出力のeigval(K) eigvec(K) の出力の見方は

>>K=[[-2, 1]
     [ 1,-1]]



>>eigval(K) //固有値
  === [ans] : ( 2, 1) ===
     [   ( 1)-Real       ( 1)-Imag ]
  ( 1)-3.81966011E-001 0.00000000E+000
  ( 2)-2.61803399E+000 0.00000000E+000

>>eigvec(K) //固有ベクトル
  === [ans] : ( 2, 2) ===
  [ ( 1)-Real               ( 1)-Imag ]  [ ( 2)-Real        ( 2)-Imag ]
  ( 1)-5.25731112E-001 0.00000000E+000 -8.50650808E-001 0.00000000E+000
  ( 2)-8.50650808E-001 0.00000000E+000  5.25731112E-001 0.00000000E+000

青字で示した組が、１つめの固有値、固有ベクトル。
緑字で示した組が、２つめの固有値、固有ベクトル。

を表しており、この場合、虚数部（Imag）は0なので省略すると

固有値：-0.382に対して

固有値：-2.618に対して

という、固有値と、それに相当する固有ベクトルの組が得られる。

これは固有値なので

となる。これは元の問題であったバネで言うと何を意味するか？というと、バネの式は

だったわけで、もしも、

を満たすならば

となり、

を満たし、固有値の定義である、を満たす。

逆に言うと、固有方程式を満たす、x1,x2 はという正弦波となり、これを 固有振動と称する。

この様に、固有振動（＝各部が正弦波の振動をする）を起こすx1,x2を固有ベクトルと言って、 その時の固有振動数ωがω=√ λ　から求まるので、行列の固有値を固有値と 言う様になった、という説を探しているのですが、見つからないので やっぱり僕の妄想かもしれない(笑)。

振動の様子

（MaTXスクリプトはこんなの）

固有値１の時 (sin(ω1t)の振動)		固有値２の時 (sin(ω2t)の振動）		固有値でないとき (sin(ω2t)+sin(ω2t)）という２つのモードが重なった振動をする
	初期位置： X1=-0.526 X2=-0.851		初期位置： X1=-0.851 X2= 0.526		初期位置： X1=1 X2=0

強制振動

ええっと、強制振動と減衰振動をあわせた方程式となると、
一般的に
-----------------------
|
+-+-+
o |
o ＝ ダンパー
o |
| |
+-+-+
|
●

なんかこんなやつですね。
バネだけだと なのが、速度()に比例して（比例定数c)
ブレーキ(cx')がかかるので
mx''= -kx - cx' となります。

この系に強制振動を加えるとします。一般的には　F(t)の力を加える説明が 多い様ですが、強制振動力じゃなくて、 強制振幅を与えてみます。

つまり実際に上の系の天井の固定部分を 距離±１、周波数ｗで、ぶんぶんと振ってやる感じですね。
この場合、 移動が x=sin(wt)　で行われますので、結局加えている力としては
F=mx''=m*w^2*sin(wt)の力を加えていることになるので、結局方程式は
mx'' + cx' + kx = m*w^2 * sin(wt)
簡単にするため、m=1とすると

x'' + cx' + kx = w^2 * sin(wt)

となります。

これがどんな感じになるかというと、
上式をラプラス変換すれば、
s^2*X(s) + c*s*X(s) + k*X(s) = F(s) となり、

→ X(s)= F(s) / (s^2 + c*s + k)

　となりますな。さてさて、F(s)はw^2* sin(wt) なので、
ボード線図→周波数応答＝伝達関数にsin(wt)を入力したときの出力
だから、単純にゲイン、位相曲線を求めればいいっすね。(後でゲインにw^2をかけておく）

また、実際の応答は
x'' + cx' + kx = w^2 * sin(wt)
を解けばいいですね。　

　これを連立微分方程式にするために
x1 = x （∴ x1'=x' )
x2 = x1' （∴ x2'=x1''=x'')
と置いて

(x1' = x2
(x2' = -c*x2 -k*x1 - w^2*sin(wt)

となるので
(x1') = ( 0 1)(x1) + ( 0 0)(0 )
(x2') ( -k -c)(x2) ( 0 1)(w^2*sin(w*t))

ってなマトリクス組めば、あとは微分方程式の処理系にかければよしと。

この問題は解析的には細かく解かれておりまして、

x'' + cx' + kx の系は　これを
x'' + 2ζωnx' + ωn^2*y と書き換えて

固有振動数: ωn=√(k)　となり、入力振動数：wに対してλ=w/ωn
とおいたとき、入力に対する振幅は
λ^2 / √{(1-λ^2)^2 + (2ζλ)^2 }倍とかになるそうで、
早い話しが、λ＝１の時、つまり固有振動数の入力に対しては
1/(2ζ) = √(k) / C 倍になるとか。

例えば 固有周波数を50Hzにして、この周波数の入力に対して振幅が３０倍
になるようにすると、
k=ωn^2=(2*π*50)^2 = 98696.0
C=√k/30 = (2*π*50)/30 = 10.47

(x1') = ( 0 1)(x1) + ( 0 0)(0 )
(x2') ( -k -c)(x2) ( 0 1)((2*π*50)^2*sin((2*π*50)*t))
の応答を見ると、

また、ラプラス変換で周波数応答を見ると、

となりますな。

HomePage(FrontPage)へ