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

作成:2000-04-20

更新:2000-10-05


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


自由度1の振動

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

 

バネ定数が[N/m]であるバネにm[kg]の重りをつけると

という微分方程式になる。これをさらに

と変形する。実は、これは単振動の式なので、

という解が判っているが、ここでは、それが判っていないとして、この微分方程式をMaTXで数値解析してみる。

 

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

において

 とおけば、 なので、

 となって


を解く事になる。

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

MaTX Interpreter (matx)
Windows95/NT(Visual C++) version 5.0.5
last modified Fri Jan 28 18:28:19 JST 2000
Copyright (C) 1989-2000, Masanobu Koga

Send bugs and comments to matx@matx.org
Type 'quit' to exit, 'help' for functions, 'demo' for demonstration.

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);

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

 

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


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

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

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

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

によりxを求めて行く。

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すると走るはず。 )

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

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


自由度2の振動

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

 

0000000

 

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

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

従ってX1において

X2においては、

の微分方程式となり、

となり、結局

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

ここでは、超簡単にするために、上の重りと下の重りは同じ重さ(m=1[kg])で、 バネは2つとも同じバネ定数(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

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

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

固有値:-0.382に対してclip_a7.png

固有値:-2.618に対して

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

これは固有値なので

clip_a9.png

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


clip_a11.png

だったわけで、もしも、

clip_a12.png

を満たすならば
clip_a13.png

となり、

clip_a14.png

 

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

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

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

振動の様子

(MaTXスクリプトはこんなの)

固有値1の時
(sin(ω1t)の振動)
固有値2の時
(sin(ω2t)の振動)

固有値でないとき

(sin(ω2t)+sin(ω2t))という2つのモードが重なった振動をする

 

初期位置:
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)の力を加える説明が 多い様ですが、強制振動力じゃなくて、 強制振幅を与えてみます。

つまり実際に上の系の天井の固定部分を 距離±1、周波数wで、ぶんぶんと振ってやる感じですね。
この場合、 移動が 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の時、つまり固有振動数の入力に対しては
1/(2ζ) = √(k) / C 倍になるとか。

例えば 固有周波数を50Hzにして、この周波数の入力に対して振幅が30倍
になるようにすると、
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)へ