Main Content

求解具有多边界条件的 BVP

以下示例说明如何求解多点边界值问题,其中关注的解满足积分区间内的条件。

对于 [0,λ] 中的 x,考虑以下方程

v=C-1n,

C=vC-min(x,1)η.

问题的已知参数是 nκλ>1η=λ2nκ2

C(x) 的方程中的项 min(x,1)x=1 处不平滑,因此该问题不能直接求解。在这种情况下,您可以将问题分成两部分:一部分在区间 [0,1] 内,另一部分在区间 [1,λ] 内。这两个区域之间的联系是在 x=1 处的解必须为连续的。解还必须满足边界条件

v(0)=0,

C(λ)=1.

每个区域的方程如下

区域 1: 0x1

v=C-1n,

C=vC-xη.

区域 2: 1xλ

v=C-1n,

C=vC-1η.

交界点 x=1 同时包含在这两个区域中。在此交界点上,求解器会产生解和解,这两个解必须相等,以确保解的连续性。

要在 MATLAB® 中对此方程组求解,您需要先编写方程组、边界条件和初始估计值的代码,然后再调用边界值问题求解器 bvp5c。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写方程代码

v(x)C(x) 的方程取决于正在求解的区域。对于多点边界值问题,导数函数必须接受第三个输入参量 region,该参量用于标识正在计算导数的区域。求解器从左到右对区域进行编号,从 1 开始。

创建一个函数以使用签名 dydx = f(x,y,region,p) 来编写方程代码,其中:

  • x 是自变量。

  • y 是因变量。

  • dydx(1) 给出 v(x) 的方程,dydx(2) 给出 C(x) 的方程。

  • region 是计算导数的区域编号(本例中问题分为两个区域,region12)。

  • p 是向量,包含常量参数 [n, κ, λ, η] 的值。

根据要求解的具体区域,使用 switch 语句返回不同方程。函数是

function dydx = f(x,y,region,p) % equations being solved
n = p(1);
eta = p(4);

dydx = zeros(2,1);
dydx(1) = (y(2) - 1)/n;

switch region
    case 1    % x in [0 1]
        dydx(2) = (y(1)*y(2) - x)/eta;
    case 2    % x in [1 lambda]
        dydx(2) = (y(1)*y(2) - 1)/eta;
end
end

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

编写边界条件代码

在两个区域中求解两个一阶微分方程需要四个边界条件。这些条件中有两个来自原始问题:

v(0)=0,

C(λ)-1=0.

另外两个条件强制交界点 x=1 处的左解和右解具备连续性:

vL(1)-vR(1)=0,

CL(1)-CR(1)=0.

对于多点 BVP,边界条件函数 YLYR 的参量会是矩阵。具体来说,第 kYL(:,k) 是第 k 个区域左边界的解。YR(:,k) 则是第 k 个区域右边界的解。

在此问题中,y(0) 通过 YL(:,1) 来逼近,而 y(λ) 通过 YR(:,end) 来逼近。解在 x=1 处的连续性要求 YR(:,1) = YL(:,2)

用于编写四个边界条件的残差值计算代码的函数是

function res = bc(YL,YR)
res = [YL(1,1)               % v(0) = 0
       YR(1,1) - YL(1,2)     % Continuity of v(x) at x=1
       YR(2,1) - YL(2,2)     % Continuity of C(x) at x=1
       YR(2,end) - 1];       % C(lambda) = 1
end

获取初始估计值

对于多点 BVP,边界条件自动应用于积分区间的开始处和结束处。但是,您必须在 xmesh 中为其他交界点分别指定双重项。满足边界条件的简单估计值是常量估计值 y = [1; 1]

xc = 1;
xmesh = [0 0.25 0.5 0.75 xc xc 1.25 1.5 1.75 2];
yinit = [1; 1];
sol = bvpinit(xmesh,yinit);

求解方程

定义常量参数的值,并将其放入向量 p 中。使用语法 @(x,y,r) f(x,y,r,p)bvp5c 提供函数,以提供参数向量。

计算 κ 的几个值的解,其中使用每个解作为下一个解的初始估计值。对于 κ 的每个值,计算渗透性 Os=1v(λ) 的值。对于循环的每次迭代,将计算值与近似解析解进行比较。

lambda = 2;
n = 5e-2;
for kappa = 2:5
   eta = lambda^2/(n*kappa^2);
   p = [n kappa lambda eta];
   sol = bvp5c(@(x,y,r) f(x,y,r,p), @bc, sol);
   K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
   approx = 1/(1 - K2);
   computed = 1/sol.y(1,end);
   fprintf('  %2i    %10.3f     %10.3f \n',kappa,computed,approx);
end
   2         1.462          1.454 
   3         1.172          1.164 
   4         1.078          1.071 
   5         1.039          1.034 

对解进行绘图

绘制 v(x)C(x) 的解分量,以及在交界点 x=1 处的垂直线。显示的 κ=5 的解是循环的最后一次迭代的结果。

plot(sol.x,sol.y(1,:),'--o',sol.x,sol.y(2,:),'--o')
line([1 1], [0 2], 'Color', 'k')
legend('v(x)','C(x)')
title('A Three-Point BVP Solved with bvp5c')
xlabel({'x', '\lambda = 2, \kappa = 5'})
ylabel('v(x) and C(x)')

Figure contains an axes object. The axes object with title A Three-Point BVP Solved with bvp5c, xlabel x lambda blank = blank 2 , blank kappa blank = blank 5, ylabel v(x) and C(x) contains 3 objects of type line. These objects represent v(x), C(x).

局部函数

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

function dydx = f(x,y,region,p) % equations being solved
n = p(1);
eta = p(4);

dydx = zeros(2,1);
dydx(1) = (y(2) - 1)/n;

switch region
    case 1    % x in [0 1]
        dydx(2) = (y(1)*y(2) - x)/eta;
    case 2    % x in [1 lambda]
        dydx(2) = (y(1)*y(2) - 1)/eta;
end
end
%-------------------------------------------
function res = bc(YL,YR) % boundary conditions
res = [YL(1,1)               % v(0) = 0
       YR(1,1) - YL(1,2)     % Continuity of v(x) at x=1
       YR(2,1) - YL(2,2)     % Continuity of C(x) at x=1
       YR(2,end) - 1];       % C(lambda) = 1
end
%-------------------------------------------

另请参阅

| |

相关主题