強制減衰振動モデルをFourier変換で解く
from ✅@2023-06-23T08:50D90 AM4-2023F-10
強制減衰振動モデルをFourier変換で解く
2階線型非斉次微分方程式
強制振動&減衰振動式を、Fourier変換で解く
$ m\ddot x+c\dot x+kx=f(t)
$ k:ばね定数
$ c:減衰係数
バネとダンパー(ダッシュポット)が繋がった台車のモデル
https://kakeru.app/f089542c8fd98e4cad0937d9c8b64c46 https://i.kakeru.app/f089542c8fd98e4cad0937d9c8b64c46.svg
$ mで割って、monicな式にする
$ \ddot x+2\omega_0h\dot x+{\omega_0}^2x=F(t)ー①
$ h:=\frac1{2\omega_0}\frac cm
$ \omega_0は正確には非減衰固有角振動数と呼ぶ
臨界減衰比$ hが$ h=0のときの固有角振動数だから、「非減衰」固有角振動数と呼ぶ
$ h=0 \text{ then }\ddot x+{\omega_0}^2x=0
$ F(t):=\frac{f(t)}{m}
$ \omega_0と$ hは、耐震設計、振動工学で重要
両辺に$ {\cal F}を作用させる
$ {\cal F}(\ddot x)+2\omega_0h{\cal F}(\dot x)+{\omega_0}^2{\cal F}(x)={\cal F}(F)
$ \iff -\omega^2{\cal F}(x)(\omega)+2i\omega\omega_0h{\cal F}(x)(\omega)+{\omega_0}^2{\cal F}(x)(\omega)={\cal F}(F)(\omega)
$ \iff -\omega^2{\cal F}(x)(\omega)+2i\omega\omega_0h{\cal F}(x)(\omega)+{\omega_0}^2{\cal F}(x)(\omega)={\cal F}(F)(\omega)
$ \iff (-\omega^2+2i\omega\omega_0h+{\omega_0}^2){\cal F}(x)(\omega)={\cal F}(F)(\omega)
Fourier変換を作用させると、代数方程式に変換されるのか!takker.icon
演算子法やラプラス変換みたい
$ \iff {\cal F}(x)(\omega)=\frac{{\cal F}(F)(\omega)}{-(\omega -i\omega_0h)^2+{\omega_0}^2(1-h^2)}
$ i\omegaをひとかたまりにして書き直すと
$ \iff ((i\omega)^2+2h\omega_0(i\omega)+{\omega_0}^2){\cal F}(x)(\omega)={\cal F}(F)(\omega)ー②
$ \iff {\cal F}(x)(\omega)=\frac{{\cal F}(F)(\omega)}{(i\omega +\omega_0h)^2+{\omega_0}^2(1-h^2)}
$ =:{\cal F}(G)(\omega){\cal F}(f)(\omega)
ここで、伝達函数$ {\cal F}(G)(\omega)を定義した
$ {\cal F}(G)(\omega):=\frac1m\frac1{(i\omega +\omega_0h)^2+{\omega_0}^2(1-h^2)}
$ = -\frac1m\frac1{\omega^2-{\omega_0}^2-2h\omega_0i\omega}
分母を$ (\omega-\omega_\alpha)(\omega-\omega_\beta)の形にする
単に$ \omegaについての2次方程式を解けばいいだけ
$ \omega=ih\omega_0\pm\sqrt{-h^2{\omega_0}^2+{\omega_0}^2}
$ = \left(ih\pm\sqrt{1-h^2}\right)\omega_0
平方完成をつかうとこう
✅@2023-06-23T08:50D90 AM4-2023F-10#6494e5c81280f0000023b335からの続き
$ \iff\omega=ih\omega_0\pm\sqrt{1-h^2}\omega_0
まあ同じことである
$ \omega_\alpha=ih\omega_0+\sqrt{1-h^2}\omega_0、$ \omega_\beta=ih\omega_0-\sqrt{1-h^2}\omega_0としておく
$ \omega_\alphaは複素固有角振動数
$ \omega_\beta=-{\omega_\alpha}^*
以上より
$ {\cal F}(x)(\omega)={\cal F}(G)(\omega){\cal F}(f)(\omega)
$ {\cal F}(G):\omega\mapsto-\frac1m\frac1{\omega-\omega_\alpha}\frac1{\omega-\omega_\beta}
となる。
臨界減衰比と振動モードとの関係
あとは、これを逆変換で戻せばいいだけ
$ x(t)=\frac1m{\cal F}^{-1}\left(\frac1{\omega-\omega_\alpha}\frac1{\omega-\omega_\beta}{\cal F}(f)(\omega)\right)(t)
$ {\cal F}^{-1}(xy)(t)の演算法則を求める必要があるのかtakker.icon
$ {\cal F}(xy)=\int_\R x(t)y(t)e^{-i\omega t}\mathrm dt
$ = -\frac{1}{i\omega}\int_{t\in\R}x(t)y(t)\mathrm de^{-i\omega t}
$ =0+\frac{1}{i\omega}\int_{t\in\R}e^{-i\omega t}(x\mathrm dy+y\mathrm dx)
うーん、だめだな。分解できない
伝達函数$ {\cal F}(G)の性質を調べる
$ {\cal F}(G)は$ C:=\{\omega_\alpha,\omega_\beta\}で特異点を持ち、$ \Complex\setminus Cで正則である
$ h=0のとき$ C=\{\pm\omega_0\}だから、伝達函数の特異点は固有振動数という事実がわかる
今回は時間がないので、$ {\cal F}(f)(\omega)=1としておく
$ f=\frac1{2\pi}\int_\R e^{i\omega t}\mathrm d\omega
これはDiracのdelta函数らしい
次週やる
$ {\cal F}^{-1}{\cal F}(G)を求めるために、Cauchyの積分定理を用いる
正則函数であることとCauchy-Riemannの方程式を満たすことは同値である
$ fが正則\iff \frac{\partial\Re f}{\partial\Re z}=\frac{\partial\Im f}{\partial\Im z}\land\frac{\partial\Re f}{\partial\Im z}=-\frac{\partial\Im f}{\partial\Re z}
ちょい違ったかもtakker.icon
例
$ r>0、$ D:=\{z\in\Complex||z|\le r\}のとき、$ I=\int_{\partial D}z\mathrm dzを求める
$ I= \int_{re^{i\theta}\in\partial D}re^{i\theta}\mathrm dre^{i\theta}
$ = i\int_{0\le\theta<2\pi}e^{2i\theta}\mathrm d\theta
$ = \frac12\int_{0\le\theta<2\pi}\mathrm de^{2i\theta}
$ =0
とくに、正則函数$ fにて次が成立する
$ \oint_{\partial D}\frac{f(z)}{z-\alpha}\mathrm dz=2\pi if(\alpha)
式の意味
$ C:=\{z\in\Complex||z-\alpha|=\varepsilon\}を考える
$ C上の微小線素$ \mathrm dzは、パラメタ$ 0\le\theta<2\piを使って以下のように表せる
$ \mathrm dz=\varepsilon e^{i\theta}i\mathrm d\theta
$ \partial D=Cとして積分すると
$ \oint_{\partial D}\frac{f(z)}{z-\alpha}\mathrm dz=\int_0^{2\pi}\frac{f(\alpha+\varepsilon e^{i\theta})}{\varepsilon e^{i\theta}}\varepsilon e^{i\theta}i\mathrm d\theta
$ =i\int_0^{2\pi}f(\alpha+\varepsilon e^{i\theta})\mathrm d\theta
領域を特異点外にずらした
$ =i\int_0^{2\pi}f(\alpha+\varepsilon e^{i\theta})\mathrm d\theta
$ \to2\pi if(\alpha)\quad(\varepsilon\to0)
???ここなにしたtakker.icon
from ✅@2023-06-30T08:50D90 AM4-2023F-11
おさらい
2階線型非斉次微分方程式$ m\ddot x+c\dot x+kx=f(t)をFourier変換で解いた結果、
$ x(t)={\cal F}^{-1}(\omega\mapsto{\cal F}(G)(\omega){\cal F}(f)(\omega))
$ {\cal F}(G):\omega\mapsto-\frac1m\frac1{\omega-\omega_\alpha}\frac1{\omega-\omega_\beta}
伝達函数と呼ぶ
伝達函数をFourier逆変換した以下の式を求めるのが今回のテーマの一つ
$ G:t\mapsto\frac1{2\pi}\int_\R{\cal F}(G)(\omega)e^{i\omega t}\mathrm d\omega
これを求めるために、Cauchyの積分定理とJordanの補題を使う
材料
Jordanの補題より$ \forall t>0;\int_{M_+(R)}F(\omega)e^{i\omega t}\mathrm d\omega\to0\quad(R\to\infty)
$ M_+(R):=\{\omega\in\Complex|\exist\theta\in[0,\pi];\omega=Re^{i\theta}\} とした
同様に$ \forall t<0;\int_{M_-(R)}F(\omega)e^{i\omega t}\mathrm d\omega\to0\quad(R\to\infty)が成り立つ
$ M_-(R):=\{\omega\in\Complex|\exist\theta\in[0,\pi];\omega=Re^{-i\theta}\} とした
Cauchyの積分定理の拡張(Cauchyの留数定理の系でもある)
$ \oint_{\partial D}\frac{f(z)}{z-\alpha}\mathrm dz=2\pi if(\alpha)
$ \alpha\in Dとする
以上より、伝達函数$ Gを求める
$ {\cal F}(G)(\omega)=-\frac1m\frac1{\omega-\omega_\alpha}\frac1{\omega-\omega_\beta}
$ = -\frac1m\frac{1}{\omega_\alpha-\omega_\beta}\left(\frac{1}{\omega-\omega_\alpha}-\frac{1}{\omega-\omega_\beta}\right)
部分分数分解で特異点$ \omega_\alpha,\omega_\betaを線型的に分離することで、Cauchyの積分定理の拡張を適用できるようにした
$ = -\frac1m\frac{1}{2\sqrt{1-h^2}\omega_0}\left(\frac{1}{\omega-\omega_\alpha}-\frac{1}{\omega-\omega_\beta}\right)
$ \implies G(t)=\frac1{2\pi}\int_\R{\cal F}(G)(\omega)e^{i\omega t}\mathrm d\omega
$ = -\frac1m\frac{1}{2\sqrt{1-h^2}\omega_0}\left(\frac1{2\pi}\int_\R\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega-\frac1{2\pi}\int_\R\frac{e^{i\omega t}}{\omega-\omega_\beta}\mathrm d\omega\right)
$ =:-\frac1m\frac{1}{2\sqrt{1-h^2}\omega_0}(J_1+J_2)
$ J_1,J_2をCauchyの積分定理の拡張より求める
$ J_1=\frac1{2\pi}\int_\R\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega
$ = \lim_{R\to\infty}\left(\frac1{2\pi}\oint_{M_+(R)\cup[-R,R]}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega-\frac1{2\pi}\int_{M_+(R)}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega\right)\llbracket t>0\rrbracket
$ +\lim_{R\to\infty}\frac1{2\pi}\int_{-R}^R\frac1{\omega-\omega_\alpha}\mathrm d\omega\llbracket t=0\rrbracket
$ +\lim_{R\to\infty}\left(\frac1{2\pi}\oint_{M_-(R)\cup[-R,R]}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega-\frac1{2\pi}\int_{M_-(R)}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega\right)\llbracket t<0\rrbracket
積分範囲を、周回路から半周を抜いたものに置き換える
これで各第一項にCauchyの積分定理を適用できる状態にできた
$ = \lim_{R\to\infty}\left(\frac{2\pi ie^{i\omega_\alpha t}}{2\pi}-\frac1{2\pi}\int_{M_+(R)}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega\right)\llbracket t>0\rrbracket+\lim_{R\to\infty}\frac1{2\pi}(\ln|R-\omega_\alpha|-\ln|-R-\omega_\alpha|)\llbracket t=0\rrbracket+\lim_{R\to\infty}\left(0-\frac1{2\pi}\int_{M_-(R)}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega\right)\llbracket t<0\rrbracket
$ \becauseCauchyの積分定理
$ \omega_\alpha は$ M_+(R)\cup[-R,R] 内にいるので、留数定理から定数項が残る
$ \omega_\alpha は$ M_-(R)\cup[-R,R] 内にいないので、0になる
$ = ie^{i\omega_\alpha t}\llbracket t>0\rrbracket+\lim_{R\to\infty}\frac1{2\pi}\int_{M_+(R)}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega\llbracket t>0\rrbracket+\lim_{R\to\infty}\frac1{2\pi}\ln\frac{|R-\omega_\alpha|}{|R+\omega_\alpha|}+\lim_{R\to\infty}\frac1{2\pi}\int_{M_-(R)}\frac{e^{i\omega t}}{\omega-\omega_\alpha}\mathrm d\omega\llbracket t<0\rrbracket
$ = ie^{i\omega_\alpha t}\llbracket t>0\rrbracket+0+0+0
$ \becauseJordanの補題
$ J_2=ie^{i\omega_\beta t}\llbracket t>0\rrbracket
導出はほぼ同じなので略
$ \omega_\alpha,\omega_\betaはともに虚軸正の位置にある
以上より、
$ G(t)=-\frac im\frac{1}{2\sqrt{1-h^2}\omega_0}(e^{i\omega_\alpha t}-e^{i\omega_\beta t})\llbracket t>0\rrbracket
$ =\frac1m\frac{1}{\sqrt{1-h^2}\omega_0}e^{-h\omega_0t}\frac1{2i}(e^{i\sqrt{1-h^2}\omega_0t}-e^{-i\sqrt{1-h^2}\omega_0t})\llbracket t>0\rrbracket
$ \underline{=\frac1m\frac{1}{\sqrt{1-h^2}\omega_0}e^{-h\omega_0t}\sin\sqrt{1-h^2}\omega_0t\llbracket t>0\rrbracket\quad}_\blacksquare
$ GをGreen函数と呼ぶ
Ramp函数と三角波を組み合わせたような函数
https://gyazo.com/4d723973ba7a9da18f79dd17f26b7a7c
$ G(t)が表す物理現象はなにか?
前回のモデルを思い出す
https://kakeru.app/f089542c8fd98e4cad0937d9c8b64c46 https://i.kakeru.app/f089542c8fd98e4cad0937d9c8b64c46.svg
このモデルは$ m\ddot x+c\dot x+k x=f(t)で表せる
Fourier変換すると、$ {\cal F}(x)(\omega)={\cal F}(G)(\omega){\cal F}(f)(\omega)になった
もし$ {\cal F}(f):\omega\mapsto1の定数函数だったら、
$ \boxed{x(t)=G(t)}
つまり、外力のFourier変換が1のときの台車の動きが$ G(t)となる
$ t=0から急に減衰振動し始める動き
一般化するには、$ f(t)を求める必要がある
$ {\cal F}(f):\omega\mapsto1のとき、$ f(t)はいくつになる?
Review
Gauss函数のFourier変換
$ \int_Re^{-at^2}e^{-i\omega t}\mathrm dt=\sqrt{\frac\pi a}e^{-\frac{\omega^2}{4a}}
$ \int_Re^{-a\omega^2}e^{i\omega t}\mathrm dt=\sqrt{\frac\pi a}e^{-\frac{t^2}{4a}}
$ {\cal F}(f)(\omega)=e^{-a\omega^2}のとき、$ a\to0で裾が広がり、
$ {\cal F}(f)(\omega)\to 1\quad(a\to0)
になる
一方、$ f(t)=\sqrt{\frac\pi a}e^{-\frac{t^2}{4a}}は$ a\to 0でDiracのdelta函数$ \delta(t)になる
以上より
delta函数のFourier変換すると定値函数$ \omega\mapsto1になる
$ {\cal F}(\delta): \omega\mapsto1
ことがわかった
強制減衰振動モデルをFourier変換で解く#64a753fc1280f0000000e4e46cはDiracのdelta函数による急激な応答だったのであった
from ✅@2023-07-07T08:50D90 AM4-2023F-12
おさらいその2
$ \ddot x+2h\omega_0\dot x+{\omega_0}^2x=\frac{f(t)}{m}
をFourier変換して
$ (-\omega^2+2h\omega_0i\omega+{\omega_0^2}){\cal F}(x)(\omega)=\frac 1m{\cal F}(f)(\omega)
$ \iff {\cal F}(x):\omega\mapsto{\cal F}(G)(\omega){\cal F}(f)(\omega)
になった
本当は$ {\cal F}(x)をFourier逆変換したいが、ひとまず$ f=\deltaとして$ Gだけ求めた
$ {\cal F}(f)(\omega)=1なので、$ x(t)=G(t)となり簡単になる
ここから今回の話
定数函数をFourier逆変換するとdelta函数になるから、$ {\cal F}^{-1}(\omega\mapsto1)=\deltaだとわかった
さて、本当は$ {\cal F}^{-1}(\omega\mapsto{\cal F}(G)(\omega){\cal F}(f)(\omega))を求めたい
さっきは$ f=\deltaで求めたが、実際には任意の函数を使いたい
函数の積のFourier逆変換ということになるが、そんなことできるのだろうか?takker.icon
準備
ダブり防止のために、記号をちょい変える
$ {\cal F}(f)(\omega)=\int_\R f(\tau)e^{-i\omega\tau}\mathrm d\tau
導出
$ x(t)=\frac1{2\pi}\int_\R{\cal F}(G)(\omega){\cal F}(f)(\omega)e^{i\omega t}\mathrm d\omega
$ = \frac1{2\pi}\int_\R{\cal F}(G)(\omega)\int_\R f(\tau)e^{-i\omega\tau}\mathrm d\tau e^{i\omega t}\mathrm d\omega
$ = \int_\R f(\tau)\int_\R\frac1{2\pi}{\cal F}(G)(\omega)e^{i\omega(t-\tau)}\mathrm d\omega\mathrm d\tau
積分の順序交換の議論は端折った
$ = \int_\R f(\tau)G(t-\tau)\mathrm d\tau
$ = \int_\R G(t-\tau)f(\tau)\mathrm d\tau
この形式を畳み込み積分(convolution integral)と呼ぶらしい
$ \tauだけ遅れが発生する
もしかして強制振動の固有振動とのずれが関係していたりする?
合成積でも表せる
https://eman-physics.net/math/fourier08.html
よって
$ x: t\mapsto\int_\R G(t-\tau)f(\tau)\mathrm d\tau\\{\cal F}^{-1}\uparrow\quad\downarrow{\cal F}\\{\cal F}(x):\omega\mapsto{\cal F}(G)(\omega){\cal F}(f)(\omega)ー①
どんなinput$ f(t)があったとしても、Fourierの世界に持ち込むことで伝達函数$ \hat G(\omega)をかけるだけでoutputが求まるようになる
入力要素を排除し、システム固有の性質をつめこんだ、まさにシステムそのものを表すのがGreen函数$ Gの特徴であり意義である
①に$ f=\deltaを代入すると何が起こるかを考える
最初にやったやつと、畳み込み積分を使って導いたものとを比較すると、
$ x(t)= \int_\R G(t-\tau)f(\tau)\mathrm d\tau=\int_\R G(t-\tau)\delta(\tau)\mathrm d\tau=G(t)
ここから、Diracのdelta函数の重要な性質がわかる
$ \boxed{\int_\R f(u)\delta(t-u)\mathrm du=f(t)}
また①を使うと、inputがDiracのdelta函数(瞬間の衝撃力)のとき、outputがGreen函数になることもわかる
講義内容はここまで
次回以降はまとめと問題演習
設計つまらなかったけど、その道に進んだ
意外と役に立った。そういうこともある
耐震設計での重要性はもう少し後で話す
#2025-06-30 13:20:18
#2023-07-12 20:42:48
#2023-07-07 08:52:16
#2023-06-30 09:28:45
#2023-06-23 10:18:21