代数方程式系の求解

このトピックでは、Symbolic Math Toolbox™ を使用して方程式系をシンボリックに求解する方法を説明します。このツールボックスでは数値とシンボリックの両方の方程式ソルバーが提供されています。数値ソルバーとシンボリック ソルバーの比較については、数値ソルバーまたはシンボリック ソルバーの選択を参照してください。

solve の出力の処理

次の方程式系があり

$$\begin{array}{l}{\mathit{x}}^2{\mathit{y}}^2=0\\ \mathit{x}-\frac{\mathit{y}}{2}=\alpha ,\end{array}$$

これを $$x$$ と $$y$$ について解くことを考えます。はじめに、必要なシンボリック オブジェクトを作成します。

syms x y a

いくつかの方法で、solve の出力先を指定できます。1 つの方法は、2 つの出力を使って呼び出すことです。この呼び出しは次を返します。

[solx,soly] = solve(x^2*y^2 == 0, x-y/2 == a)

solx = 
$$\left(\begin{array}{c}a\\ 0\end{array}\right)$$

soly = 
$$\left(\begin{array}{c}0\\ -2\hspace{0.17em}a\end{array}\right)$$

最初の方程式を $$x^2{y}^2=1$$ に変更します。新しい方程式系には、より多くの解があります。4 つの異なる解が生成されました。

[solx,soly] = solve(x^2*y^2 == 1, x-y/2 == a)

solx = 
$$\left(\begin{array}{c}\frac{a}{2}-\frac{\sqrt{a^2-2}}{2}\\ \frac{a}{2}-\frac{\sqrt{a^2+2}}{2}\\ \frac{a}{2}+\frac{\sqrt{a^2-2}}{2}\\ \frac{a}{2}+\frac{\sqrt{a^2+2}}{2}\end{array}\right)$$

soly = 
$$\left(\begin{array}{c}-a-\sqrt{a^2-2}\\ -a-\sqrt{a^2+2}\\ \sqrt{a^2-2}-a\\ \sqrt{a^2+2}-a\end{array}\right)$$

従属変数を指定していないので、solve は symvar を使って変数を特定します。

solve からの出力を代入するこの方法は、"小さな" 方程式系に対しては、かなり有効です。たとえば、10 元 10 次の方程式系がある場合、次のように入力するのは面倒で時間が掛かります。

[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10] = solve(...)

この問題を解決するために、solve は、そのフィールドがそれぞれの解である構造体を出力できます。たとえば、方程式系 u^2 - v^2 = a^2、u + v = 1、a^2 - 2*a = 3 を解きます。ソルバーは、その結果を構造体に格納して返します。

syms u v a
S = solve(u^2 - v^2 == a^2, u + v == 1, a^2 - 2*a == 3)

S = struct with fields:
    a: [2x1 sym]
    u: [2x1 sym]
    v: [2x1 sym]

a の解は、S の "a フィールド" に帰属します。

S.a

ans = 
$$\left(\begin{array}{c}-1\\ 3\end{array}\right)$$

u と v の解も同様に解釈されます。これで、このフィールドとインデックスによる構造体 S の操作が可能になり、特定範囲の解にアクセスできます。たとえば、2 番目の解を調べるためには、次のステートメントを使用して、各フィールドの 2 番目の成分を取り出します。

s2 = [S.a(2),S.u(2),S.v(2)]

s2 = $$\left(\begin{array}{ccc}3& 5& -4\end{array}\right)$$

次のステートメントは、各行が個々の系の解である、解の行列 M を作成します。

M = [S.a,S.u,S.v]

M = 
$$\left(\begin{array}{ccc}-1& 1& 0\\ 3& 5& -4\end{array}\right)$$

solx と solyを消去して後で使えるようにします。

clear solx soly

線形方程式系の求解

線形方程式系は、行列の除算を使用して解くこともできます。たとえば、次の方程式系を解きます。

clear u v x y
syms u v x y
eqns = [x + 2*y == u, 4*x + 5*y == v];
S = solve(eqns);
sol = [S.x;S.y]

sol = 
$$\left(\begin{array}{c}\frac{2\hspace{0.17em}v}{3}-\frac{5\hspace{0.17em}u}{3}\\ \frac{4\hspace{0.17em}u}{3}-\frac{v}{3}\end{array}\right)$$

[A,b] = equationsToMatrix(eqns,x,y);
z = A\b

