メインコンテンツまでスキップ

数値積分手法の比較シミュレーション

前進差分、後退差分、双一次変換をそれぞれ用いて、例題1.2.1の電気回路のステップ応答を数値計算で求め、数値計算の結果と真の応答をプロットし、比較する。

電気回路モデル​

RL回路の微分方程式：

Ldi(t)dt+Ri(t)=u(t)L\frac{di(t)}{dt} + Ri(t) = u(t)

真の解（ステップ応答）：

i(t)=1R(1−e−RLt)i(t) = \frac{1}{R}(1 - e^{-\frac{R}{L}t})

数値積分手法​

1. 前進差分（Forward Difference）​

微分の近似：

sf(n)=f(n+1)−f(n)Tsf(n) = \frac{f(n+1)-f(n)}{T}

よって：

Li(n+1)−i(n)T+Ri(n)=u(n)L\frac{i(n+1)-i(n)}{T} + Ri(n) = u(n)

整理すると：

i(n+1)=1L(Tu(n)+(L−RT)i(n))i(n+1) = \frac{1}{L}(Tu(n) + (L-RT)i(n))

時間をシフトして：

i(n)=1L(Tu(n−1)+(L−RT)i(n−1))i(n) = \frac{1}{L}(Tu(n-1) + (L-RT)i(n-1))

2. 後退差分（Backward Difference）​

微分の近似：

sf(n)=f(n)−f(n−1)Tsf(n) = \frac{f(n)-f(n-1)}{T}

よって：

Li(n)−i(n−1)T+Ri(n)=u(n)L\frac{i(n)-i(n-1)}{T} + Ri(n) = u(n)

整理すると：

i(n)=1RT+L(Tu(n)+Li(n−1))i(n) = \frac{1}{RT+L}(Tu(n) + Li(n-1))

3. 双一次変換（Bilinear Transform）​

台形積分による近似：

sf(n)=2T1−z−11+z−1f(n)sf(n) = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}f(n)

微分方程式に適用：

L2T(i(n)−i(n−1))+R(i(n)+i(n−1))=u(n)+u(n−1)L\frac{2}{T}(i(n)-i(n-1)) + R(i(n)+i(n-1)) = u(n)+u(n-1)

整理すると：

i(n)=12L+RT(Tu(n)+Tu(n−1)+(2L−RT)i(n−1))i(n) = \frac{1}{2L+RT} \left( Tu(n)+Tu(n-1) + (2L-RT)i(n-1) \right)

MATLABシミュレーションコード​

初期設定​

clear
close all
set(0,'DefaultAxesFontName','Times New Roman')
set(0,'DefaultAxesFontSize',12)

%--------------------------------------------------------------------
%  同定ゼミ宿題 電気回路応答のシミュレーション
%                                                4/21 by YANG
%---------------------------------------------------------------------

パラメータ設定​

Tmax = 10;           % 最大時間
samp = 0.5;          % サンプリング周期
T = samp;
n = Tmax/samp;       % サンプル数
u = ones(1,n+1);     % ステップ入力

% 回路パラメータ
R = 1;               % 抵抗値
L = 1;               % インダクタンス値

配列初期化​

% 数列
t  = zeros(1, n+1);  % 時間配列
i  = zeros(1, n+1);  % 真の信号
i1 = zeros(1, n+1);  % 前進差分
i2 = zeros(1, n+1);  % 後退差分
i3 = zeros(1, n+1);  % 双一次変換

各手法による数値計算​

真の解析解​

% 真の信号
for k = 2:n+1
    t(k) = (k-1) * samp;
    i(k) = (1 - exp((-t(k)*R)/L)) / R;
end

前進差分​

% 前進差分
for k = 2:n+1
    t(k) = (k-1) * samp;
    i1(k) = (samp*u(k-1) + (L-R*samp)*i1(k-1)) / L;
end

後退差分​

% 後退差分
for k = 2:n+1
    t(k) = (k-1) * samp;
    i2(k) = (T*u(k) + L*i2(k-1)) / (R*samp+L);
end

双一次変換​

% 双一次変換
for k = 2:n+1
    t(k) = (k-1) * samp;
    i3(k) = (T*u(k) + T*u(k-1) + (2*L-R*T)*i3(k-1)) / (2*L+R*samp);
end

結果のプロット​

% 比較プロット
subplot(3,1,1); 
plot(t,i,t,i1); 
title('前進差分');
legend('真の解','前進差分');  
xlabel('t'); 
ylabel('i(t)');
grid on

subplot(3,1,2); 
plot(t,i,t,i2); 
title('後退差分');
legend('真の解','後退差分');      
xlabel('t'); 
ylabel('i(t)');
grid on

subplot(3,1,3); 
plot(t,i,t,i3); 
title('双一次変換');
legend('真の解','双一次変換');  
xlabel('t'); 
ylabel('i(t)');
grid on

数値積分手法の特性比較​

1. 前進差分の特性​

特徴：

  • 明示的（explicit）手法
  • 計算が簡単
  • 安定性に制限がある

安定性条件： ∣1−RTL∣≤1|1 - \frac{RT}{L}| \leq 1

これより： T≤2LRT \leq \frac{2L}{R}

安定性の制限

前進差分は大きなサンプリング周期で不安定になる可能性がある。

2. 後退差分の特性​

特徴：

  • 暗示的（implicit）手法
  • 無条件安定
  • 計算がやや複雑

安定性： 分母がRT+L>0RT + L > 0で常に正のため、無条件安定。

安定性の優位性

後退差分は大きなサンプリング周期でも安定である。

3. 双一次変換の特性​

特徴：

  • 台形積分による近似
  • 優れた周波数特性
  • バランスの取れた精度

安定性： 連続時間システムが安定なら、離散化後も安定性が保持される。

周波数特性の保持

双一次変換は周波数応答特性を最もよく保持する。

精度と安定性の評価​

サンプリング周期の影響​

手法精度安定性計算コスト
前進差分低〜中条件付き低
後退差分中無条件中
双一次変換高無条件中

適用場面​

前進差分​

  • 高速な計算が必要
  • 小さなサンプリング周期
  • リアルタイム処理

後退差分​

  • 安定性が重要
  • 大きなサンプリング周期
  • ロバスト性が要求される場合

双一次変換​

  • 高精度が必要
  • 周波数特性の保持が重要
  • 制御系設計

まとめ​

本シミュレーションにより以下のことが確認できる：

  1. 精度比較: 双一次変換が最も高精度
  2. 安定性: 後退差分と双一次変換が優秀
  3. 計算効率: 前進差分が最も簡単
  4. 実用性: 双一次変換が総合的に優れている
手法選択の指針
  • 高精度が必要 → 双一次変換
  • 安定性重視 → 後退差分
  • 高速計算 → 前進差分（安定性条件下）

これらの数値積分手法の理解は、システム同定や制御系設計において重要な基礎知識となる。