二次システムの自由振動の一般解

#概要
二次システムの自由振動，過減衰の応答を計算します．

#二次システム
```math
G(s) = A \frac{w_n^2}{s^2 +2 \zeta wn s+w_n^2} 
```
が二次システムです．$A$は振幅で一般的なゲインです．
ステップ入力した場合，最終的に$A$に収束します．

```math
\frac{Y(s)}{U(s)} = A \frac{w_n^2}{s^2 +2 \zeta w_n s+w_n^2} \\
\therefore (s^2 +2 \zeta w_n s+w_n^2)Y(s) = A w_n^2 U(s)\\
\therefore y''(t) +2 \zeta w_n y(t)'+w_n^2y(t) = A w_n^2 u(t)\\
```
を必ず満たす．

#自由振動
u=0,で初期値が0でないときの応答を求めます．

```math
y''(t) +2 \zeta w_n y(t)'+w_n^2y(t) = 0\\
```
を満たす．

```math
y(t) =  Be^{-\zeta w_n t}sin(\sqrt{1-\zeta^2}w_nt+\phi)
```
が解になる，$ B,\phi$が任意定数

###検証

```math
y'(t) =  Be^{-\zeta w_n t}sin(\sqrt{1-\zeta^2}w_nt+\phi) (-\zeta w_n)\\
+  Be^{-\zeta w_n}cos(\sqrt{1-\zeta^2}w_nt+\phi) \sqrt{1-\zeta^2}w_n

```

```math
y''(t) =  Be^{-\zeta w_n t}sin(\sqrt{1-\zeta^2}w_nt+\phi) (-\zeta w_n)^2\\
-  Be^{-\zeta w_n}sin(\sqrt{1-\zeta^2}w_nt+\phi) (1-\zeta^2)w_n^2\\

+ Be^{-\zeta w_n t}cos(\sqrt{1-\zeta^2}w_nt+\phi) (-\zeta w_n)\sqrt{1-\zeta^2}w_n\\
+ Be^{-\zeta w_n}cos(\sqrt{1-\zeta^2}w_nt+\phi) (-\zeta w_n)\sqrt{1-\zeta^2}w_n\\

= Be^{-\zeta w_n t}sin(\sqrt{1-\zeta^2}w_nt+\phi) 
\{\zeta^2 w_n^2 - (1-\zeta^2)w_n^2

\}
+\\
 Be^{-\zeta w_n t}cos(\sqrt{1-\zeta^2}w_nt+\phi) 
\{
2(-\zeta w_n)\sqrt{1-\zeta^2}w_n
\}

```
微分方程式に代入し

```math
y''(t) +2 \zeta w_n y(t)'+w_n^2y(t) \\
=
 Be^{-\zeta w_n t}sin(\sqrt{1-\zeta^2}w_nt+\phi) 
\{ \zeta^2 w_n^2 - (1-\zeta^2)w_n^2
+2 \zeta w_n (-\zeta w_n) +w_n^2
\}\\
+
 Be^{-\zeta w_n t}cos(\sqrt{1-\zeta^2}w_nt+\phi) 
\{
2(-\zeta w_n)\sqrt{1-\zeta^2}w_n
+2 \zeta w_n \sqrt{1-\zeta^2}w_n
\}\\
=0
```
となり，成立する

#初期値

```math
y(0) = Bsin(\phi)\\
y'(0)  =  B w_n sin(\phi + \pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )
```
この二つの式を$ B,\phi$に対し解くことで，解が得られる．

```math
B=\frac{y(0)}{sin \phi}
```
代入し

```math
y'(0)  =  \frac{y(0)}{sin \phi}　 w_n sin(\phi + \pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )\\

\therefore
\frac{y'(0)}{y(0)w_n}=cos(\pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )+\frac{sin(\pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )}{tan \phi}
```
よって，

```math
tan \phi = \frac{sin(\pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )}{\frac{y'(0)}{y(0)w_n}-cos(\pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )}
```
で計算できる

```math
B=\frac{y(0)}{sin \phi}
```
で完了

ただし，$y(0) = 0$だと0割が発生するので，場合分けがいる
まぁ，$\phi = 0$に確定するので簡単

```math
B= \frac{y'(0)}{   w_n sin( \pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )}
```
で完了


#Scilabで検証
csimと理論解の比較を行った

![hikaku.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/fc4253dd-2952-5932-8b3a-8335f98b144c.png)
ちゃんと同じになった．

```c

clear();
for i = 1:10 do
    close
end

zeta = 0.5
wn = 3
a = 1

A=[-2*zeta*wn -1*wn^2
1 0]

B=[a*wn^2
0]

C=[0 1]
D=[0]

x0=[1
3]

sl = syslin('c',A,B,C,D,x0);


t=0:0.01:10;
u=zeros(1,length(t))
y=csim(u,t,sl)

beta1 = %pi - atan(sqrt(1-zeta^2)/zeta);

if x0(2)==0 then
    phi = 0
    b=x0(1)/wn/sin(beta1)
else
    phi = atan(sin(beta1)/(x0(1)/x0(2)/wn - cos(beta1) ));
    b=x0(2)/sin(phi);
end


y2= diag(b*exp(-zeta*wn*t')*sin(sqrt(1-zeta^2)*wn*t +phi) )
clf();
plot(t,y);
figure();
plot(t,y2);
```

#自由振動まとめ

```math

y(t) =  Be^{-\zeta w_n t}sin(\sqrt{1-\zeta^2}w_nt+\phi)
```
初期状態に対し$B,\phi$を次の式で調整

```math
tan \phi = \frac{sin(\pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )}{\frac{y'(0)}{y(0)w_n}-cos(\pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )}\\
B=\frac{y(0)}{sin \phi}
```
$y(0) = 0$の時のみ例外的に

```math
\phi = 0\\
B= \frac{y'(0)}{   w_n sin( \pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta} )}
```

さらに，$y'(t)$の方を観測できる場合は

```math
y'(t) =  Bw_ne^{-\zeta w_n t}sin(\sqrt{1-\zeta^2}w_nt+\phi + \pi -tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta}) \\
```

となる．
こちらも，scilabで検証したが問題なしだった