最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。
このワークフローでは、Symbolic Math
Toolbox™ でラプラス変換を使用して微分方程式を解きます。ラプラス変換の単純な例については、laplace
および ilaplace
を参照してください。
関数 f(t) のラプラス変換は以下になります。
シンボリック ワークフローは、計算を数値形式の代わりにもともとのシンボリック形式で保持します。この手法は、解のプロパティについて理解し、厳密なシンボリック値を使用するのに役立ちます。数値結果が必要な場合や、シンボリック的に続行できない場合に限り、シンボリック変数に数値を代入します。詳細は、数値演算またはシンボリック演算の選択を参照してください。通常、手順は以下になります。
方程式を宣言する。
方程式を解く。
値を代入する。
結果をプロットする。
結果を解析する。
ラプラス変換は、初期条件が指定された微分方程式の求解に利用できます。たとえば、以下のような抵抗-インダクタ-コンデンサ (RLC) 回路を解くことができます。
抵抗 (オーム): R1、R2、R3
電流 (アンペア): I1、I2、I3
インダクタンス (ヘンリー): L
静電容量 (ファラド): C
起電力 (ボルト): E(t)
電荷 (クーロン): Q(t)
キルヒホッフの電圧則および電流則を適用すると、RLC 回路の微分方程式系が得られます。
変数を宣言します。物理量は正の値であるため、対応する仮定を変数に設定します。このときの 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)
初期値の電流と電圧 I0 と Q0 はどちらも 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
の逆ラプラス変換を計算して、I1 と 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、I1(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]