bvp4c

境界値問題の求解 — 4 次法

構文

sol = bvp4c(odefun,bcfun,solinit)

sol = bvp4c(odefun,bcfun,solinit,options)

説明

sol = bvp4c(odefun,bcfun,solinit) は、odefun で指定された y′ = f(x,y) の形式の連立微分方程式を、bcfun で表される境界条件と解の初期推定 solinit に従って積分します。関数 bvpinit を使用して初期推定 solinit を作成します。これにより、bcfun の境界条件が適用される点も定義されます。

sol = bvp4c(odefun,bcfun,solinit,options) は options (関数 bvpset を使用して作成された引数) で定義された積分設定も使用します。たとえば、AbsTol オプションおよび RelTol オプションを使用して絶対許容誤差と相対許容誤差を指定したり、FJacobian オプションを使用して odefun の解析的な偏導関数を指定したりできます。

例

2 次 BVP の求解

MATLAB® で関数を使用して 2 次 BVP を解きます。この例では、次の 2 次方程式を使用します。

${{\mathit{y}}^{\prime }}^{\prime }+\mathit{y}=0$.

この方程式は、区間 $\left[0,\pi /2\right]$ で次の境界条件に従って定義されます。

$\mathit{y}\left(0\right)=0$,

$\mathit{y}\left(\pi /2\right)=2$.

この方程式を MATLAB で解くには、この方程式を連立 1 階方程式として表す関数、境界条件の関数および初期推定の関数を作成する必要があります。次に、BVP ソルバーがこれら 3 つの入力を使用して方程式を解きます。

方程式のコード化

方程式をコード化する関数を作成します。代入 ${\mathit{y}}_1=\mathit{y}$ および ${\mathit{y}}_2={\mathit{y}}^{\prime }$ を使用して、方程式を連立 1 階方程式として書き換えます。

${{\mathit{y}}_1}^{\prime }={\mathit{y}}_2$,

${{\mathit{y}}_2}^{\prime }=-{\mathit{y}}_1$.

対応する関数は次のようになります。

function dydx = bvpfcn(x,y)
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end

メモ: 関数はすべて、この例の最後にローカル関数として含まれています。

境界条件のコード化

境界条件を $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$ の形式でコード化する関数を記述します。この形式では、境界条件は次となります。

$\mathit{y}\left(0\right)=0$,

$\mathit{y}\left(\pi /2\right)-2=0$.

対応する関数は次のようになります。

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-2];
end

初期推定の作成

関数 bvpinit を使用して、方程式の解の初期推定を作成します。この方程式は ${{\mathit{y}}^{\prime }}^{\prime }$ を $\mathit{y}$ に関連付けるため、解には三角関数が含まれるという推定が妥当です。積分区間で 5 点からなるメッシュを使用します。メッシュの最初と最後の値に、ソルバーが境界条件を適用します。

初期推定の関数は、$\mathit{x}$ を入力として受け入れ、${\mathit{y}}_1$ および ${\mathit{y}}_2$ の推定値を返します。関数は次のようになります。

function g = guess(x)
g = [sin(x)
     cos(x)];
end

xmesh = linspace(0,pi/2,5);
solinit = bvpinit(xmesh, @guess);

方程式の求解

導関数、境界条件関数および初期推定を指定して bvp4c を使用し、問題を解きます。

sol = bvp4c(@bvpfcn, @bcfcn, solinit);

解のプロット

plot(sol.x, sol.y, '-o')

ローカル関数

ここでは、bvp4c が方程式を解くために使用するローカル関数のリストを示します。

function dydx = bvpfcn(x,y) % equation to solve
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end
%--------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-2];
end
%--------------------------------
function g = guess(x) % initial guess for y and y'
g = [sin(x)
     cos(x)];
end
%--------------------------------

bvp4c ソルバーと bvp5c ソルバーの比較

2 つの異なるソルバーにより粗い許容誤差で BVP を解き、その結果を比較します。

次の 2 次 ODE を考えます。

${{\mathit{y}}^{\prime }}^{\prime }+\frac{2}{\mathit{x}}{\mathit{y}}^{\prime }+\frac{1}{{\mathit{x}}^4}\mathit{y}=0$.

この方程式は、区間 $\left[\frac{1}{3\pi },1\right]$ で次の境界条件に従って定義されます。

$\mathit{y}\left(\frac{1}{3\pi }\right)=0$,

$\mathit{y}\left(1\right)=\sin \left(1\right)$.

