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

Использование замен ${\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

Код якобиан

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

Для функции ODE, матрица Якобия является

${\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

Начальное предположение формы

Используйте постоянное предположение для решения на mesh из пяти точек в $\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
%-------------------------------------------