このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
bvp5c
境界値問題の求解 — 5 次法
説明
例
2 次 BVP の求解
MATLAB® で関数を使用して 2 次 BVP を解きます。この例では、次の 2 次方程式を使用します。
.
この方程式は、区間 で次の境界条件に従って定義されます。
,
.
この方程式を MATLAB で解くには、この方程式を 1 次方程式系として表す関数、境界条件の関数および初期推定の関数を作成する必要があります。次に、BVP ソルバーがこれら 3 つの入力を使用して方程式を解きます。
方程式のコード化
方程式をコード化する関数を作成します。代入 および を使用して、方程式を 1 次方程式系として書き換えます。
,
.
対応する関数は次のようになります。
function dydx = bvpfcn(x,y) dydx = zeros(2,1); dydx = [y(2) -y(1)]; end
メモ: 関数はすべて、この例の最後にローカル関数として含まれています。
境界条件のコード化
境界条件を の形式でコード化する関数を記述します。この形式では、境界条件は次となります。
,
.
対応する関数は次のようになります。
function res = bcfcn(ya,yb) res = [ya(1) yb(1)-2]; end
初期推定の作成
関数 bvpinit
を使用して、方程式の解の初期推定を作成します。この方程式は を に関連付けるため、解には三角関数が含まれるという推定が妥当です。積分区間で 5 点からなるメッシュを使用します。メッシュの最初と最後の値に、ソルバーが境界条件を適用します。
初期推定の関数は、 を入力として受け入れ、 および の推定値を返します。関数は次のようになります。
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 を考えます。
.
この方程式は、区間 で次の境界条件に従って定義されます。
,
.
この方程式を MATLAB® で解くには、この方程式を 1 次方程式系として表す関数と境界条件の関数を作成し、いくつかのオプション値を設定して、初期推定を作成する必要があります。次に、BVP ソルバーはこれら 4 つの入力を使用して方程式を解きます。
方程式のコード化
代入 および を使用して、ODE を 1 次方程式系として書き換えることができます。
,
.
対応する関数は次のようになります。
function dydx = bvpfcn(x,y) dydx = [y(2) -2*y(2)/x - y(1)/x^4]; end
メモ: 関数はすべて、この例の最後にローカル関数として含まれています。
境界条件のコード化
境界条件関数では、境界条件が の形式でなければなりません。この形式では、境界条件は次のようになります。
,
.
対応する関数は次のようになります。
function res = bcfcn(ya,yb) res = [ya(1) yb(1)-sin(1)]; end
オプションの設定
bvpset
を使用してソルバー統計の表示を有効にし、粗い許容誤差を指定して、ソルバー間の誤差制御の差異を強調表示します。また、効率化のために、次の解析ヤコビアンを指定します。
.
これに対応し、ヤコビアンの値を返す関数は次のようになります。
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
を使用して、解の初期推定を作成します。区間 に 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.
結果のプロット
の 2 通りの計算結果を、比較のため解析解と共にプロットします。解析解は次のとおりです。
,
.
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)
を使用します。戻り値 dydt
は、データ型が single
または double
で f(x,y) に対応する列ベクトルです。odefun
は、x
と y
のいずれかの引数が関数で使用されない場合でも、両方の入力引数を受け入れなければなりません。
たとえば、 を解くには、次の関数を使用します。
function dydt = odefun(t,y) dydt = 5*y-3; end
方程式系の場合、odefun
の出力はベクトルです。ベクトルの各要素は 1 つの方程式の解です。たとえば、
を解くには、次の関数を使用します。
function dydt = odefun(t,y) dydt = zeros(2,1); dydt(1) = y(1)+2*y(2); dydt(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
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
— 解の初期推定
構造体
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
には次のフィールドがあります。
フィールド | 説明 |
---|---|
|
|
|
|
|
|
|
|
|
|
| 解に関連する計算コスト統計 (メッシュ点の数、残差誤差、 |
詳細
多点境界値問題
多点境界値問題では、積分区間のいくつかの点に境界条件が適用されます。
bvp5c
は、a = a0 < a1 < a2 < ...< an = b が区間 [a,b] 内の境界点である多点境界値問題を解くことができます。点 a1,a2,...,an−1 は、[a,b] を領域に分割する接点です。bvp5c
は、1 から始まるインデックスで左から右へ (a から b へ) と領域を列挙します。領域 k、すなわち [ak−1,ak], において、bvp5c
は以下のように微分係数を評価します。
yp = odefun(x,y,k)
境界条件関数 bcfun(yleft,yright)
において、yleft(:,k)
は、[ak−1,ak] の左境界における解です。同様に、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
を伴う問題も含め、特異境界値問題のクラスを解きます。
区間は、[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 次の精度をもつ C1 連続解を与えます。この公式は、陰的ルンゲ・クッタ公式として実装されます。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 で導入
MATLAB コマンド
次の MATLAB コマンドに対応するリンクがクリックされました。
コマンドを MATLAB コマンド ウィンドウに入力して実行してください。Web ブラウザーは MATLAB コマンドをサポートしていません。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)