MATLABでもっとも単純な数値積分と数値微分を実装する

伝達関数

\begin{align}
P=\frac{1}{s^2+1}
\end{align}

の単純な数値積分と数値微分を考える。

始めに確認のため、連続時間での結果を確認する。連続時間での微分と積分はラプラス演算子を用いて

\begin{align}
P_I=P C_2=\frac{1}{s(s^2+1)} \hspace{5mm} P_D=P C_1=\frac{s}{s^2+1}
\end{align}

結果

コード

figure;
subplot(2,1,1)
impulse(P*C1,10,'k')
grid on
subplot(2,1,2)
impulse(P*C2,10,'k')
grid on

次にこの問題を数値積分と数値微分をで考える。離散時間での１ステップの信号の遅延を\(z^{-1}\)と定義し、これを使って遅延演算子\(\delta=1-z^{-1}\)定義する。遅延演算子を使えば連続時間の微分と積分はそれぞれ

\begin{align}
\frac{1}{s}= \frac{Ts}{\Delta} \hspace{5mm} s= \frac{\Delta}{Ts}
\end{align}

と表現できる。今回は連続時間での伝達関数を離散化することで同様の結果を得ている。

結果

ここで理想的なインパルス応答を得るにはデルタ関数の定義

\begin{align}
\delta(x)=\begin{cases}
\infty \hspace{10mm} &(x=0) \\
0 &(x \neq 0)
\end{cases}
\end{align}

に基づき与えるのが好ましいが、今回は十分大きな値を与えることでデルタ関数とした。

全体のコード

Ts=0.001;
s=tf('s');

P=1/(s^2+1);
Pd=c2d(1/(s^2+1),Ts,'zoh')

C1=s;
C2=1/s;

t=linspace(0,10,10/Ts);
r=zeros(1,10/Ts);
r(1,10)=1000;

figure;
subplot(2,1,1)
impulse(P*C1,10,'k')
grid on
subplot(2,1,2)
impulse(P*C2,10,'k')
grid on

y1=lsim(P*C1,r,t);
y2=lsim(P*C2,r,t);

figure;
subplot(2,1,1)
plot(t,y1,'k')
grid on
subplot(2,1,2)
plot(t,y2,'k')
grid on