方程式の数値的な求解

Symbolic Math Toolbox™ では数値とシンボリックの両方の方程式ソルバーが提供されています。数値ソルバーとシンボリック ソルバーの比較については、数値ソルバーまたはシンボリック ソルバーの選択を参照してください。方程式または連立方程式には複数の解がある場合があります。これらの解を数値的に求めるには、関数 vpasolve を使用します。多項方程式の場合、vpasolve はすべての解を返します。非多項方程式では、vpasolve は最初に見つかった解を返します。この例では vpasolve を使用して多項方程式および非多項方程式の解を求め、これらの解を任意の精度で取得する方法を説明します。

多項式関数のすべての根を求める

vpasolve を使用して関数 $$f(x)=6x^7-2x^6+3x^3-8$$ のすべての解を求めます。

syms f(x)
f(x) = 6*x^7-2*x^6+3*x^3-8;
sol = vpasolve(f)

sol = 
$$\left(\begin{array}{c}1.0240240759053702941448316563337\\ -0.88080620051762149639205672298326+0.50434058840127584376331806592405\hspace{0.17em}\mathrm{i}\\ -0.88080620051762149639205672298326-0.50434058840127584376331806592405\hspace{0.17em}\mathrm{i}\\ -0.22974795226118163963098570610724+0.96774615576744031073999010695171\hspace{0.17em}\mathrm{i}\\ -0.22974795226118163963098570610724-0.96774615576744031073999010695171\hspace{0.17em}\mathrm{i}\\ 0.7652087814927846556172932675903+0.83187331431049713218367239317121\hspace{0.17em}\mathrm{i}\\ 0.7652087814927846556172932675903-0.83187331431049713218367239317121\hspace{0.17em}\mathrm{i}\end{array}\right)$$

vpasolve は予想どおり関数の 7 つの根を返します。これは関数が 7 次の多項式であるためです。

検索範囲と開始点を使用して非多項式関数の零点を求める

関数 $$f(x)=e^{(x/7)}\cos (2x)$$ のプロットを見ると零点が周期的に存在し、$$x$$ が増加するにつれて零点における傾きが増加していることが分かります。

syms x
h = fplot(exp(x/7)*cos(2*x),[-2 25]);
grid on

vpasolve を使用して関数 f の零点を求めます。vpasolve は複数の解が存在する場合でも、非多項式方程式の 1 つの解のみを返す点に注意してください。繰り返しの vpasolve 呼び出しにおいて、同じ結果が返されます。

f = exp(x/7)*cos(2*x);
for k = 1:3
  vpasolve(f,x)
end

ans = $$-7.0685834705770347865409476123789$$

ans = $$-7.0685834705770347865409476123789$$

ans = $$-7.0685834705770347865409476123789$$

複数の解を求めるには、オプション 'Random' を true に設定します。これにより vpasolve は開始点をランダムに選択するようになります。ランダムな開始点の選択にかかわるアルゴリズムについての詳細は vpasolve ページのアルゴリズムを参照してください。

for k = 1:3
  vpasolve(f,x,'Random',true)
end

ans = $$-226.98006922186256147892598444194$$

ans = $$98.174770424681038701957605727484$$

ans = $$52.621676947629036744249276669932$$

$$x=10$$ に近い零点を求めるには、開始点を 10 にします。

vpasolve(f,x,10)

ans = $$10.210176124166828025003590995658$$

$$x=1000$$ に近い零点を求めるには、開始点を 1000 にします。

vpasolve(f,x,1000)

ans = $$999.8118620049516981407362567287$$

$$15\le x\le 25$$ の範囲にある零点を求めるには、検索範囲を [15 25] に設定します。

vpasolve(f,x,[15 25])

ans = $$21.205750411731104359622842837137$$

範囲 [15 25] にある複数の零点を求める場合、vpasolve を繰り返し呼び出すことはできません。前述のように毎回同じ結果を返すためです。代わりに、検索範囲を設定し、'Random' を true に設定します。

