本文讲述了如何使用延拓将 BVP 的一个解逐渐扩展到更大的区间。
Falkner-Skan 边界值问题源于为平板粘性不可压缩层流问题求取相似解的过程。示例方程是
f ′ ′ ′ + f f ′ ′ + β ( 1 − f ′ 2 ) = 0 f'^{'^{'}} + f f'^{'} + β(1 - f'^2) = 0 f′′′+ff′′+β(1−f′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_1 = f f1=f、 f 2 = f ′ f_2 = f′ f2=f′ 和 f 3 = f ′ ′ f_3 = f′^{′} f3=f′′ 将三阶方程重写为一阶方程组。方程变为
f 1 ′ = f 2 f'_1 = f_2 f1′=f2，
f 2 ′ = f 3 f'_2 = f_3 f2′=f3，
f 3 ′ = − f 1 f 3 − β ( 1 − f 2 2 ) f'_3 = -f_1f_3 - β(1 - f_2^{2}) f3′=−f1f3−β(1−f22)。
对应的函数是
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 ′ ′ ( 0 ) f^{'^{'}}(0) f′′(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 ′ ′ ( 0 ) f^{'^{'}}(0) f′′(0) 的计算值，并将解的绘图叠加到之前的解上。当 infinity = 6 时，解的一致行为变得明显， f ′ ′ ( 0 ) f^{'^{'}}(0) f′′(0) 的值非常接近预测值。
for Bnew = infinity+1:maxinfinitysolinit = bvpinit(sol,[0 Bnew]); % Extend solution to new intervalsol = 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');drawnowendValue 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 solveddfdeta = [ f(2)f(3)-f(1)*f(3) - 0.5*(1 - f(2)^2) ];end %-------------------------------------------function res = fsbc(f0,finf) % boundary conditionsres = [f0(1)f0(2)finf(2) - 1];end%-------------------------------------------
