bvp5c

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

構文

sol = bvp5c(odefun,bcfun,solinit)

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

説明

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

sol = bvp5c(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);

方程式の求解

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

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

解のプロット

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

ローカル関数

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

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

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

例: sol = bvp5c(@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 = bvp5c(@odefun, @bcfun, solinit)

未知パラメーター

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

データ型: function_handle

solinit — 解の初期推定

構造体

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

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

例: sol = bvp5c(@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	bvp5c で選択されたメッシュ。このメッシュは通常、初期メッシュ solinit.x とは異なる点を含みます。
y	sol.x のメッシュ点における y(x) の近似。
yp	sol.x のメッシュ点における y′(x) の近似。
parameters	solinit.parameters で指定された未知パラメーターの最終値。
solver	'bvp5c'
stats	解に関連する計算コスト統計 (メッシュ点の数、残差誤差、odefun と bcfun の呼び出しの数)。

詳細

多点境界値問題

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

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

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 を解くを参照してください。

特異境界値問題

関数 bvp5c は、以下の形式の未知のパラメーター 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 を解くを参照してください。

アルゴリズム

関数 bvp5c は、4 段 Lobatto IIIa 公式を実装する有限差分コードです [1]。これは、選点公式で、選点多項式が [a,b] 間で一様に 5 次の精度をもつ C¹ 連続解を与えます。この公式は、陰的ルンゲ・クッタ公式として実装されます。bvp5c と bvp4c の差異をいくつか示します。

• bvp5c は代数方程式を直接解き、bvp4c は解析的な凝縮を使用します。

• bvp4c は未知パラメーターを直接処理します。bvp5c は未知パラメーターに関する自明の微分方程式により連立方程式を増大させます。

参照

[1] Shampine, L.F., and J. Kierzenka. "A BVP Solver that Controls Residual and Error." J. Numer. Anal. Ind. Appl. Math. Vol. 3(1-2), 2008, pp. 27–41.

拡張機能

スレッドベースの環境

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

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

バージョン履歴

R2006b で導入

参考

bvp4c | bvpget | bvpinit | bvpset | bvpxtend | deval

トピック

• 境界値問題の求解