Решите BVP Используя продолжение

В этом примере показано, как решить численно трудную краевую задачу с помощью продолжения, которое эффективно разбивает проблему в последовательность более простых проблем.

Для $0 < e ≪ 1$ , рассмотрите дифференциальное уравнение

$e {y^{'}}^{'} + {x, y}^{'} = - e π^{2} потому что (π x) - π xsin (π x)$ .

Проблема создана на интервале $[- 1, 1]$ и подчиняется граничным условиям

$y (- 1) = - 2$ ,

$y (1) = 0$ .

Когда $e = 10^{- 4}$ , решение уравнения подвергается быстрому переходу рядом $x = 0$ , таким образом, это затрудняет, чтобы решить численно. Вместо этого этот пример использует продолжение, чтобы выполнить итерации через несколько значений $e$ до $e = 10^{- 4}$ . Промежуточные решения каждый используются в качестве исходного предположения для следующей проблемы.

Чтобы решить эту систему уравнений в MATLAB, необходимо закодировать уравнения, граничные условия и исходное предположение прежде, чем вызвать решатель для краевой задачи bvp4c. Вы любой может включать необходимые функции как локальные функции в конце файла (как сделано здесь), или можно сохранить их как отдельные, именованные файлы в директории на пути MATLAB.

Уравнение кода

Используя замены $y_{1} = y$ и $y_{2} = y^{'}$ , можно переписать уравнение как систему уравнений первого порядка

${y_{1}}^{'} = y_{2}$ ,

${y_{2}}^{'} = - \frac{x}{e} y^{'} - π^{2} потому что (π x) - \frac{π x}{e} \sin (π x)$ .

Запишите функцию, чтобы закодировать уравнения с подписью dydx = shockode(x,y), где:

x независимая переменная.
y зависимая переменная.
dydx(1) дает уравнение для ${y_{1}}^{'}$ , и dydx(2) дает уравнение для ${y_{2}}^{'}$ .

Сделайте функцию векторизованной, так, чтобы shockode([x1 x2 ...],[y1 y2 ...]) возвращает [shockode(x1,y1) shockode(x2,y2) ...]. Этот подход улучшает производительность решателя.

Соответствующая функция

function dydx = shockode(x,y)
pix = pi*x;
dydx = [y(2,:)
       -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)];
end

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

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

Решатель BVP требует, чтобы граничные условия были в форме $g (y (a), y (b)) = 0$ . В этой форме граничные условия:

$y (- 1) + 2 = 0$ ,

$y (1) = 0$ .

Запишите функцию, чтобы закодировать граничные условия с подписью res = shockbc(ya,yb), где:

ya значение граничного условия в начале интервала $[a, b]$ .
yb значение граничного условия в конце интервала $[a, b]$ .

Соответствующая функция

function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end

Якобианы кода

Аналитические Якобианы для функции ОДУ и граничных условий могут быть вычислены легко в этой проблеме. Обеспечение Якобианов делает решатель более эффективным, поскольку решатель больше не должен аппроксимировать их конечными разностями.

Для функции ОДУ якобиевская матрица

$J_{ОДУ} = \frac{\partial f}{\partial y} = [\begin{matrix} \frac{\partial f_{1}}{\partial y_{1}} & \frac{\partial f_{1}}{\partial y_{2}} \\ \frac{\partial f_{2}}{\partial y_{1}} & \frac{\partial f_{2}}{\partial y_{2}} \end{matrix}] = [\begin{matrix} 0 & 1 \\ 0 & - \frac{x}{e} \end{matrix}]$ .

Соответствующая функция

function jac = shockjac(x,y,e)
jac = [0   1
       0  -x/e];
end

Точно так же для граничных условий, якобиевские матрицы

$J_{y (a)} = [\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}]$ , $J_{y (b)} = [\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix}]$ .

Соответствующая функция

function [dBCdya,dBCdyb] = shockbcjac(ya,yb)
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end

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

Используйте постоянное предположение в решении на сетке пяти точек в $[- 1, 1]$ .

sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);

Решите уравнение

При попытке решить уравнение непосредственно с помощью $e = 10^{- 4}$ , затем решатель не может преодолеть плохое создание условий проблемы около $x = 0$ точка перехода. Вместо этого чтобы получить решение для $e = 10^{- 4}$ , этот пример использует продолжение путем решения последовательности проблем для $10^{- 2}$ , $10^{- 3}$ , и $10^{- 4}$ . Выход от решателя в каждой итерации действует как предположение для решения в следующей итерации (это то, почему переменная для исходного предположения от bvpinit sol, и выход от решателя также называют sol).

Поскольку значение якобиана зависит от значения $e$ , установите опции в цикле, задающем shockjac и shockbcjac функции для Якобианов. Кроме того, включите векторизацию начиная с shockode закодирован, чтобы обработать векторы значений.

e = 0.1;
for i = 2:4
   e = e/10;
   options = bvpset('FJacobian',@(x,y) shockjac(x,y,e),'BCJacobian',@shockbcjac,'Vectorized','on');
   sol = bvp4c(@(x,y) shockode(x,y,e),@shockbc, sol, options);
end

Постройте решение

Постройте результаты bvp4c для mesh $x$ и решение $y (x)$ . С продолжением решатель может обработать разрыв в $x = 0$ .

plot(sol.x,sol.y(1,:),'-o');
axis([-1 1 -2.2 2.2]);
title(['There Is a Shock at x = 0 When e =' sprintf('%.e',e) '.']);
xlabel('x');
ylabel('solution y');

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

Перечисленный здесь локальные функции что решатель BVP bvp4c вызовы, чтобы вычислить решение. В качестве альтернативы можно сохранить эти функции как их собственные файлы в директории на пути MATLAB.

function dydx = shockode(x,y,e) % equation to solve
pix = pi*x;
dydx = [y(2,:)
       -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)];
end
%-------------------------------------------
function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end
%-------------------------------------------
function jac = shockjac(x,y,e) % jacobian of shockode
jac = [0   1
       0  -x/e];
end
%-------------------------------------------
function [dBCdya,dBCdyb] = shockbcjac(ya,yb) % jacobian of shockbc
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end
%-------------------------------------------

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

bvp4c | bvpinit | bvpset

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