for k = 1:3
  vpasolve(f,x,[15 25],'Random',true)
end

ans = $$21.205750411731104359622842837137$$

ans = $$21.205750411731104359622842837137$$

ans = $$16.493361431346414501928877762217$$

'Random' は開始点をランダムに選択するため、連続的な呼び出しで同じ解が得られる場合があります。

指定した検索範囲内のすべての零点を求める

関数 findzeros を作成し、指定された検索範囲および許容誤差における f のすべての零点を体系的に求めます。この関数は入力された検索範囲から開始して vpasolve を呼び出すことにより零点を求めます。そして、検索範囲を零点値の周囲で 2 つに分割し、新しい検索範囲を入力として関数自身を再帰的に呼び出し、他の零点を求めます。

ここでは、この関数をセクションごとに説明します。

3 つの入力と 1 つの出力をもつ関数を宣言します。1 番目の入力は関数、2 番目の入力は範囲、オプションである 3 番目の入力には零点と零点から生成された上限および下限の間のエラーを指定できます。

function sol = findzeros(f,range,err)

オプションの引数でエラー許容値を指定しない場合、findzeros は err を 0.001 に設定します。

if nargin < 2
  err = 1e-3;
end

vpasolve を使用して検索範囲内の零点を求めます。

sol = vpasolve(f,range);

vpasolve が零点を見つけられない場合は終了します。

if(isempty(sol))
  return

vpasolve が零点を見つけた場合は、その零点の上と下の 2 つの検索範囲に検索範囲を分割します。

else
  lowLimit = sol-err;
  highLimit = sol+err;

