未知パラメーターをもつ BVP を解く

この例では bvp4c を使用して未知パラメーターをもつ境界値問題を解く方法を説明します。

Mathieu 方程式は、区間 $\left[0,\pi \right]$ で次により定義されます。

${{\mathit{y}}^{\prime }}^{\prime }+\left(\lambda -2\mathit{q}\cos \left(2\mathit{x}\right)\right)\mathit{y}=0$.

パラメーターが $\mathit{q}=5$ の場合、境界条件は次のようになります。

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

${\mathit{y}}^{\prime }\left(\pi \right)=0$.

ただし、これは $\mathit{y}\left(\mathit{x}\right)$ を定数の倍数までしか決定しないため、次の特殊解を指定する 3 番目の条件が必要となります。

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

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

方程式のコード化

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

• x は独立変数です。

• y は従属変数です。

• lambda は固有値を表す未知パラメーターです。

代入 ${\mathit{y}}_1=\mathit{y}$ および ${\mathit{y}}_2={\mathit{y}}^{\prime }$ を使用して、Mathieu 方程式を 1 次系として記述できます。

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

${{\mathit{y}}_2}^{\prime }=-\left(\lambda -2\mathit{q}\cos \left(2\mathit{x}\right)\right){\mathit{y}}_1$.

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

function dydx = mat4ode(x,y,lambda) % equation being solved
dydx = [y(2)
      -(lambda - 2*q*cos(2*x))*y(1)];
end

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

境界条件のコード化

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

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

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

• lambda は固有値を表す未知パラメーターです。

この問題には区間 $\left[0,\pi \right]$ に 3 つの境界条件があります。残差値を計算するには、境界条件を $\mathit{g}\left(\mathit{x},\mathit{y}\right)=0$ の形式にしなければなりません。この形式では、境界条件は次となります。

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

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

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

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

function res = mat4bc(ya,yb,lambda) % boundary conditions
res = [ya(2)
       yb(2)
       ya(1)-1];
end

初期推定の作成

最後に、解の初期推定を作成します。両方の解要素 ${\mathit{y}}_1=\mathit{y}\left(\mathit{x}\right)$ および ${\mathit{y}}_2={\mathit{y}}^{\prime }\left(\mathit{x}\right)$ に対する初期推定と、未知パラメーター $\lambda$ を指定しなければなりません。計算できるのは初期推定に近い固有値および固有関数のみです。計算済みの固有関数が 4 番目の固有値に対応する可能性を高めるには、正確な定性的動作をもつ初期推定を選択する必要があります。

この問題の場合、余弦関数は 3 つの境界条件を満たすため、適切な初期推定を行います。${\mathit{y}}_1$ と ${\mathit{y}}_2$ の推定を返す関数を使用して、$\mathit{y}$ の初期推定をコード化します。

function yinit = mat4init(x) % initial guess function
yinit = [cos(4*x)
        -4*sin(4*x)];
end

区間 $\left[0,\pi \right]$ の 10 点のメッシュ、初期推定関数、$\lambda$ の値に対する 15 の推定を使用して bvpinit を呼び出します。

lambda = 15;
solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);

方程式の求解

ODE 関数、境界条件関数、初期推定を使用して bvp4c を呼び出します。

sol = bvp4c(@mat4ode, @mat4bc, solinit);

パラメーターの値

bvp4c を使って求められる、未知パラメーター $\lambda$ の値を表示します。この値は Mathieu 方程式の 4 番目の固有値 ($\mathit{q}=5$) です。

fprintf('Fourth eigenvalue is approximately %7.3f.\n',...
            sol.parameters)

Fourth eigenvalue is approximately  17.097.

解のプロット

deval を使用して、区間 $\left[0,\pi \right]$ の 100 点で bvp4c によって計算される解を評価します。

xint = linspace(0,pi);
Sxint = deval(sol,xint);

両方の解要素をプロットします。プロットは、4 番目の固有値 ${\lambda }_4=17.097$ に関連付けられている固有関数 (およびその導関数) を示しています。

plot(xint,Sxint)
axis([0 pi -4 4])
title('Eigenfunction of Mathieu''s Equation.') 
xlabel('x')
ylabel('y')
legend('y','y''')

ローカル関数

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

function dydx = mat4ode(x,y,lambda) % equation being solved
q = 5;
dydx = [y(2)
      -(lambda - 2*q*cos(2*x))*y(1)];
end
%-------------------------------------------
function res = mat4bc(ya,yb,lambda) % boundary conditions
res = [ya(2)
       yb(2)
       ya(1)-1];
end
%-------------------------------------------
function yinit = mat4init(x) % initial guess function
yinit = [cos(4*x)
        -4*sin(4*x)];
end
%-------------------------------------------

参考

bvp4c | bvp5c | bvpinit

関連するトピック

• 境界値問題の求解