力学

【力学:強制振動】2階非斉次線形微分方程式の解法をわかりやすく解説する

こんにちは(@t_kun_kamakiri)。

2020年7月頃から一ヶ月ほど、物理の質問サービスをしています。

 

カマキリ

一ヶ月で30件くらい対応しました!受験勉強のようです_(._.)_

特に多かった質問として、以下の微分方程式の解法に関する質問です。

ひとつが、↓これ。

ダンパーの振動の運動方程式

\begin{align}m\frac{d^2x}{dt^2}=-D\frac{dx}{dt}-kx\tag{1}\end{align}

※本記事は、↓こちらの内容を読んでいることが前提となります_(._.)_

前回の記事はこちら

 

もうひとつが、↓これ。

強制振動の運動方程式(←本記事の内容

\begin{align*}
m\frac{d^2x}{dt^2}=-kx-D\frac{dx}{dt}+F_{0}\sin\Omega t \tag{2}
\end{align*}

この2つの微分方程式の解法についての質問が結構ありました。

僕が学生時代にもこの手の問題というのは、理解に苦しみました。
なぜなら、大学に入ってはじめて触れる微分方程式というものだからです。

前回の記事で(1)の解法は解説済みなので、本日は(2)の解法の解説をしたいと思います。

本記事の内容はこちら
本日はこちらの微分方程式、
\begin{align*}
m\frac{d^2x}{dt^2}=-kx-D\frac{dx}{dt}+F_{0}\sin\Omega t \tag{2}
\end{align*}

の解法を理解する。
 
微分方程式を理解するだけではなく、
 
  • (2)の物理的な意味は何かを理解する
  • (2)の微分方程式の解法を理解する
 

こちらの2点を詳しく解説します。

では、さっそくやっていきましょう!

スポンサーリンク

強制振動の運動方程式

 

今回考えるのは以下のような問題です。

バネにつながれた質量\(m\)の質点が、速度に比例する抵抗(空気抵抗や摩擦)を受けている場合に加えて、振動数\(\Omega\)で周期的に力\(F_{0}\sin\Omega t\)が加わっている場合の運動方程式です。
※\(x\)は自然長からのバネの伸びです。

運動方程式

\begin{align*}
m\frac{d^2x}{dt^2}=-kx-D\frac{dx}{dt}+F_{0}\sin\Omega t \tag{2}
\end{align*}

となります。

そして、これを以下のように整理します。

  • \(\gamma=\frac{D}{2m}\)
  • \(\omega_{0}=\sqrt{\frac{k}{m}}\)
  • \(f_{0}=\frac{F_{0}}{m}\sin\Omega t\)

強制振動の運動方程式

\begin{align*}\ddot{x}+2\gamma\dot{x}+\omega_{0}^2 x=f_{0}\sin\Omega t\tag{3}\end{align*}


※\(\dot{x}=\frac{dx(t)}{dt}\)

※\(\ddot{x}=\frac{d^2x(t)}{dt^2}\)

本日、解く微分方程式は(3)です。

2階非斉次線形微分方程式とは何か

 

前回の記事では、「2階斉次線形微分方程式」についての話をしましたが、今回考える微分方程式は「2階斉次線形微分方程式」です。

以下のにその特徴をまとめておきます。

最大の特徴は、右辺が0ではなく\(f(t)\)となっていることです。
右辺が0ではないだけで、ものすごくややこしくなる・・・ように感じるのですが、この手の微分方程式の解法は前回の記事の解法の応用を効かせるだけですのでそれほど難しくありません。

カマキリ

この記事を読めば慌てずに対処できるはずです!

2階非斉次線形微分方程式の解法

 

では、これから(3)の微分方程式を解くのですが・・・・もう一度式を書いておきます。

強制振動の運動方程式

\begin{align*}\ddot{x}+2\gamma\dot{x}+\omega_{0}^2 x=f_{0}\sin\Omega t\tag{3}\end{align*}

  • \(\gamma=\frac{D}{2m}\)
  • \(\omega_{0}=\sqrt{\frac{k}{m}}\)
  • \(f_{0}=\frac{F_{0}}{m}\sin\Omega t\)

解き方は定石というのがあります。

以下の「解法のフロー」を用意しました。

このフローにしたがって(何だったら何も考えずに)、計算を進めていけば良いです。
微分方程式を十分に習っていない状態であれば、機械的にこれにしたがってやりましょう!!

では、順番にやっていきましょう!

\(f_{0}=0\)の場合の一般解を求める

 

\(f_{0}=0\)についての微分方程式↓

\begin{align*}\ddot{x}+2\gamma\dot{x}+\omega_{0}^2 x=0\tag{4}\end{align*}

この微分方程式の解は、前回の記事で求めましたので結果だけ記載しておきます。

\begin{align*}
x(t)=A e^{(-\gamma+\sqrt{\gamma^2-\omega_{0}^2})t}+B e^{(-\gamma-\sqrt{\gamma^2-\omega_{0}^2})t}\tag{4}
\end{align*}

\(A,B\)が積分定数です。
(3)は2階の微分を含む微分方程式ですので、(3)を2回積分して得られる解\(x(t)\)は必ず2つの積分定数が含まれています。

であれば、(4)で求めた一般解(7)を利用するのが上手い方法である・・・と考えるのが自然だと思います。

では、(3)の一般解を以下の形であると仮定します。

(3)の一般解

\begin{align*}
x(t)&=A e^{(-\gamma+\sqrt{\gamma^2-\omega_{0}^2})t}+B e^{(-\gamma-\sqrt{\gamma^2-\omega_{0}^2})t}+C\sin\Omega t+D\cos\Omega t \\
&=x_{1}(t)+x_{2}(t)\tag{5}
\end{align*}

※(4)の一般解の項を\(x_{1}(t)=A e^{(-\gamma+\sqrt{\gamma^2-\omega_{0}^2})t}+B e^{(-\gamma-\sqrt{\gamma^2-\omega_{0}^2})t}\)でまとめておきました。
※(3)の特殊解の項を\(x_{2}(t)=C\sin\Omega t+D\cos\Omega t\)でまとめておきました。

(3)の左辺は「(5)の\(x_{1}(t)\)で0になる」として、(3)の右辺を出すためには同じ振動数\(\Omega\)を使った周期関数を持ってくるのがセンスがいいんじゃないかって思いますよね。
しかも、積分定数はすでに2つ使っているので、これ以上積分定数を増やす必要もありません。

なので、注力すべきは「(5)が(3)を満たしてくれるように\(C,D\)を決定する」ということになります。

強制振動の特殊解\(x_{2}(t)\)を求める

「(5)が(3)を満たすように\(C\)を決定すればよい」という言いましたが、少し急ぎすぎたかもしれません。
 
まず、(5)の一般解は「\(x(t)=x_{1}(t)+x_{2}(t)\)」と書きました。
  • \(x_{1}(t)\)は\(f_{0}=0\)での一般解
  • \(x_{2}(t)\)は(5)の特殊解(積分定数を含まない)」

ということを意識する必要があります。

だから、そもそも\(x_{1}(t)\)は以下を満たしているという事を頭に置いておかなくてはいけません。
 
\begin{align*}\ddot{x_{1}}+2\gamma\dot{x_{1}}+\omega_{0}^2 x_{1}=0\tag{6}\end{align*}
 
しつこいくらい強調しているのは、自分自身がはじめて強制振動を習ったときに十分理解していなかったので、全く何をやっているのか理解ができなかったからです。
 
では、\(x(t)=x_{1}(t)+x_{2}(t)\)を(5)に代入してみましょう。
 
\begin{align*}\big(\ddot{x_{1}}+2\gamma\dot{x_{1}}+\omega_{0}^2 x_{1}\big)+\big(\ddot{x_{2}}+2\gamma\dot{x_{2}}+\omega_{0}^2 x_{2}\big)=f_{0}\sin\Omega t\tag{7}\end{align*}
これより(6)を用いると、\(x_{1}\)が消えてくれます。
 
\begin{align*}
\ddot{x_{2}}+2\gamma\dot{x_{2}}+\omega_{0}^2 x_{2}=f_{0}\sin\Omega t\tag{8}
\end{align*}
となります。
 
なんだ全然さっきと微分方程式が変わってないじゃないかと思うかもしれませんが、ちょっと違います。
(8)の\(x_{2}(t)\)は、積分定数が必要ない特殊解で良いという事です。つまり、とりあえず(8)を満たしてくれるような解を探してこれば良いという事です。
 
となると、先ほども言ったように右辺が振動数\(\Omega\)の周期関数であるなら、特殊解も\(\Omega\)を持つ周期関数にしておくのがセンスが良いのではないかと考えますよね。
 
ならば、(8)の特殊解を、
\begin{align*}
x_{2}=C\sin\Omega t+D\cos\Omega t\tag{9}
\end{align*}
とおいてみてはどうかとなります。
 
こうすることで、ここでも目的がはっきりしてきます。
 
【ここでの目的】
(8)を満たすような特殊解を「\(x_{2}=C\sin\Omega t+D\cos\Omega t\)」と仮定して、(8)を満たしてくれるように\(C\)を定める。

 

目的は一言で言うと「\(C\)」を決めることですね。

では、(9)を(8)に代入してみましょう。

\begin{align*}
x_{2}=C\sin\Omega t+D\cos\Omega t\tag{9}
\end{align*}

(9)を微分したものを用意しておいて・・・

  • \(\frac{dx_{2}}{dt}=C\Omega\cos\Omega t-D\omega\sin\omega t\)
  • \(\frac{d^2x_{2}}{dt^2}=-C\omega^2\sin\Omega t-D\Omega^2\cos\Omega t\)

なので、(8)に代入すると、

\begin{align*}
&\bigg(-C\Omega^2\sin\Omega t-D\Omega^2\cos\Omega t\bigg)\\
&+\gamma\bigg(C\Omega\cos\Omega-D\Omega\sin\Omega t\bigg)\\
&+\Omega_{0}^2\bigg(C\sin\Omega t+D\cos\Omega t\bigg)=f_{0}\sin\Omega t\tag{10}
\end{align*}

これを\(\sin,\cos\)にまとめて、

\begin{align*}
&\bigg(-C\Omega^2-\gamma D\Omega+C\omega_{0}^2\bigg)\sin\Omega t\\
&+\bigg(-D\Omega^2+\gamma C\Omega +\omega_{0}^2 D\bigg)\cos\omega t \\
&=f_{0}\sin\Omega t\tag{11}
\end{align*}

これより、

  • \(C=\frac{(-\omega^2+\omega_{0}^2)f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)
  • \(D=\frac{-\gamma\omega f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)

\(C,D\)を計算する方法として中には、\(C\)と\(D\)の連立方程式を愚直に計算した人もいることでしょう。
その計算方法は個々人の方法で良いかと思いますが、今回の場合は行列計算を考える方が手っ取り早いですし、ミスが少ないと思いますので少し紹介をしておきます。

ず、\(\sin\)と\(\cos\)の係数が両辺で等しいとし、\(C\)と\(D\)については以下のように連立方程式が解きやすいように整理します。

\begin{align*}
\left\{\begin{matrix}
\big(-\Omega^2+\omega_{0}^{2}\big)C-\gamma\Omega D=f_{0}\\
\gamma\Omega C+\big(-\Omega^2+\omega_{0}^{2}\big)D=0
\end{matrix}\right.
\end{align*}

これを行列になおします。

\begin{align*}
\begin{pmatrix}
-\Omega^2+\omega_{0}^{2} & -\gamma\Omega\\
\gamma\Omega & -\Omega^2+\omega_{0}^{2}
\end{pmatrix}
\begin{pmatrix}
C \\
D
\end{pmatrix}=
\begin{pmatrix}
f_{0} \\
0
\end{pmatrix}
\end{align*}

2行2列の行列計算は以下のように簡単に行えます。

\begin{align*}
\begin{pmatrix}
C \\
D
\end{pmatrix}=
\frac{1}{\big(-\Omega^2+\omega_{0}^2\big)^2+\big(\gamma\Omega\big)^2}
\begin{pmatrix}
-\Omega^2+\omega_{0}^{2} & \gamma\Omega\\
-\gamma\Omega & -\Omega^2+\omega_{0}^{2}
\end{pmatrix}
\begin{pmatrix}
f_{0} \\
0
\end{pmatrix}
\end{align*}

これより、\(C,D\)が求まります。

  • \(C=\frac{(-\omega^2+\omega_{0}^2)f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)
  • \(D=\frac{-\gamma\omega f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)
カマキリ

この方法なら簡単にC,Dが求めることができます。

2階非斉次線形微分方程式の一般解は?

以上のような流れを確認しつつ、(3)の一般解を求めることができました。

(3)の一般解

\begin{align*}
x(t)=A e^{(-\gamma+\sqrt{\gamma^2-\omega_{0}^2})t}+B e^{(-\gamma-\sqrt{\gamma^2-\omega_{0}^2})t}+C\sin\Omega t+D\cos\Omega t \tag{12}
\end{align*}

※(4)の一般解の項を\(x_{1}(t)\)でまとめておきました。
※(3)の特殊解の項を\(x_{2}(t)\)でまとめておきました。

  • \(C=\frac{(-\omega^2+\omega_{0}^2)f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)
  • \(D=\frac{-\gamma\omega f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)

\(C,D\)は代入せずに別で記述しておきました。

代入して書くと、

\begin{align*}
x(t)&=A e^{(-\gamma+\sqrt{\gamma^2-\omega_{0}^2})t}+B e^{(-\gamma-\sqrt{\gamma^2-\omega_{0}^2})t}\\
&+\frac{(-\omega^2+\omega_{0}^2)f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\sin\Omega t\\
&+\frac{-\gamma\omega f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\cos\Omega t \tag{13}
\end{align*}

↑めちゃくちゃごちゃごちゃしてますね。

(13)が本当に(3)の一般解になっているのか、(3)の左辺に代入して確認してみてください。
ちゃんと満たしていることが確認できます。

まとめ

 

今回は「2階非斉次線形微分方程式」の解法について解説をしました。

まず、以下のフローを今一度確認しておきましょう!

 

(3)の一般解

\begin{align*}
x(t)=A e^{(-\gamma+\sqrt{\gamma^2-\omega_{0}^2})t}+B e^{(-\gamma-\sqrt{\gamma^2-\omega_{0}^2})t}+C\sin\Omega t+D\cos\Omega t \tag{12}
\end{align*}

※(4)の一般解の項を\(x_{1}(t)\)でまとめておきました。
※(3)の特殊解の項を\(x_{2}(t)\)でまとめておきました。

  • \(C=\frac{(-\omega^2+\omega_{0}^2)f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)
  • \(D=\frac{-\gamma\omega f_{0}}{(-\omega^2+\omega_{0}^2)^2+(\gamma\omega)^2}\)

内容は以上となります。

この「2階非斉次線形微分方程式」は大学一年生の力学で一番初めに習う微分方程式だと思います。

大学院試験やレポート問題でも最頻出のひとつですので必ずマスターしておきましょう(^^)/

おまけ:空気抵抗は本当に速度に比例するのか

 

本記事の内容で、速度に比例する項を「空気抵抗」と呼びましたが、実は空気抵抗は速度に比例するものではありません。

空気抵抗がどのように速度に依存するかは、正確には流れによるという事です。

「空気抵抗」に関する記事はこちら

ここでは、簡単にまとめておきましょう!

  • 層流状態➡抵抗は速度に比例
  • 乱流状態➡抵抗は速度の2乗に比例

 

おすすめの参考書

 

大学一年生がはじめて力学を学ぶときに、高校物理と大学物理の難易度のギャップに戸惑う人も多いでしょう。
ここで、紹介するのは高校物理と大学物理を埋めてくれるわかりやすい参考書ですので、どうしれもつまづいて先に進めないという方は是非読んでみてください。

力学キャンパス・ゼミ 改訂6

力学キャンパス・ゼミ 改訂6

馬場 敬之
2,783円(11/27 19:56時点)
発売日: 2020/02/05
Amazonの情報を掲載しています
演習 力学キャンパス・ゼミ

演習 力学キャンパス・ゼミ

馬場 敬之, 高杉 豊
2,299円(11/28 08:22時点)
発売日: 2018/12/13
Amazonの情報を掲載しています

高校数学と物理を教えている予備校講師陣のわかりやすい参考書と演習書を読めば大学物理の触りはOKです。

難易度として同じくらいの以下の参考書もお勧めです。

力学 (講談社基礎物理学シリーズ)

力学 (講談社基礎物理学シリーズ)

副島雄児, 杉山忠男
2,750円(11/28 08:22時点)
発売日: 2009/09/25
Amazonの情報を掲載しています
カマキリ

では(^^)/

【プロフィール】

カマキリ
(^^)

大学の専攻は物性理論で、Fortranを使って数値計算をしていました。
CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。

プロフィール記事はこちら

 

大学学部レベルの物理の解説をします 大学初学者で物理にお困りの方にわかりやすく解説します。

このブログでは主に大学以上の物理を勉強して記事にわかりやすくまとめていきます。

  • ・解析力学
  • ・流体力学
  • ・熱力学
  • ・量子統計
  • ・CAE解析(流体解析)
  • note
    noteで内容は主に「プログラミング言語」の勉強の進捗を日々書いています。また、「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。
  • youtube
    youtubeではオープンソースの流体解析、構造解析、1DCAEの操作方法などを動画にしています。
    (音声はありません_(._.)_)
  • Qiita
    Qiitaではプログラミング言語の基本的な内容をまとめています。