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

微分方程式系の求解

複数の変数をもつ複数の常微分方程式系を、dsolve を用いて初期条件の指定に関係なく解きます。一元微分方程式の解を求めるには、微分方程式の求解を参照してください。

微分方程式系の求解

この線形 1 階微分方程式系を解きます。

dudt=3u+4v,dvdt=4u+3v.

はじめに、u と v を syms を用いて表現し、シンボリック関数 u(t)v(t) を作成します。

syms u(t) v(t)

== を使用して方程式を定義し、関数 diff を使用して微分を表します。

ode1 = diff(u) == 3*u + 4*v;
ode2 = diff(v) == -4*u + 3*v;
odes = [ode1; ode2]
odes(t) =
 diff(u(t), t) == 3*u(t) + 4*v(t)
 diff(v(t), t) == 3*v(t) - 4*u(t)

構造体の要素として解を返す関数 dsolve を使用して系を解きます。

S = dsolve(odes)
S = 
  struct with fields:

    v: [1×1 sym]
    u: [1×1 sym]

dsolve を使用して方程式を解くことができない場合、方程式を数値的に解いてみてください。2 階微分方程式の数値的な求解を参照してください。

u(t)v(t) にアクセスするため、構造体 S でインデックスを指定します。

uSol(t) = S.u
vSol(t) = S.v
uSol(t) =
C2*cos(4*t)*exp(3*t) + C1*sin(4*t)*exp(3*t)
vSol(t) =
C1*cos(4*t)*exp(3*t) - C2*sin(4*t)*exp(3*t)

または、複数の出力引数を指定して、u(t)v(t) を直接格納します。

[uSol(t), vSol(t)] = dsolve(odes)
uSol(t) =
C2*cos(4*t)*exp(3*t) + C1*sin(4*t)*exp(3*t)
vSol(t) =
C1*cos(4*t)*exp(3*t) - C2*sin(4*t)*exp(3*t)

条件が指定されていないため、定数 C1C2 が出現します。初期条件 u(0) == 0v(0) == 0 で系を解きます。これらの条件を満たす定数の値を、関数 dsolve が求めます。

cond1 = u(0) == 0;
cond2 = v(0) == 1;
conds = [cond1; cond2];
[uSol(t), vSol(t)] = dsolve(odes,conds)
uSol(t) =
sin(4*t)*exp(3*t)
vSol(t) =
cos(4*t)*exp(3*t)

fplot を使用して解を可視化します。

fplot(uSol)
hold on
fplot(vSol)
grid on
legend('uSol','vSol','Location','best')

行列形式での微分方程式の求解

dsolve を使用して、微分方程式を行列形式で解きます。

この微分方程式系を考えます。

dxdt=x+2y+1,dydt=x+y+t.

この系の行列形式は次になります。

[x'y']=[1211][xy]+[1t].

以下により

Y=[xy],A=[1211],B=[1t].

この系は、Y′ = AY + B となります。

これらの行列と行列方程式を定義します。

syms x(t) y(t)
A = [1 2; -1 1];
B = [1; t];
Y = [x; y];
odes = diff(Y) == A*Y + B
odes(t) =
  diff(x(t), t) == x(t) + 2*y(t) + 1
   diff(y(t), t) == t - x(t) + y(t)

dsolve を使用して行列方程式を解きます。関数 simplify を使用して解を単純化します。

[xSol(t), ySol(t)] = dsolve(odes);
xSol(t) = simplify(xSol(t))
ySol(t) = simplify(ySol(t))
xSol(t) =
(2*t)/3 + 2^(1/2)*C2*exp(t)*cos(2^(1/2)*t) + 2^(1/2)*C1*exp(t)*sin(2^(1/2)*t) + 1/9
ySol(t) =
C1*exp(t)*cos(2^(1/2)*t) - t/3 - C2*exp(t)*sin(2^(1/2)*t) - 2/9

条件が指定されていないため、定数 C1C2 が出現します。

初期条件 u(0) = 2v(0) = –1 で系を解きます。方程式を行列形式で指定する場合、初期条件も行列形式で指定しなければなりません。dsolve で、これらの条件を満たす定数の値を求められます。

C = Y(0) == [2; -1];
[xSol(t), ySol(t)] = dsolve(odes,C)
xSol(t) =
(2*t)/3 + (17*exp(t)*cos(2^(1/2)*t))/9 - (7*2^(1/2)*exp(t)*sin(2^(1/2)*t))/9 + 1/9
ySol(t) =
- t/3 - (7*exp(t)*cos(2^(1/2)*t))/9 - (17*2^(1/2)*exp(t)*sin(2^(1/2)*t))/18 - 2/9

fplot を使用して解を可視化します。

clf
fplot(ySol)
hold on
fplot(xSol)
grid on
legend('ySol','xSol','Location','best')

参考