下の検索範囲に対して findzeros を呼び出します。findzeros が零点を返した場合は、値を解の配列にコピーし並べ替えます。

  temp = findzeros(f,[range(1) lowLimit],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end

上の検索範囲に対して findzeros を呼び出します。findzeros が零点を返した場合は、値を解の配列にコピーし並べ替えます。

  temp = findzeros(f,[highLimit range(2)],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end
  return
end
end

関数 findzeros の全体は以下のとおりです。この関数に findzeros.m という名前を付けて現在のフォルダーに保存します。

function sol = findzeros(f,range,err)
if nargin < 3
  err = 1e-3;
end
sol = vpasolve(f,range);
if(isempty(sol))
  return
else
  lowLimit = sol-err;
  highLimit = sol+err;
  temp = findzeros(f,[range(1) lowLimit],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end
  temp = findzeros(f,[highLimit range(2)],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end
  return
end
end

検索範囲 [15 25] に対して findzeros を呼び出し、f(x) = exp(x/7)*cos(2*x) の同範囲に存在する零点すべてを、既定の許容誤差内で求めます。

syms f(x)
f(x) = exp(x/7)*cos(2*x);
sol = findzeros(f,[15 25])'

sol = 
$$\left(\begin{array}{c}16.493361431346414501928877762217\\ 18.064157758141311121160199453857\\ 19.634954084936207740391521145497\\ 21.205750411731104359622842837137\\ 22.776546738526000978854164528776\\ 24.347343065320897598085486220416\end{array}\right)$$

任意の精度による解の取得

digits を使用して、vpasolve によって返される解の精度を設定します。既定では、vpasolve は有効桁数 32 桁の精度で解を返します。

f = exp(x/7)*cos(2*x);
vpasolve(f)

ans = $$-7.0685834705770347865409476123789$$

digits を使用して精度を有効桁数 64 桁に高めます。digits の値を変更する際は、必ず現在の値を保存し復元できるようにしておきます。

digitsOld = digits;
digits(64)
vpasolve(f)

ans = $$-7.068583470577034786540947612378881489443631148593988097193625333$$

次に、解の精度を有効桁数 16 桁に変更します。

digits(16)

検索範囲を使用した多変数方程式の求解

次の連立方程式を考えます。

$$\begin{array}{l}\mathit{z}=10\left(\cos \left(\mathit{x}\right)+\cos \left(\mathit{y}\right)\right)\\ \mathit{z}=\mathit{x}+{\mathit{y}}^2-0.1{\mathit{x}}^2\mathit{y}\\ \mathit{x}+\mathit{y}-2.7=0\end{array}$$

方程式 $$0\le x\le 2.5$$ と $$0\le x\le 2.5$$ のプロットを見ると 3 つの表面が 2 つの点で交差していることがわかります。プロットを見やすくするため、viewを使用して視線を変更します。

syms x y z
eqn1 = z == 10*(cos(x) + cos(y));
eqn2 = z == x+y^2-0.1*x^2*y;
eqn3 = x+y-2.7 == 0;
equations = [eqn1 eqn2 eqn3];
fimplicit3(equations)
axis([0 2.5 0 2.5 -20 10])
title('System of Multivariate Equations')
view(69,28)

vpasolve を使用して表面が交差する点を求めます。関数 vpasolve が構造体を返します。解の値 x、y、z にアクセスするため、構造体でインデックスを指定します。

sol = vpasolve(equations);
[sol.x sol.y sol.z]

ans = $$\left(\begin{array}{ccc}2.369747722454798& 0.3302522775452021& 2.293354376823228\end{array}\right)$$

解空間の範囲検索を行うため、変数の検索範囲を指定します。範囲 $$0\le x\le 1.5$$ と $$1.5\le y\le 2.5$$ を指定すると、関数 vpasolve は以下に示す区域を検索します。

vpasolve を使用してこの検索範囲における解を求めます。$$z$$ の検索範囲を除外するには、3 番目の検索範囲を [NaN NaN] に設定します。

vars = [x y z];
range = [0 1.5; 1.5 2.5; NaN NaN];
sol = vpasolve(equations,vars,range);
[sol.x sol.y sol.z]

ans = $$\left(\begin{array}{ccc}0.9106266172563336& 1.789373382743666& 3.964101572135625\end{array}\right)$$

複数の解を求めるには、オプション 'Random' を true に設定します。これにより、vpasolve は連続的な実行においてランダムな開始点を使用します。'Random' オプションを検索範囲と併せて使用することにより、vpasolve で検索範囲内のランダムな開始点を使用できます。'Random' は開始点をランダムに選択するため、連続的な呼び出しで同じ解が得られる場合があります。vpasolve を繰り返し呼び出して必ず両方の解が見つかるようにします。

clear sol
range = [0 3; 0 3; NaN NaN];
for k = 1:5
  temp = vpasolve(equations,vars,range,'Random',true);
  sol(k,1) = temp.x;
  sol(k,2) = temp.y;
  sol(k,3) = temp.z;
end
sol

sol = 
$$\left(\begin{array}{ccc}2.369747722454798& 0.3302522775452021& 2.293354376823228\\ 2.369747722454798& 0.3302522775452021& 2.293354376823228\\ 2.369747722454798& 0.330252277545202& 2.293354376823228\\ 0.9106266172563336& 1.789373382743666& 3.964101572135625\\ 0.9106266172563336& 1.789373382743666& 3.964101572135625\end{array}\right)$$

方程式をプロットします。scatter3を使用して、解を黄色の X マーカーが付いた点から成る散布図として重ね合わせます。プロットを見やすくするため、alphaを使用して 2 つの表面を透明にし、viewを使用して視線を変更します。

図のように、vpasolve は方程式によって形成される表面の交点で解を見つけます。

clf
ax = axes;
h = fimplicit3(equations);
h(2).FaceAlpha = 0;
h(3).FaceAlpha = 0;
axis([0 2.5 0 2.5 -20 10])
hold on
scatter3(sol(:,1),sol(:,2),sol(:,3),600,'red','X','LineWidth',2)
title('Randomly Found Solutions in Specified Search Range')
cz = ax.Children;
view(69,28)
hold off

最後に、digits の古い値に戻して計算を続けます。

digits(digitsOld)