ラプラス変換を使用した微分方程式の求解

最新のリリースでは、このページがまだ翻訳されていません。このページの最新版は英語でご覧になれます。

このワークフローでは、Symbolic Math Toolbox™ でラプラス変換を使用して微分方程式を解きます。ラプラス変換の単純な例については、laplace および ilaplace を参照してください。

定義: ラプラス変換

関数 f(t) のラプラス変換は以下になります。

$F (s) = \int_{0}^{\infty} f (t) e^{- t s} d t .$

概念: シンボリックワークフローの使用

シンボリックワークフローは、計算を数値形式の代わりにもともとのシンボリック形式で保持します。この手法は、解のプロパティについて理解し、厳密なシンボリック値を使用するのに役立ちます。数値結果が必要な場合や、シンボリック的に続行できない場合に限り、シンボリック変数に数値を代入します。詳細は、数値演算またはシンボリック演算の選択を参照してください。通常、手順は以下になります。

方程式を宣言する。
方程式を解く。
値を代入する。
結果をプロットする。
結果を解析する。

ワークフロー: ラプラス変換を使用した RLC 回路の求解

方程式の宣言

ラプラス変換は、初期条件が指定された微分方程式の求解に利用できます。たとえば、以下のような抵抗-インダクタ-コンデンサ (RLC) 回路を解くことができます。

抵抗 (オーム): R₁、R₂、R₃
電流 (アンペア): I₁、I₂、I₃
インダクタンス (ヘンリー): L
静電容量 (ファラド): C
起電力 (ボルト): E(t)
電荷 (クーロン): Q(t)

キルヒホッフの電圧則および電流則を適用すると、RLC 回路の微分方程式系が得られます。

$\frac{d I_{1}}{d t} + \frac{R_{2}}{L} \frac{d Q}{d t} = \frac{R_{2} - R_{1}}{L} I_{1} .$

$\frac{d Q}{d t} = \frac{1}{R_{3} + R_{2}} (E (t) - \frac{1}{C} Q (t)) + \frac{R_{2}}{R_{3} + R_{2}} I_{1} .$

変数を宣言します。物理量は正の値であるため、対応する仮定を変数に設定します。このときの E(t) を 1 V の交流電圧とします。

syms L C I1(t) Q(t) s
R = sym('R%d',[1 3]);
assume([t L C R] > 0)
E(t) = 1*sin(t);        % Voltage = 1 V

微分方程式を宣言します。

dI1 = diff(I1,t);
dQ = diff(Q,t);
eqn1 = dI1 + (R(2)/L)*dQ == (R(2)-R(1))/L*I1
eqn2 = dQ == (1/(R(2)+R(3))*(E-Q/C)) + R(2)/(R(2)+R(3))*I1

eqn1(t) =
diff(I1(t), t) + (R2*diff(Q(t), t))/L == -(I1(t)*(R1 - R2))/L
eqn2(t) =
diff(Q(t), t) == (sin(t) - Q(t)/C)/(R2 + R3) + (R2*I1(t))/(R2 + R3)

初期値の電流と電圧 I₀ と Q₀ はどちらも 0 であると仮定します。これらの初期条件を宣言します。

cond1 = I1(0) == 0
cond2 = Q(0) == 0

cond1 =
I1(0) == 0
cond2 =
Q(0) == 0

方程式の求解

eqn1 と eqn2 のラプラス変換を計算します。

eqn1LT = laplace(eqn1,t,s)
eqn2LT = laplace(eqn2,t,s)

eqn1LT =
s*laplace(I1(t), t, s) - I1(0) - (R2*(Q(0) - s*laplace(Q(t), t, s)))/L == ...
-((R1 - R2)*laplace(I1(t), t, s))/L
eqn2LT =
s*laplace(Q(t), t, s) - Q(0) == (R2*laplace(I1(t), t, s))/(R2 + R3) + ...
(C/(s^2 + 1) - laplace(Q(t), t, s))/(C*(R2 + R3))

関数 solve は、シンボリック変数のみを解きます。したがって、solve を使用するには、最初に laplace(I1(t),t,s) と laplace(Q(t),t,s) を変数 I1_LT と Q_LT に置き換えます。

syms I1_LT Q_LT
eqn1LT = subs(eqn1LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])

eqn1LT =
I1_LT*s - I1(0) - (R2*(Q(0) - Q_LT*s))/L == -(I1_LT*(R1 - R2))/L

eqn2LT = subs(eqn2LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])

eqn2LT =
Q_LT*s - Q(0) == (I1_LT*R2)/(R2 + R3) - (Q_LT - C/(s^2 + 1))/(C*(R2 + R3))

I1_LT と Q_LT について方程式を解きます。

eqns = [eqn1LT eqn2LT];
vars = [I1_LT Q_LT];
[I1_LT, Q_LT] = solve(eqns,vars)

