連立常微分方程式の解き方が分かりません

問題は
X＝x1－x2・・・(1)
m1＊x1"=Tr-m1*g・・・(2)
m2*x2"+c2*x'+k*x2=-Tr+c2*X'+k2*X・・・(3)
X'=G*(Tr-To)・・・(4)
という連立微分方程式です。
変数がx1,x2,X,Trでそれぞれ時間tの関数です。
私は普段数値計算にMatlabを用いておりますが、数学、Matlabの知識ともに未熟でして、解くことができません。
私はそれぞれの式を差分化して強引に計算したのですが、上手くいきません。どなたかこの連立微分方程式の解き方を教えていただけないでしょうか？
よろしくお願いいたします。

No.3 回答者： inara1 回答日時： 2008/06/22 10:06

＞解析手法はどういったものか教えていただけないでしょうか？

特に変わった方法ではありません（最後に書いておきます）。

ご質問の連立微分方程式の (3) ですが、両辺に c2*X' があるので、この項は消えてしまいます。そうするとこれは空気抵抗は働かないモデルになってしまいます。また、方程式には小文字の x と大文字の X が混ざっていますが、これらは全部同じ（台車の位置）ですね。

これは以前の質問（http://okwave.jp/qa4110689.html）の連立微分方程式ですが、このモデルはおかしいと思います。
台車を強制的に運動させているのなら、台車にかかる空気抵抗も台車に加える力（T0）も関係ないことになりませんか？ 釣竿部分の取扱いも変です。それと、釣竿と錘との間が単なる糸なのか曲がらない線なのかで錘の運動が異なります（台車が減速させたとき）。これは振動工学の課題でしょうか？釣竿の質量や釣竿に働く空気抵抗を考慮しないといけないのか、釣竿と錘との間のところの扱いも確認してください。

【方程式の解き方】
問題の連立微分方程式で、x を X に統一して、式(3)の両辺にある c2 の項を消してしまうと

　　　X = x1 - x2 --- (1)
　　　m1*x1'' = Tr - m1*g --- (2)
　　　m2*x2'' + k*x2 = -Tr + k2*X --- (3)
　　　X' = G*( Tr - To ) --- (4)

