特異項をもつ BVP を解く

この例では、Emden 方程式を解く方法について説明します。Emden 方程式は、ガスの球体をモデリングするときに使われる特異項をもつ境界値問題です。

対称性を使用してモデルの PDE を減らすと、方程式は区間 $\left[0,1\right]$ で定義された 2 次 ODE になります。

$$y^{\prime \prime }+\frac{2}{x}{y}^{\prime }+y^5=0$$.

$\mathit{x}=0$ では、項 $\left(2/\mathit{x}\right)$ は特異項ですが、対称性は境界条件 ${\mathit{y}}^{\prime }\left(0\right)=0$ を意味しています。この境界条件では、項 $\left(2/\mathit{x}\right){\mathit{y}}^{\prime }$ は $\mathit{x}\to 0$ において適切に定義されています。境界条件 $\mathit{y}\left(1\right)=\sqrt{3}/2$ の場合、BVP には MATLAB® で計算される数値解と比較できる解析解があります。

$$y\left(x\right)={\left[\sqrt{1+\frac{x^2}{3}}\right]}^{-1}$$.

BVP ソルバー bvp4c は、次の形式をもつ特異な BVP を解くことができます。

$$y^{\prime }=S\frac{y}{x}+f(x,y)$$.

行列 $\mathit{S}$ は定数でなければならず、$\mathit{x}=0$ での境界条件は必要条件 $\mathit{S}\cdot \mathit{y}\left(0\right)=0$ と整合的でなければなりません。bvpset の 'SingularTerm' オプションを使用して、S 行列をソルバーに渡します。

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

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

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

境界条件は以下になります。

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

${\mathit{y}}_1\left(1\right)=\sqrt{3}/2$.

必須の行列形式では、方程式系は次のようになります。

$\left[\begin{array}{c}{{\mathit{y}}_1}^{\prime }\\ {{\mathit{y}}_2}^{\prime }\end{array}\right]=\frac{1}{\mathit{x}}\left[\begin{array}{cc}0& 0\\ 0& -2\end{array}\right]\left[\begin{array}{c}{\mathit{y}}_1\\ {\mathit{y}}_2\end{array}\right]+\left[\begin{array}{c}{\mathit{y}}_2\\ -{\mathit{y}}_1^5\end{array}\right]$.

行列形式では、$\mathit{S}=\left[\begin{array}{cc}0& 0\\ 0& -2\end{array}\right]$ および $\mathit{f}\left(\mathit{x},\mathit{y}\right)=\left[\begin{array}{c}{\mathit{y}}_2\\ -{\mathit{y}}_1^5\end{array}\right]$ であることが明らかです。

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

方程式のコード化

$\mathit{f}\left(\mathit{x},\mathit{y}\right)$ の方程式をコード化する関数を作成します。この関数にはシグネチャ dydx = emdenode(x,y) がなければなりません。ここでは以下のようになります。

• x は独立変数です。

• y は従属変数です。

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

これらの入力はソルバーによって自動的に関数に渡されますが、変数名によって方程式のコード化方法が決まります。この場合、

function dydx = emdenode(x,y)
dydx = [y(2) 
       -y(1)^5];
end

S を含む項はオプションを使用してソルバーに渡されるため、その項は関数に含まれません。

境界条件のコード化

次に、境界点で境界条件の残差値を返す関数を記述します。この関数にはシグネチャ res = emdenbc(ya,yb) がなければなりません。ここでは以下のようになります。

• ya は、区間の始まりにある境界条件の値です。

• yb は、区間の終わりにある境界条件の値です。

この問題の場合、境界条件の 1 つは ${\mathit{y}}_1$ に用いられ、もう 1 つは ${\mathit{y}}_2$ に用いられます。残差値を計算するには、境界条件を $\mathit{g}\left(\mathit{x},\mathit{y}\right)=0$ の形式にしなければなりません。

この形式では、境界条件は次となります。

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

${\mathit{y}}_1\left(1\right)-\sqrt{3}/2=0$.

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

function res = emdenbc(ya,yb)
res = [ya(2)
       yb(1) - sqrt(3)/2];
end

初期推定の作成

最後に、解の初期推定を作成します。この問題では、境界条件を満たす定数の初期推定と、0 ～ 1 の間にある 5 つの点のシンプルなメッシュを使用します。BVP ソルバーは解法プロセス中にこれらの点を適応させるため、多くのメッシュ点を使用する必要はありません。

${\mathit{y}}_1=\sqrt{3}/2$,

${\mathit{y}}_2=0$.

guess = [sqrt(3)/2; 0];
xmesh = linspace(0,1,5);
solinit = bvpinit(xmesh, guess);

方程式の求解

S の行列を作成して、それを 'SingularTerm' オプションの値として bvpset に渡します。最後に、ODE 関数、境界条件関数、初期推定、オプション構造体を使用して bvp4c を呼び出します。

S = [0 0; 0 -2];
options = bvpset('SingularTerm',S);
sol = bvp4c(@emdenode, @emdenbc, solinit, options);

解のプロット

$\left[0,1\right]$ における解析解の値を計算します。

x = linspace(0,1);
truy = 1 ./ sqrt(1 + (x.^2)/3);

比較のために、解析解と bvp4c によって計算された解をプロットします。

plot(x,truy,sol.x,sol.y(1,:),'ro');
title('Emden Problem -- BVP with Singular Term.')
legend('Analytical','Computed');
xlabel('x');
ylabel('solution y');

ローカル関数

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

function dydx = emdenode(x,y) % equation being solved 
dydx = [y(2) 
       -y(1)^5];
end
%-------------------------------------------
function res = emdenbc(ya,yb) % boundary conditions
res = [ya(2)
       yb(1) - sqrt(3)/2];
end
%-------------------------------------------

参考

bvp4c | bvp5c | bvpinit | bvpset

関連するトピック

• 境界値問題の求解