В этом примере показано, как решить задачу многоточечного граничного значения, где интересующее решение удовлетворяет условиям внутри интервала интегрирования.
Для в ] рассмотрим уравнения
,
) start.
Известные параметры проблемы - , , , и .
Член 1) в уравнении (x) не является гладким x = 1, поэтому задача не может быть решена напрямую. Вместо этого можно разбить задачу на два: один набор в [0,1], а другой набор в [1, λ]. Связь между двумя областями заключается в том, что решения должны быть при x = 1. Решение должно также удовлетворять граничным условиям
= 0,
= 1.
Уравнения для каждой области:
Регион 1:
,
.
Регион 2:
,
.
Точка интерфейса 1 включена в обе области. В этот момент решатель создает как левое, так и правое решения, которые должны быть равными для обеспечения непрерывности решения.
Чтобы решить эту систему уравнений в MATLAB, необходимо кодировать уравнения, граничные условия и начальное предположение, прежде чем вызывать решатель задач граничных значений bvp5c. Требуемые функции можно либо включить в качестве локальных функций в конце файла (как здесь сделано), либо сохранить их как отдельные именованные файлы в каталоге по пути MATLAB.
Уравнения для x) (x) зависят от решаемой области. Для многоточечных граничных задач производная функция должна принимать третий входной аргументregion, который используется для идентификации области, где оценивается производная. Решатель нумерует области слева направо, начиная с 1.
Создание функции для кодирования уравнений с сигнатурой dydx = f(x,y,region,p), где:
x является независимой переменной.
y является зависимой переменной.
dydx(1) дает уравнение для x), иdydx(2) уравнение для x).
region - номер области, в которой вычисляется производная (в этой двухобластной задаче, region является 1 или 2).
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
Примечание.В конце примера все функции включаются как локальные.
Решение двух дифференциальных уравнений первого порядка в двух областях требует четырех граничных условий. Два из этих условий вытекают из первоначальной проблемы:
= 0,
= 0.
Два других условия обеспечивают непрерывность левого и правого решений в точке интерфейса 1:
) = 0,
) = 0.
Для многоточечных BVP, аргументы граничных условий функции YL и YR стать матрицами. В частности, kВ-й колонне YL(:,k) - решение на левой границе kВ-й регион. Аналогично, YR(:,k) - решение на правой границе kВ-й регион.
В этой задаче ) аппроксимируется YL(:,1), в то время как ) аппроксимируется YR(:,end). Непрерывность решения при 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) для обеспечения вектора параметров.
Вычислите решение для нескольких значений, используя каждое решение в качестве начального предположения для следующего. Вычислите для каждого значения, , значение осмолярности (λ). Для каждой итерации цикла сравните вычисленное значение с приблизительным аналитическим решением.
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
Постройте график компонентов решения для ) и x) и вертикальной линии в сопряжения
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 %-------------------------------------------