Кодовые уравнения

Уравнения для ${\mathit{v}}^{\prime }\left(\mathit{x}\right)$ и ${\mathit{C}}^{\prime }\left(\mathit{x}\right)$ зависят от решаемой области. Для многоточечных граничных задач значения производная функция должна принять третий входной параметр region, который используется для идентификации области, в которой оценивается производная. Решатель нумерует области слева направо, начиная с 1.

Создайте функцию, чтобы кодировать уравнения сигнатурой dydx = f(x,y,region,p), где:

Используйте оператор 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

Примечание.Все функции включены в качестве локальных функций в конце примера.

Граничные условия кода

Решение двух дифференциальных уравнений первого порядка в двух областях требует четырех граничных условий. Два из этих условий являются результатом первоначальной задачи:

$\mathit{v}\left(0\right)=0$,

$\mathit{C}\left(\lambda \right)-1=0$.

Другие два условия обеспечивают непрерывность левых и правых решений в точке интерфейса $\mathit{x}=1$:

${\mathit{v}}_{\mathit{L}}\left(1\right)-{\mathit{v}}_{\mathit{R}}\left(1\right)=0$,

${\mathit{C}}_{\mathit{L}}\left(1\right)-{\mathit{C}}_{\mathit{R}}\left(1\right)=0$.

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

В этой задаче, $\mathit{y}\left(0\right)$ аппроксимируется YL(:,1), в то время как $\mathit{y}\left(\lambda \right)$ аппроксимируется YR(:,end). Непрерывность решения при $\mathit{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. Предоставьте функцию для bvp5c с синтаксисом @(x,y,r) f(x,y,r,p) для предоставления вектора параметров.

Вычислим решение для нескольких значений $\kappa$, использование каждого решения в качестве начального предположения для следующего. Для каждого значения $\kappa$, вычислите значение осмолярности $\mathrm{Os}=\frac{1}{\mathit{v}\left(\lambda \right)}$. Для каждой итерации цикла сравните вычисленное значение с приблизительным аналитическим решением.

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 

Решение для построения графика

Постройте график компонентов решения для $\mathit{v}\left(\mathit{x}\right)$ и $\mathit{C}\left(\mathit{x}\right)$ и вертикальную линию в точке интерфейса $\mathit{x}=1$. Отображаемое решение для $\kappa =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)')

Локальные функции

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