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

Используя замены ${\mathit{y}}_1=\mathit{y}$ и ${\mathit{y}}_2={\mathit{y}}^{\prime }$, можно переписать уравнение как систему уравнений первого порядка

${{\mathit{y}}_1}^{\prime }={\mathit{y}}_2$,

${{\mathit{y}}_2}^{\prime }=-\frac{\mathit{x}}{\mathit{e}}{\mathit{y}}^{\prime }-{\pi }^2\cos \left(\pi \mathit{x}\right)-\frac{\pi \mathit{x}}{\mathit{e}}\sin \left(\pi \mathit{x}\right)$.

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

Сделайте функцию векторизованной, так, чтобы 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 требует, чтобы граничные условия были в форме $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$. В этой форме граничные условия:

$\mathit{y}\left(-1\right)+2=0$,

$\mathit{y}\left(1\right)=0$.

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

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

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

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

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

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

${\mathit{J}}_{\mathrm{ODE}}=\frac{\partial \mathit{f}}{\partial \mathit{y}}=\left[\begin{array}{cc}\frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_2}\\ \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_2}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& -\frac{\mathit{x}}{\mathit{e}}\end{array}\right]$.

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

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

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

${\mathit{J}}_{\mathit{y}\left(\mathit{a}\right)}=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$ , ${\mathit{J}}_{\mathit{y}\left(\mathit{b}\right)}=\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right]$.

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

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

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

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

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

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

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

Поскольку значение якобиана зависит от значения $\mathit{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 $\mathit{x}$ и решение $\mathit{y}\left(\mathit{x}\right)$. С продолжением решатель может обработать разрыв в $\mathit{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
%-------------------------------------------