I1_LT =
(R2*Q(0) + L*I1(0) - C*R2*s + L*s^2*I1(0) + R2*s^2*Q(0) + C*L*R2*s^3*I1(0) + ...
C*L*R3*s^3*I1(0) + C*L*R2*s*I1(0) + C*L*R3*s*I1(0))/((s^2 + 1)*(R1 - R2 + L*s + ...
C*L*R2*s^2 + C*L*R3*s^2 + C*R1*R2*s + C*R1*R3*s - C*R2*R3*s))
Q_LT =
(C*(R1 - R2 + L*s + L*R2*I1(0) + R1*R2*Q(0) + R1*R3*Q(0) - R2*R3*Q(0) + ...
L*R2*s^2*I1(0) + L*R2*s^3*Q(0) + L*R3*s^3*Q(0) + R1*R2*s^2*Q(0) + R1*R3*s^2*Q(0) - ...
R2*R3*s^2*Q(0) + L*R2*s*Q(0) + ...
L*R3*s*Q(0)))/((s^2 + 1)*(R1 - R2 + L*s + C*L*R2*s^2 + C*L*R3*s^2 + ...
                C*R1*R2*s + C*R1*R3*s - C*R2*R3*s))

I1_LT と Q_LT の逆ラプラス変換を計算して、I₁ と Q を計算します。結果を単純化します。出力は長いため非表示にします。

I1sol = ilaplace(I1_LT,s,t);
Qsol = ilaplace(Q_LT,s,t);
I1sol = simplify(I1sol);
Qsol = simplify(Qsol);

値の代入

結果をプロットする前に、回路の各要素の数値でシンボリック変数を置き換えます。このとき、R1 = 4 Ω、R2 = 2 Ω、R3 = 3 Ω、C = 1/4 F、L = 1.6 H、I₁(0) = 15 A、Q(0) = 2 C とします。

vars = [R L C I1(0) Q(0)];
values = [4 2 3 1.6 1/4 15 2];
I1sol = subs(I1sol,vars,values)
Qsol = subs(Qsol,vars,values)

I1sol =
15*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) - ...
        (293*1001^(1/2)*sinh((1001^(1/2)*t)/40))/21879) - (5*sin(t))/51
Qsol =
(4*sin(t))/51 - (5*cos(t))/51 + (107*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) + ...
        (2039*1001^(1/2)*sinh((1001^(1/2)*t)/40))/15301))/51

結果のプロット

電流 I1sol と電荷 Qsol をプロットします。2 つの異なる時間間隔 0 ≤ t ≤ 10 および 5 ≤ t ≤ 25 を使用して、過渡と定常状態の両方の動作を示します。

subplot(2,2,1)
fplot(I1sol,[0 10])                      
title('Current')
ylabel('I1(t)')
xlabel('t')

subplot(2,2,2)
fplot(Qsol,[0 10])                        
title('Charge')
ylabel('Q(t)')
xlabel('t')

subplot(2,2,3)
fplot(I1sol,[5 25])                  
title('Current')
ylabel('I1(t)')
xlabel('t')
text(7,0.25,'Transient')
text(16,0.125,'Steady State')

subplot(2,2,4)
fplot(Qsol,[5 25])                        
title('Charge')
ylabel('Q(t)')
xlabel('t')
text(7,0.25,'Transient')
text(15,0.16,'Steady State')

結果の解析

最初は、電流と電荷は指数関数的に減衰します。しかし、長期的にはこれらは振動性です。これらの動作はそれぞれ "過渡" および "定常状態" と呼ばれます。シンボリックな結果からは結果の特性が解析できます。数値結果ではこれは不可能です。

I1sol と Qsol を目視で検査します。これらは、項の和です。children を使用して項を求めます。次に、[0 15] にわたって項をプロットして項の寄与を求めます。プロットは、過渡および定常状態の項を示します。

I1terms = children(I1sol);
Qterms = children(Qsol);

subplot(1,2,1)
fplot(I1terms,[0 15])
ylim([-2 2])
title('Current terms')

subplot(1,2,2)
fplot(Qterms,[0 15])
ylim([-2 2])
title('Charge terms')

このプロットは、I1sol が 1 つの過渡と 1 つの定常状態の項を持ち、Qsol が 1 つの過渡と 2 つの定常状態の項を持つことを示します。目視の検査により、I1sol と Qsol が exp 関数を含む項を持つことがわかります。この項によって過渡の指数減衰が発生すると推定します。has を使用して exp の項を確認し、I1sol と Qsol 内の過渡と定常状態の項を分離します。

I1transient = I1terms(has(I1terms,'exp'))
I1steadystate = I1terms(~has(I1terms,'exp'))

I1transient =
15*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) - (293*1001^(1/2)*sinh((1001^(1/2)*t)/40))/21879)
I1steadystate =
-(5*sin(t))/51

同様に、Qsol を過渡および定常状態の項に分離します。この結果は、シンボリック計算が問題の解析にどのように役立つかを示しています。

Qtransient = Qterms(has(Qterms,'exp'))
Qsteadystate = Qterms(~has(Qterms,'exp'))

Qtransient =
(107*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) + (2039*1001^(1/2)*sinh((1001^(1/2)*t)/40))/15301))/51
Qsteadystate =
[ -(5*cos(t))/51, (4*sin(t))/51]