MATLAB 数学应用 微分方程 边界值问题 使用延拓验证BVP一致性

本文讲述了如何使用延拓将 BVP 的一个解逐渐扩展到更大的区间。

Falkner-Skan 边界值问题源于为平板粘性不可压缩层流问题求取相似解的过程。示例方程是

$f^{\prime {}^{\prime {}^{\prime }}}+f{f}^{\prime {}^{\prime }}+\beta (1-f^{\prime 2})=0$。

此问题位于无限区间 [0,∞] 和 β=0.5 上，并且需要满足边界条件

f(0)=0，

f’(0)=0，

f′(∞)=1。

在无限区间上无法求解 BVP，在非常大的有限区间上求解 BVP 也不切实际。在这种情况下，此示例转而求解位于较小区间 [0,a] 上的一系列问题，从而验证解具有与 a→∞ 一致的行为。这种做法称为延拓，它是将一个问题分解成多个更简单的问题，将每个小问题所反馈的解作为下一个小问题的初始估计值。

要在 MATLAB 中对此方程组求解，您需要先编写方程组、边界条件和选项的代码，然后再调用边界值问题求解器 bvp4c。您可以将所需的函数作为局部函数包含在文件末尾，或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写方程代码

创建一个函数以编写方程代码。此函数应具有签名 dfdx = fsode(x,f)，其中：

x 是自变量。

f 是因变量。

您可以使用代换法 $f_1=f$、$f_2=f\prime$ 和 $f_3=f{\prime }^{\prime }$ 将三阶方程重写为一阶方程组。方程变为

$f_1^{\prime }=f_2$，

$f_2^{\prime }=f_3$，

$f_3^{\prime }=-f_1{f}_3-\beta (1-f_2^2)$。

对应的函数是

function dfdeta = fsode(x,f)
b = 0.5;
dfdeta = [ f(2)
           f(3)
          -f(1)*f(3) - b*(1 - f(2)^2) ];
end

注意：所有函数都作为局部函数包含在示例的末尾。

编写边界条件代码

现在，编写一个函数，该函数返回在边界点处的边界条件的残差值。此函数应具有签名 res = fsbc(f0,finf)，其中：

f0 是在区间的开始处的边界条件的值。

finf 是在区间的结束处的边界条件的值。

要计算残差值，您需要将边界条件设置为 g(x,y)=0 形式。在此形式中，边界条件是

f(0)=0，

f′(0)=0，

f′(∞)−1=0。

对应的函数是

function res = fsbc(f0,finf)
res = [f0(1)
       f0(2)
       finf(2) - 1];
end

创建初始估计值

最后，您必须提供解的初始估计值。具有五个点的粗网格以及满足边界条件的常量估计值便可在区间 [0,3] 上进行收敛。变量 infinity 表示积分区间的右侧极限。随着 infinity 的值在随后的迭代中从 3 增加到其最大值 6，每个先前迭代给出的解充当下一次迭代的初始估计值。

infinity = 3;
maxinfinity = 6;
solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);

求解方程并绘制解

在初始区间 [0, 3] 中求解问题。绘制解，并将 $f^{{}^{\prime {}^{\prime }}}(0)$ 的值与解析值 [1] 进行比较。

sol = bvp4c(@fsode,@fsbc,solinit);
x = sol.x;
f = sol.y;
plot(x,f(2,:),x(end),f(2,end),'o');
axis([0 maxinfinity 0 1.4]);
title('Falkner-Skan Equation, Positive Wall Shear, \beta = 0.5.')
xlabel('x')
ylabel('df/dx')
hold on

fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n')
Cebeci & Keller report that f''(0) = 0.92768.
fprintf('Value computed using infinity = %g is %7.5f.\n', ...
         infinity,f(3,1))
Value computed using infinity = 3 is 0.92915.

现在，通过为每次迭代增加 infinity 的值，在逐步增大的区间上求解问题。bvpinit 函数将每个解外推至新区间，以作为 infinity 的下一个值的初始估计值。每次迭代都会打印 $f^{{}^{\prime {}^{\prime }}}(0)$ 的计算值，并将解的绘图叠加到之前的解上。当 infinity = 6 时，解的一致行为变得明显，$f^{{}^{\prime {}^{\prime }}}(0)$ 的值非常接近预测值。

for Bnew = infinity+1:maxinfinity
  solinit = bvpinit(sol,[0 Bnew]); % Extend solution to new interval
  sol = bvp4c(@fsode,@fsbc,solinit);
  x = sol.x;
  f = sol.y;
  
  fprintf('Value computed using infinity = %g is %7.5f.\n', ...
           Bnew,f(3,1))
  plot(x,f(2,:),x(end),f(2,end),'o');
  drawnow
end
Value computed using infinity = 4 is 0.92774.
Value computed using infinity = 5 is 0.92770.
Value computed using infinity = 6 is 0.92770.
hold off

局部函数

此处列出了 BVP 求解器 bvp4c 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。

function dfdeta = fsode(x,f) % equation being solved
dfdeta = [ f(2)
           f(3)
          -f(1)*f(3) - 0.5*(1 - f(2)^2) ];
end 
%-------------------------------------------
function res = fsbc(f0,finf) % boundary conditions
res = [f0(1)
       f0(2)
       finf(2) - 1];
end
%-------------------------------------------