この方程式を MATLAB® で解くには、この方程式を連立 1 階方程式として表す関数と境界条件の関数を作成し、いくつかのオプション値を設定して、初期推定を作成する必要があります。次に、BVP ソルバーはこれら 4 つの入力を使用して方程式を解きます。

方程式のコード化

代入 ${\mathit{y}}_1=\mathit{y}$ および ${\mathit{y}}_2={\mathit{y}}^{\prime }$ を使用して、ODE を連立 1 階方程式として書き換えることができます。

${{\mathit{y}}_1}^{\prime }={\mathit{y}}_2$,

${{\mathit{y}}_2}^{\prime }=-\frac{2}{\mathit{x}}{\mathit{y}}_2-\frac{1}{{\mathit{x}}^4}{\mathit{y}}_1$.

対応する関数は次のようになります。

function dydx = bvpfcn(x,y)
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end

メモ: 関数はすべて、この例の最後にローカル関数として含まれています。

境界条件のコード化

境界条件関数では、境界条件が $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$ の形式でなければなりません。この形式では、境界条件は次のようになります。

$\mathit{y}\left(\frac{1}{3\pi }\right)=0$,

$\mathit{y}\left(1\right)-\sin \left(1\right)=0$.

対応する関数は次のようになります。

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-sin(1)];
end

オプションの設定

bvpset を使用してソルバー統計の表示を有効にし、粗い許容誤差を指定して、ソルバー間の誤差制御の差異を強調表示します。また、効率化のために、次の解析ヤコビアンを指定します。

$\mathit{J}=\frac{\partial {\mathit{f}}_{\mathit{i}}}{\partial \mathit{y}}=\left[\begin{array}{cc}\frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_2}\\ \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_2}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -\frac{1}{{\mathit{x}}^4}& -\frac{2}{\mathit{x}}\end{array}\right]$.

これに対応し、ヤコビアンの値を返す関数は次のようになります。

function dfdy = jac(x,y)
dfdy = [0      1
       -1/x^4 -2/x];
end

opts = bvpset('FJacobian',@jac,'RelTol',0.1,'AbsTol',0.1,'Stats','on');

初期推定の作成

bvpinit を使用して、解の初期推定を作成します。区間 $\left[1/3\pi ,1\right]$ に 10 点からなる初期メッシュをもつ初期推定として、定数関数を指定します。

xmesh = linspace(1/(3*pi), 1, 10);
solinit = bvpinit(xmesh, [1; 1]);

方程式の求解

bvp4c と bvp5c の両方を使用して方程式を解きます。

sol4c = bvp4c(@bvpfcn, @bcfcn, solinit, opts);

The solution was obtained on a mesh of 9 points.
The maximum residual is  9.794e-02. 
There were 157 calls to the ODE function. 
There were 28 calls to the BC function. 

sol5c = bvp5c(@bvpfcn, @bcfcn, solinit, opts);

The solution was obtained on a mesh of 11 points.
The maximum error is  6.742e-02. 
There were 244 calls to the ODE function. 
There were 29 calls to the BC function. 

結果のプロット

${\mathit{y}}_1$ の 2 通りの計算結果を、比較のため解析解と共にプロットします。解析解は次のとおりです。

${\mathit{y}}_1=\sin \left(\frac{1}{\mathit{x}}\right)$,

${\mathit{y}}_2=-\frac{1}{{\mathit{x}}^2}\cos \left(\frac{1}{\mathit{x}}\right)$.

xplot = linspace(1/(3*pi),1,200);
yplot = [sin(1./xplot); -cos(1./xplot)./xplot.^2];

plot(xplot,yplot(1,:),'k',sol4c.x,sol4c.y(1,:),'r*',sol5c.x,sol5c.y(1,:),'bo')
title('Comparison of BVP Solvers with Crude Error Tolerance')
legend('True','BVP4C','BVP5C')
xlabel('x')
ylabel('solution y')

このプロットにより、bvp5c は計算の真の誤差を直接制御し、bvp4c は間接的にのみ制御していることが確認されます。許容誤差がより厳密である場合、ソルバー間の差異はそれほど明確にはなりません。

ローカル関数

ここでは、BVP ソルバーが問題を解くために使用するローカル関数のリストを示します。

function dydx = bvpfcn(x,y) % equation to solve
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end
%---------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-sin(1)];
end
%---------------------------------
function dfdy = jac(x,y) % analytical jacobian for f
dfdy = [0      1
       -1/x^4 -2/x];
