Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

使用延拓验证 BVP 一致性

以下示例说明如何使用延拓将 BVP 的一个解逐渐扩展到更大的区间。

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

f+f f+β(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 是因变量。

您可以使用代换法 f1=ff2=ff3=f 将三阶方程重写为一阶方程组。方程变为

f1=f2

f2=f3

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) 的值与解析值 [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

Figure contains an axes. The axes with title Falkner-Skan Equation, Positive Wall Shear, \beta = 0.5. contains 2 objects of type line.

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) 的计算值,并将解的绘图叠加到之前的解上。当 infinity = 6 时,解的一致行为变得明显,f(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

Figure contains an axes. The axes with title Falkner-Skan Equation, Positive Wall Shear, \beta = 0.5. contains 8 objects of type line.

局部函数

此处列出了 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
%-------------------------------------------

参考

[1] Cebeci, T. and H. B. Keller."Shooting and Parallel Shooting Methods for Solving the Falkner-Skan Boundary-layer Equation."J. Comp.Phys., Vol. 7, 1971, pp. 289-300.

另请参阅

| |

相关主题