z = 
$$\left(\begin{array}{c}\frac{2\hspace{0.17em}v}{3}-\frac{5\hspace{0.17em}u}{3}\\ \frac{4\hspace{0.17em}u}{3}-\frac{v}{3}\end{array}\right)$$

このように、sol と z は、結果が異なる変数に代入されていますが、同じ解です。

方程式系の完全解を返す

solve は自動的には方程式のすべての解を返しません。すべての解を解のパラメーターおよび条件と共に返すには、オプション ReturnConditions を true に設定します。

次の方程式系を考えます。

$$\begin{array}{l}\sin \left(\mathit{x}\right)+\cos \left(\mathit{y}\right)=\frac{4}{5}\\ \sin \left(\mathit{x}\right)\cos \left(\mathit{y}\right)=\frac{1}{10}\end{array}$$

fimplicit を使用して方程式系を可視化します。x 軸と y 軸の値を pi に換算して指定するため、a の axes を使用して座標軸のハンドルを取得します。-2*pi から 2*pi までの値のシンボリック配列 S を pi/2 の間隔で作成します。目盛りを S に設定するため、a のプロパティ XTick と YTick を使用します。x 軸と y 軸のラベルを設定するため、S を文字ベクトルに変換します。arrayfun を使用して S のすべての要素に char を適用し、T を返します。a のプロパティ XTickLabel と YTickLabel を T に設定します。

syms x y
eqn1 = sin(x)+cos(y) == 4/5;
eqn2 = sin(x)*cos(y) == 1/10;
a = axes;
fimplicit(eqn1,[-2*pi 2*pi],'b');
hold on
grid on
fimplicit(eqn2,[-2*pi 2*pi],'m');
L = sym(-2*pi:pi/2:2*pi);
a.XTick = double(L);
a.YTick = double(L);
M = arrayfun(@char, L, 'UniformOutput', false);
a.XTickLabel = M;
a.YTickLabel = M;
title('Plot of System of Equations')
legend('sin(x)+cos(y) == 4/5','sin(x)*cos(y) == 1/10',...
    'Location','best','AutoUpdate','off')

解は 2 つのプロットの交点にあります。これは、この系には反復する周期解があることを示しています。この方程式系を解いてすべての解を求めるには、solve を使用し、ReturnConditions オプションを true に設定します。

S = solve(eqn1,eqn2,'ReturnConditions',true)

S = struct with fields:
             x: [2x1 sym]
             y: [2x1 sym]
    parameters: [z    z1]
    conditions: [2x1 sym]

solve は、x に対する解をもつ S.x、y に対する解をもつ S.y、解のパラメーターをもつ S.parameters、解の条件をもつ S.conditions の各フィールドから成る構造体 S を返します。1 つの解は S.x、S.y および S.conditions の同じインデックス番号をもつ要素により構成されます。したがって、S.x(1)、S.y(1) および S.conditions(1) により方程式系の解の 1 つが構成されます。S.parameters のパラメーターはすべての解に表れる場合があります。

S にインデックス付けして、解、パラメーターおよび条件を返します。

S.x

ans = 
$$\left(\begin{array}{c}{z}_1\\ z_1\end{array}\right)$$

S.y

ans = 
$$\left(\begin{array}{c}z\\ z\end{array}\right)$$

S.parameters

ans = $$\left(\begin{array}{cc}z& z_1\end{array}\right)$$

S.conditions

