求解具有奇异项的 BVP

以下示例说明如何求解埃姆登方程，埃姆登方程是一个具有奇异项的边界值问题，源于对气体球体建模的过程。

在使用对称性法简化模型的 PDE 后，该方程变为在区间 $\left[0,1\right]$ 上定义的二阶 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 }$ 将埃姆登方程重写为一阶方程组

${{\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 路径上的目录中。

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

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

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

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

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

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

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

相关主题

• 求解边界值问题