連続性を使って BVP を解く

この例では、数値的に複雑な境界値問題を連続性によって解く方法を説明し、問題を比較的単純な一連の問題に効果的に分解します。

$0<\mathit{e}\ll 1$ の場合について、次の微分方程式を考えます。

$\mathit{e}{{\mathit{y}}^{\prime }}^{\prime }+{\mathrm{xy}}^{\prime }=-\mathit{e}{\pi }^2\cos \left(\pi \mathit{x}\right)-\pi \mathrm{xsin}\left(\pi \mathit{x}\right)$.

問題は区間 $\left[-1,1\right]$ に対して提示され、次の境界条件に従います。

$\mathit{y}\left(-1\right)=-2$,

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

$\mathit{e}={10}^{-4}$ の場合、方程式に対する解は $\mathit{x}=0$ の近くで急激に遷移するため、数値的に解くのは困難です。代わりに、この例では連続性を使用して、$\mathit{e}={10}^{-4}$ になるまで $\mathit{e}$ のいくつかの値を反復します。中間解はそれぞれ次の問題の初期推定として使用されます。

MATLAB® でこの方程式系を解くには、方程式、境界条件、および初期推定をコード化してから、境界値問題ソルバー bvp4c を呼び出す必要があります。(ここで行ったように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パスのディレクトリに保存することもできます。

方程式のコード化

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

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

${{\mathit{y}}_2}^{\prime }=-\frac{\mathit{x}}{\mathit{e}}{\mathit{y}}^{\prime }-{\pi }^2\cos \left(\pi \mathit{x}\right)-\frac{\pi \mathit{x}}{\mathit{e}}\sin \left(\pi \mathit{x}\right)$.

シグネチャ dydx = shockode(x,y) を使用して方程式をコード化する関数を記述します。ここでは以下のようになります。

• x は独立変数です。

• y は従属変数です。

• dydx(1) によって ${{\mathit{y}}_1}^{\prime }$ の方程式が与えられ、dydx(2) によって ${{\mathit{y}}_2}^{\prime }$ の方程式が与えられます。

shockode([x1 x2 ...],[y1 y2 ...]) が [shockode(x1,y1) shockode(x2,y2) ...] を返すように、関数をベクトル化します。この方法により、ソルバーのパフォーマンスが向上します。

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

function dydx = shockode(x,y)
pix = pi*x;
dydx = [y(2,:)
       -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)];
end

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

境界条件のコード化

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

$\mathit{y}\left(-1\right)+2=0$,

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

シグネチャ res = shockbc(ya,yb) を使用して境界条件をコード化する関数を記述します。ここでは以下のようになります。

• ya は、区間 $\left[\mathit{a},\mathit{b}\right]$ の始点における境界条件の値です。

• yb は、区間 $\left[\mathit{a},\mathit{b}\right]$ の終点における境界条件の値です。

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

function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end

ヤコビアンのコード化

この問題では、ODE 関数の解析的なヤコビアンと境界条件を簡単に計算できます。ヤコビアンを提供することで、ソルバーは有限差分を使用してヤコビアンを近似する必要がなくなるため、ソルバーの効率性が向上します。

ODE 関数の場合、ヤコビ行列は次のようになります。

${\mathit{J}}_{\mathrm{ODE}}=\frac{\partial \mathit{f}}{\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\\ 0& -\frac{\mathit{x}}{\mathit{e}}\end{array}\right]$.

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

function jac = shockjac(x,y,e)
jac = [0   1
       0  -x/e];
end

同様に、境界条件の場合、ヤコビ行列は次のようになります。

${\mathit{J}}_{\mathit{y}\left(\mathit{a}\right)}=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$ , ${\mathit{J}}_{\mathit{y}\left(\mathit{b}\right)}=\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right]$.

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

function [dBCdya,dBCdyb] = shockbcjac(ya,yb)
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end

初期推定の作成

$\left[-1,1\right]$ にある 5 つの点のメッシュで、解に対して定数の推定値を使用します。

sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);

方程式の求解

$\mathit{e}={10}^{-4}$ を使用して方程式の解を直接求めようとする場合、ソルバーは、遷移点 $\mathit{x}=0$ の近くで問題を適切に調整できないという状況を克服することができません。代わりに、$\mathit{e}={10}^{-4}$ に対する解を得るために、この例では ${10}^{-2}$、${10}^{-3}$、および ${10}^{-4}$ に対する一連の問題を解くことで連続性を使用します。各反復におけるソルバーからの出力は、次の反復で解の推定として機能します (このため、bvpinit からの初期推定の変数は sol であり、ソルバーからの出力にも sol という名前が付けられます)。

ヤコビアンの値は $\mathit{e}$ の値によって異なるため、ヤコビアンに対して関数 shockjac と関数 shockbcjac を指定するオプションをループ内に設定します。また、shockode が値のベクトルを処理するようにコード化されているため、ベクトル化をオンにします。

e = 0.1;
for i = 2:4
   e = e/10;
   options = bvpset('FJacobian',@(x,y) shockjac(x,y,e),...
                    'BCJacobian',@shockbcjac,...
                    'Vectorized','on');
   sol = bvp4c(@(x,y) shockode(x,y,e),@shockbc, sol, options);
end

解のプロット

メッシュ $\mathit{x}$ および解 $\mathit{y}\left(\mathit{x}\right)$ に対して、bvp4c の結果をプロットします。連続性を使用することで、ソルバーは $\mathit{x}=0$ での不連続性を処理することができます。

plot(sol.x,sol.y(1,:),'-o');
axis([-1 1 -2.2 2.2]);
title(['There Is a Shock at x = 0 When e =' sprintf('%.e',e) '.']);
xlabel('x');
ylabel('solution y');

ローカル関数

ここでは、BVP ソルバー bvp4c が解を計算するために呼び出すローカル関数を紹介しています。あるいは、これらの関数を独自のファイルとして MATLAB パスのディレクトリに保存することもできます。

function dydx = shockode(x,y,e) % equation to solve
pix = pi*x;
dydx = [y(2,:)
       -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)];
end
%-------------------------------------------
function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end
%-------------------------------------------
function jac = shockjac(x,y,e) % jacobian of shockode
jac = [0   1
       0  -x/e];
end
%-------------------------------------------
function [dBCdya,dBCdyb] = shockbcjac(ya,yb) % jacobian of shockbc
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end
%-------------------------------------------

参考

bvp4c | bvpinit | bvpset

関連するトピック

• 境界値問題の求解
• 連続性を使用した BVP の整合性の検証