end
%---------------------------------

入力引数

odefun — 解を求める関数

関数ハンドル

解を求める関数。被積分関数を定義する関数ハンドルとして指定します。odefun と bcfun は同じ数の入力引数を受け入れなければなりません。

odefun をコード化するには、スカラー x と列ベクトル y の関数シグネチャ dydx = odefun(x,y) を使用します。戻り値 dydx は、データ型が single または double で f(x,y) に対応する列ベクトルです。odefun は、x と y のいずれかの引数が関数で使用されない場合でも、両方の入力引数を受け入れなければなりません。

たとえば、$$y\text{'}=5y-3$$ を解くには、次の関数を使用します。

function dydx = odefun(x,y)
dydx = 5*y-3;
end

連立方程式の場合、odefun の出力はベクトルです。ベクトルの各要素は 1 つの方程式の解です。たとえば、

$$\begin{array}{l}y{\text{'}}_1=y_1+2y_2\\ y{\text{'}}_2=3y_1+2y_2\end{array}$$

を解くには、次の関数を使用します。

function dydx = odefun(x,y)
dydx = zeros(2,1);
dydx(1) = y(1)+2*y(2);
dydx(2) = 3*y(1)+2*y(2);
end

bvp4c は、解の特異性または多点境界条件をもつ問題を解くこともできます。

例: sol = bvp4c(@odefun, @bcfun, solinit)

未知パラメーター

解く BVP に未知パラメーターが含まれる場合、関数シグネチャ dydx = odefun(x,y,p) を代わりに使用することができます。ここで、p はパラメーター値のベクトルです。この関数シグネチャを使用する場合、BVP ソルバーは、solinit で指定されたパラメーター値の初期推定から順に、未知パラメーターを計算します。

データ型: function_handle

bcfun — 境界条件

関数ハンドル

境界条件。境界条件の残差誤差を計算する関数ハンドルとして指定します。odefun と bcfun は同じ数の入力引数を受け入れなければなりません。

bcfun をコード化するには、列ベクトル ya および yb の関数シグネチャ res = bcfun(ya,yb) を使用します。戻り値 res は、データ型が single または double で、境界点における境界条件の残差値に対応する列ベクトルです。

たとえば、y(a) = 1 と y(b) = 0 の場合、境界条件関数は次のようになります。

function res = bcfun(ya,yb)
res = [ya(1)-1
       yb(1)];
end

y(a) = 1 であるため、点 x = a における ya(1)-1 の残差値は 0 になります。同様に、y(b) = 0 であるため、点 x = b における yb(1) の残差値は 0 になります。

境界条件の適用される境界点 x = a と x = b は、初期推定構造体 solinit で定義されます。2 点境界値問題の場合、a = solinit.x(1) および b = solinit.x(end) になります。

例: sol = bvp4c(@odefun, @bcfun, solinit)

未知パラメーター

解く BVP に未知パラメーターが含まれる場合、関数シグネチャ res = bcfun(ya,yb,p) を代わりに使用することができます。ここで、p はパラメーター値のベクトルです。この関数シグネチャを使用する場合、BVP ソルバーは、solinit で指定されたパラメーター値の初期推定から順に、未知パラメーターを計算します。

データ型: function_handle

solinit — 解の初期推定

構造体

解の初期推定。構造体として指定します。solinit の作成には bvpinit を使用します。

初期値問題とは異なり、境界値問題には、解がない、有限個の解がある、または解が無限にある、ということがあります。BVP を解くプロセスでは、必要な解の推定を与えることが重要な部分になります。この推定の質はソルバーのパフォーマンスや、ときには問題が解けるかどうかにも大きく影響します。良好な初期推定を作成するためのガイドラインについては、解の初期推定を参照してください。

例: sol = bvp4c(@odefun, @bcfun, solinit)

データ型: struct

options — オプション構造体

構造体

オプション構造体。オプション構造体の作成または変更を行うには、関数 bvpset を使用します。

例: options = bvpset('RelTol',1e-5,'Stats','on') は相対許容誤差 1e-5 を指定し、ソルバー統計の表示をオンにします。

データ型: struct

出力引数

sol — 解の構造体

構造体

解の構造体。sol のフィールドにはドット インデックス (sol.field1 など) を使用してアクセスできます。解 (sol.x,sol.y) は初期メッシュ solinit.x で定義された積分区間上において連続で、かつその 1 階微分も連続です。sol を関数 deval で使用して、区間内の他の点における解を評価できます。

構造体 sol には次のフィールドがあります。

フィールド	説明
x	bvp4c で選択されたメッシュ。このメッシュは通常、初期メッシュ solinit.x とは異なる点を含みます。
y	sol.x のメッシュ点における y(x) の近似。
yp	sol.x のメッシュ点における y′(x) の近似。
parameters	solinit.parameters で指定された未知パラメーターの最終値。
solver	'bvp4c'
stats	解に関連する計算コスト統計 (メッシュ点の数、残差誤差、odefun と bcfun の呼び出しの数)。

詳細

多点境界値問題

多点境界値問題では、積分区間のいくつかの点に境界条件が適用されます。

bvp4c は、a = a₀ < a₁ < a₂ < ...< a_n = b が区間 [a,b] 内の境界点である多点境界値問題を解くことができます。点 a₁,a₂,...,a_n−1 は、[a,b] を領域に分割する接点です。bvp4c は、1 から始まるインデックスで左から右へ (a から b へ) と領域を列挙します。領域 k、すなわち [a_k−1,a_k], において、bvp4c は以下のように微分係数を評価します。

yp = odefun(x,y,k) 

境界条件関数 bcfun(yleft,yright) において、yleft(:,k) は、[a_k−1,a_k] の左境界における解です。同様に、yright(:,k) は領域 k の右境界における解です。特に、yleft(:,1) = y(a) であり、また yright(:,end) = y(b) です。

bvpinit で初期推定を作成するときには、個々の境界点について xinit で double 要素を使用します。詳細については、bvpinit のリファレンス ページを参照してください。

yinit が関数の場合、関数 bvpinit は y = yinit(x,k) を呼び出し、領域 k の x における解に対する初期推定を得ます。関数 bpv4c が返す解の構造体 sol では、sol.x が個々の境界面を表す double 要素です。sol.y の対応する列には、境界面における左側および右側の解がそれぞれ格納されます。

3 点境界値問題の求解例については、複数の境界条件をもつ BVP を解くを参照してください。

特異境界値問題

関数 bvp4c は、以下の形式の未知のパラメーター p を伴う問題も含め、特異境界値問題のクラスを解きます。

$$\begin{array}{l}y\text{'}=S\frac{y}{x}+f\left(x,y,p\right),\\ 0=bc\left(y\left(0\right),y\left(b\right),p\right).\end{array}$$

区間は、[0, b] で、b > 0 でなければなりません。このような問題は、円柱あるいは球面対称による偏微分方程式 (PDE) の結果である ODE の滑らかな解を計算する場合にしばしば発生します。特異問題の場合は、関数 bvpset の 'SingularTerm' オプションの値として (定数) 行列 S を指定します。odefun は、f(x,y,p) のみを計算します。境界条件と初期推定は、平滑化に関する必要条件 S·y(0) = 0 と一致しなければなりません。

特異境界値問題の求解例については、特異項をもつ BVP を解くを参照してください。

アルゴリズム

関数 bvp4c は、3 段 Lobatto IIIa 公式を実装する有限差分コードです[1]、[2]。これは選点法の方程式で選点多項式は、シミュレーション区間で一様な 4 次精度の C¹ 連続解を与えます。メッシュの選択とエラー制御は、連続解の残差をベースにしています。

この選点法ではメッシュ点を使って、積分区間を部分区間に分割します。ソルバーは、すべての部分区間に適用される選点条件と境界条件から求まるグローバルな連立代数方程式を数値的に解いて解を求めます。そしてソルバーは、各部分区間における数値解の誤差を計算します。解が許容誤差の基準を満たさない場合、ソルバーはメッシュの適合を行ってから、再度処理を繰り返します。ユーザーは、初期メッシュの点と、それらのメッシュ点における解の初期近似を "指定しなければなりません"。

参照

[1] Shampine, L.F., and J. Kierzenka. "A BVP Solver based on residual control and the MATLAB PSE." ACM Trans. Math. Softw. Vol. 27, Number 3, 2001, pp. 299–316.

[2] Shampine, L.F., M.W. Reichelt, and J. Kierzenka. "Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c." MATLAB File Exchange, 2004.

拡張機能

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

この関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

バージョン履歴

R2006a より前に導入

参考

bvp5c | bvpget | bvpinit | bvpset | bvpxtend | deval

トピック

• 境界値問題の求解