ans = 
$$\begin{array}{l}\left(\begin{array}{c}\left(\frac{z+\mathrm{acos}\left({\sigma }_3\right)}{2\hspace{0.17em}\pi }\in \mathbb{Z}\vee \frac{z-\mathrm{acos}\left({\sigma }_3\right)}{2\hspace{0.17em}\pi }\in \mathbb{Z}\right)\wedge \left(-\frac{\pi -z_1+{\sigma }_1}{2\hspace{0.17em}\pi }\in \mathbb{Z}\vee \frac{z_1+{\sigma }_1}{2\hspace{0.17em}\pi }\in \mathbb{Z}\right)\\ \left(\frac{z_1-\pi +\mathrm{asin}\left({\sigma }_3\right)}{2\hspace{0.17em}\pi }\in \mathbb{Z}\vee \frac{z_1-\mathrm{asin}\left({\sigma }_3\right)}{2\hspace{0.17em}\pi }\in \mathbb{Z}\right)\wedge \left(\frac{z+{\sigma }_2}{2\hspace{0.17em}\pi }\in \mathbb{Z}\vee \frac{z-{\sigma }_2}{2\hspace{0.17em}\pi }\in \mathbb{Z}\right)\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\mathrm{asin}\left(\frac{\sqrt{6}}{10}-\frac{2}{5}\right)\\ \\ \mathrm{ }{\sigma }_2=\mathrm{acos}\left(\frac{2}{5}-\frac{\sqrt{6}}{10}\right)\\ \\ \mathrm{ }{\sigma }_3=\frac{\sqrt{6}}{10}+\frac{2}{5}\end{array}$$

条件下における方程式系の求解

条件下で方程式系を求解するには、solve への入力に条件を指定します。

前述の方程式系において -2*pi から 2*pi の区間で x および y を求めます。scatter を使用して解をプロットに重ね合わせます。

Srange = solve(eqn1, eqn2, -2*pi < x, x < 2*pi, -2*pi < y, y < 2*pi, 'ReturnConditions', true);
scatter(Srange.x,Srange.y,'k')

solve によって返される解、パラメーター、条件の操作

solve によって返される解、パラメーター、条件を使用して、ある区間内またはその他の条件下における解を求めることができます。この節の目標は前節同様、検索範囲内で方程式系を解くことですが、アプローチが異なります。ここでは条件を直接指定する代わりに、solve によって返されるパラメーターと条件を使用する方法を説明します。

方程式系の完全解 S を求めるには、区間 -2*pi から 2*pi における x と y の値を求めます。そのためには、パラメーター S.parameters に対する解 S.x と S.y を、その区間内および条件 S.conditions 下において求めます。

区間内の x および y について求解する前に、assume を使用して S.conditions の条件を仮定し、返される解が条件を満たすようにします。1 番目の解に対する条件を仮定します。

assume(S.conditions(1))

S.x と S.y のパラメーターを求めます。

paramx = intersect(symvar(S.x),S.parameters)

paramx = $$z_1$$

paramy = intersect(symvar(S.y),S.parameters)

paramy = $$z$$

パラメーター paramx に対する x の 1 番目の解を求めます。

solparamx(1,:) = solve(S.x(1) > -2*pi, S.x(1) < 2*pi, paramx)

solparamx = 
$$\left(\begin{array}{cccc}\pi +\mathrm{asin}\left(\frac{\sqrt{6}}{10}-\frac{2}{5}\right)& \mathrm{asin}\left(\frac{\sqrt{6}}{10}-\frac{2}{5}\right)-\pi & -\mathrm{asin}\left(\frac{\sqrt{6}}{10}-\frac{2}{5}\right)& -2\hspace{0.17em}\pi -\mathrm{asin}\left(\frac{\sqrt{6}}{10}-\frac{2}{5}\right)\end{array}\right)$$

同様に、paramy に対する y の 1 番目の解を求めます。

solparamy(1,:) = solve(S.y(1) > -2*pi, S.y(1) < 2*pi, paramy)

solparamy = 
$$\left(\begin{array}{cccc}\mathrm{acos}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)& \mathrm{acos}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)-2\hspace{0.17em}\pi & -\mathrm{acos}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)& 2\hspace{0.17em}\pi -\mathrm{acos}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)\end{array}\right)$$

assume を使用して S.conditions(1) によって設定された仮定を消去します。asumptions を呼び出して仮定が消去されていることを確認します。

assume(S.parameters,'clear')
assumptions

ans = $$\left(\begin{array}{ccccccc}n\in \mathbb{Z}& 0<C& 0<L& 0<R_1& 0<R_2& 0<R_3& 0\le n\end{array}\right)$$

2 番目の解に対する条件を仮定します。

assume(S.conditions(2))

パラメーター paramx および paramy に対する x と y の 2 番目の解を求めます。

solparamx(2,:) = solve(S.x(2) > -2*pi, S.x(2) < 2*pi, paramx)

solparamx = 
$$\begin{array}{l}\left(\begin{array}{cccc}\pi +{\sigma }_1& {\sigma }_1-\pi & -{\sigma }_1& -2\hspace{0.17em}\pi -{\sigma }_1\\ {\sigma }_2& \pi -{\sigma }_2& {\sigma }_2-2\hspace{0.17em}\pi & -\pi -{\sigma }_2\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\mathrm{asin}\left(\frac{\sqrt{6}}{10}-\frac{2}{5}\right)\\ \\ \mathrm{ }{\sigma }_2=\mathrm{asin}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)\end{array}$$

solparamy(2,:) = solve(S.y(2) > -2*pi, S.y(2) < 2*pi, paramy)

