Решите BVP с несколькими граничными условиями

Этот пример показывает, как решить многоточечную краевую задачу, где решение интереса удовлетворяет условия в интервале интегрирования.

Для $x$ \in $[0, λ]$ , рассмотрите уравнения

$v^{'} = \frac{C - 1}{n}$ ,

$C^{'} = \frac{vC - \min (x, 1)}{η}$ .

Известные параметры проблемы $n$ , $κ$ , $λ > 1$ , и $η = \frac{λ^{2}}{n \cdot κ^{2}}$ .

Термин $\min (x, 1)$ в уравнении для $C^{'} (x)$ не сглаженно в $x = 1$ , таким образом, проблема не может быть решена непосредственно. Вместо этого можно повредить проблему в два: один набор в интервале $[0, 1]$ , и другой набор в интервале $[1, λ]$ . Связь между этими двумя областями состоит в том, что решения должны быть непрерывными в $x = 1$ . Решение должно также удовлетворить граничные условия

$v (0) = 0$ ,

$C (λ) = 1$ .

Уравнения для каждой области

Область 1: $0 \leq x \leq 1$

$v^{'} = \frac{C - 1}{n}$ ,

$C^{'} = \frac{vC - x}{η}$ .

Область 2: $1 \leq x \leq λ$

$v^{'} = \frac{C - 1}{n}$ ,

$C^{'} = \frac{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 является количеством области, где производная вычисляется (в этой проблеме 2D области, region является 1 или 2).
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$ :

$v_{L} (1) - v_{R} (1) = 0$ ,

$C_{L} (1) - C_{R} (1) = 0$ .

Для многоточечного BVPs функционируют аргументы граничных условий, YL и YR становятся матрицами. В частности, k th столбец YL(:,k) является решением на левом контуре k th область. Точно так же YR(:,k) является решением на правильном контуре k th область.

В этой проблеме, $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

Сформируйте исходное предположение

Для многоточечного BVPs граничные условия автоматически применяются вначале и конец интервала интегрирования. Однако необходимо задать двойные записи в 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. Предоставьте функцию bvp5c с синтаксисом @(x,y,r) f(x,y,r,p), чтобы обеспечить вектор параметров.

Вычислите решение для нескольких значений $κ$ , использование каждого решения как исходное предположение для следующего. Для каждого значения $κ$ , вычислите значение osmolarity $OS = \frac{1}{v (λ)}$ . Для каждой итерации цикла сравните вычисленное значение с аппроксимированным аналитическим решением.

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)')

Localfunctions

Перечисленный здесь локальные функции помощника, которые решатель 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
%-------------------------------------------

Смотрите также

bvp4c | bvp5c | bvpinit

Документация