となります。式(1) から
　　　x2 = x1 - X --- (5)
これを(3)に代入すると
　　　m2*( x1'' - X'' ) + k*( x1 - X ) = -Tr + k2*X --- (6)
式(2)から
　　　Tr = m1*x1'' + m1*g --- (7)
これを式(6)に代入すると
　　　m2*( x1'' - X'' ) + k*( x1 - X ) = -m1*x1'' - m1*g + k2*X --- (8)
一方、式(7)を式(4)に代入すると
　　　X' = G*( m1*x1'' + m1*g - To )
両辺を時間 t で積分すれば
　　　X = G*{ m1*x1' + ( m1*g - To )*t + C4 } （ C4 は積分定数） --- (9)
式(9)を式(8)に代入すれば、以下のように x1 だけの微分方程式となります。
　　　m2*( x1'' - G*m1*x1''' ) + k*[ x1 - G*{ m1*x1' + ( m1*g - To )*t + C4 } ] = -m1*x1'' - m1*g + k2*G*{ m1*x1' + ( m1*g - To )*t + C4 } --- (10)
この解は
　　　x1 = A + B*t + C*exp( D*t ) --- (11)
となります。式(11)を式(10)に代入して、その結果が t によらず成り立つように A, B, C, D を求めます。A と B はすぐに出てきますが、D は3次方程式の解になるので、その解（一般に複素数）を D1, D2, D3 としたとき、一般解は
　　　x1 = A + B*t + C1*exp( D1*t ) + C2*exp( D2*t )+ C3*exp( D3*t )
となります。式(10)はMatlabでも解けるはずです。x1 が分かったら、X と x2 と Tr はそれぞれ、式(9)、式(5)、式(7)から求められます。式(9)の積分定数 C4 は、t = 0 のとき X = 0 なら C4 = 0 となります。X や x1、x2 は絶対位置でなく、t = 0 （全部静止？）からの変位と考えたほうが良いです。D1～D4 が実数なのか純虚数なのか、それ以外の複素数なのかで場合分けしないと、x1 の具体的な形は決まりません。

この回答へのお礼

ご回答まことにありがとうございます。
ご指摘された点についてですが、(3)式については
m2*x2"+c*x2'+k*x2=-Tr+c*X'+k*X・・・(3)
であり、私の記述ミスです。
そして、c、kはそれぞれ竿をばねと見立てた時のばね定数、減衰係数
として考えております。また、Trは釣竿が発生させる力ですが、Toは台車に加える力ではなく、釣竿に発生させたい力となります。そして、台車は速度制御ができ、(4)で表されるように、ToとTrが等しくなった時に台車の速度を0とすることで、任意の力Toを釣竿に発生させる事を目的としております。

また、釣竿と錘との間の糸は曲がらない線で、台車が減速させたとき
も糸は決してたるまないものとして考えております。

これは振動工学の問題と言うよりも、私が現在行っている研究に関する問題であり、実際に釣竿（片持ち梁）とリニアアクチュエータを用いて、任意の力を発生させる機構を作成したものを、モデル化しようと考えているものです。

自分の説明が非常に足りなかったので大変申し訳なく思っております。

また、回答していただいた方程式の解き方については、これから参考にしながら自分でやってみようと考えております。

お礼日時：2008/06/23 17:15

No.2 回答者： inara1 回答日時： 2008/06/21 16:51

これは x1(t) に関する3階微分方程式になりますが解析的に解けます（Mapleという数式処理ソフトで解けました）。

x1 は以下の形になりましたが、大変複雑な式です。
　　　x1(t) = A + B*t + C1*exp( D1*t ) + C2*exp( D2*t ) + C3*exp( D3*t )
A と B は実数ですが、C1～C3 と D1～D3 は複素数です。x1 が分かれば、他の関数は以下の関係で求められます。
　　　X = G*{ m1*x1' + ( m1*g - T0 )*t + C4 }
　　　x2 = x1 - X
　　　tr = X'/G + T0
ところで、質問文の g と G は同じですか？

この回答への補足

ご回答まことにありがとうございます。
Mapleというソフトは使ったことがないのですが、inara1様が解いてくださった解析手法はどういったものか教えていただけないでしょうか？
また、複素数が含まれているとのことで、大変結果がイメージしづらいのですが、Trはtが増加するにつれて振動しながらToに近づいていくような感じになりますでしょうか？
ご質問ばかりさせていただき大変申し訳ございません。

また、ご指摘にありますgとGは別のものになります。

補足日時：2008/06/22 04:40

No.1 回答者： rabbit_cat 回答日時： 2008/06/19 00:52

とりあえず

x1' = x11 ・・・(5)
x11' = x12 ・・・(6)
X' = X1 ・・・(7)
みたいな式をつけ加えて、変数が x1,x2,X,Tr,x11,x12,X1 の７元連立線形１次微分方程式にします。
で、X↑ = (x1,x2,X,Tr,x11,x12,X1) という縦ベクトルを考えると、
A*X↑' = B*X↑
という形に整理できるはずです。（A,Bは７×７の正方行列）
そうしたら、
X↑' = A^(-1)*B*X↑
で、
X↑ = C*exp(t*A^(-1)*B)
ただし、Cは７×７の正方行列（積分定数）、
exp(t＊A^(-1)*B) は行列のexponential、て形で解けます。

matlabで逆行列を求めるのは、inv()
http://dl.cybernet.co.jp/matlab/support/manual/r …

行列のexponentialを求めるのは、expm()
http://dl.cybernet.co.jp/matlab/support/manual/r …

です。