solparamy = 
$$\begin{array}{l}\left(\begin{array}{cccc}{\sigma }_2& {\sigma }_2-2\hspace{0.17em}\pi & -{\sigma }_2& 2\hspace{0.17em}\pi -{\sigma }_2\\ {\sigma }_1& {\sigma }_1-2\hspace{0.17em}\pi & -{\sigma }_1& 2\hspace{0.17em}\pi -{\sigma }_1\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\mathrm{acos}\left(\frac{2}{5}-\frac{\sqrt{6}}{10}\right)\\ \\ \mathrm{ }{\sigma }_2=\mathrm{acos}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)\end{array}$$

paramx と paramy の 1 番目の行が方程式系の 1 番目の解を構成し、2 番目の行が 2 番目の解を構成します。

paramx および paramy のこれらの値に対する x と y の値を求めるため、subs を使用して S.x および S.y の paramx と paramy に代入します。

solx(1,:) = subs(S.x(1), paramx, solparamx(1,:));
solx(2,:) = subs(S.x(2), paramx, solparamx(2,:))

solx = 
$$\begin{array}{l}\left(\begin{array}{cccc}\pi +{\sigma }_1& {\sigma }_1-\pi & -{\sigma }_1& -2\hspace{0.17em}\pi -{\sigma }_1\\ {\sigma }_2& \pi -{\sigma }_2& {\sigma }_2-2\hspace{0.17em}\pi & -\pi -{\sigma }_2\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\mathrm{asin}\left(\frac{\sqrt{6}}{10}-\frac{2}{5}\right)\\ \\ \mathrm{ }{\sigma }_2=\mathrm{asin}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)\end{array}$$

soly(1,:) = subs(S.y(1), paramy, solparamy(1,:));
soly(2,:) = subs(S.y(2), paramy, solparamy(2,:))

soly = 
$$\begin{array}{l}\left(\begin{array}{cccc}{\sigma }_2& {\sigma }_2-2\hspace{0.17em}\pi & -{\sigma }_2& 2\hspace{0.17em}\pi -{\sigma }_2\\ {\sigma }_1& {\sigma }_1-2\hspace{0.17em}\pi & -{\sigma }_1& 2\hspace{0.17em}\pi -{\sigma }_1\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\mathrm{acos}\left(\frac{2}{5}-\frac{\sqrt{6}}{10}\right)\\ \\ \mathrm{ }{\sigma }_2=\mathrm{acos}\left(\frac{\sqrt{6}}{10}+\frac{2}{5}\right)\end{array}$$

solx および soly が x および y に対する 2 組の解である点に注目してください。方程式系のすべての解は、solx と soly において可能な値すべての組み合わせからなる、2 組の点集合です。

scatter を使用してこれら 2 組の点集合をプロットします。方程式のプロットの上に重ね合わせます。予想どおり、解は 2 つの方程式のプロットの交点に表れます。

for i = 1:length(solx(1,:))
    for j = 1:length(soly(1,:))
        scatter(solx(1,i), soly(1,j), 'k')
        scatter(solx(2,i), soly(2,j), 'k')
    end
end

シンボリックな結果の数値への変換

シンボリック計算は厳密な解を出しますが、数値計算は近似値を出します。この精度低下にもかかわらず、数値計算で使用するためにシンボリックな解を近似する数値に変換しなければならない場合があります。高精度な変換を行うには、関数vpaで提供される可変精度演算を使用します。標準的な精度でパフォーマンスを向上させるには、doubleを使用して倍精度値に変換します。

vpa を使用してシンボリック解 solx および soly を数値形に変換します。

vpa(solx)

ans = 
$$\left(\begin{array}{cccc}2.9859135500977407388300518406219& -3.2972717570818457380952349259371& 0.15567910349205249963259154265761& -6.1275062036875339772926952239014\\ 0.70095651347102524787213653614929& 2.4406361401187679905905068471302& -5.5822287937085612290531502304097& -3.8425491670608184863347799194288\end{array}\right)$$

vpa(soly)

ans = 
$$\left(\begin{array}{cccc}0.86983981332387137135918515549046& -5.4133454938557151055661016110685& -0.86983981332387137135918515549046& 5.4133454938557151055661016110685\\ 1.4151172233028441195987301489821& -4.8680680838767423573265566175769& -1.4151172233028441195987301489821& 4.8680680838767423573265566175769\end{array}\right)$$

複雑な結果の単純化によるパフォーマンスの改善

結果が複雑で solve が止まった場合、またはパフォーマンスを改善したい場合は、関数 solve で求めた方程式の解のトラブルシューティングを参照してください。