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

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

${{\mathit{}}_{\mathrm{}}}^{y1\; }={\mathit{}}_{\mathrm{y2}}$,

${{\mathit{}}_{\mathrm{}}}^{}\frac{\mathit{}}{\mathit{}}{\mathit{y2\prime =-xey}}^{}\prime {}^{\mathrm{-}}\mathrm{}\mathrm{ls2cos}\mathit{}(\frac{\mathrm{¼\; x}\mathit{)}}{\mathit{}}\mathrm{-}\mathit{}āxesin$(¼ x).

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

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

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

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

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

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

Кодовые якобинцы

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

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

${\mathit{}}_{\mathrm{}}\frac{\mathit{}}{\mathit{}}\begin{array}{cc}\frac{{\mathit{}}_{\mathrm{}}}{{\mathit{}}_{\mathrm{}}}& \frac{{\mathit{}}_{\mathrm{}}}{{\mathit{}}_{\mathrm{}}}\\ \frac{{\mathit{}}_{\mathrm{}}}{{\mathit{}}_{\mathrm{}}}& \frac{{\mathit{}}_{\mathrm{}}}{{\mathit{}}_{\mathrm{}}}\end{array}\begin{array}{cc}\mathrm{}& \mathrm{}\\ \mathrm{}& \frac{\mathit{}}{\mathit{JODE=\partial f\partial y=[\partial f1\partial y1\partial f1\partial y2\partial f2\partial y1\partial f2\partial y2]=[010-xe}}\end{array}$].

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

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

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

${\mathit{}}_{\mathit{Jy}\mathit{(}a})\begin{array}{cc}\mathrm{=}& \mathrm{}\\ \mathrm{[}& \mathrm{}\end{array}$1000]${\mathit{,}}_{\mathit{}\mathit{Jy}}(b\begin{array}{cc}\mathrm{)}& \mathrm{}\\ \mathrm{=}& \mathrm{}\end{array}[$0010].

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

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

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

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

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

Уравнение решения

Если попытаться решить уравнение непосредственно с помощью $\mathit{e}{\mathrm{=}}^{\mathrm{}}$10-4, то решатель не сможет преодолеть плохое кондиционирование задачи вблизи $\mathit{}\mathrm{}$точки перехода x = 0. Вместо этого для получения решения $\mathit{}{\mathrm{для}}^{e\mathrm{}}$= 10-4 в этом примере используется продолжение путем решения последовательности ${\mathrm{задач}}^{\mathrm{}}$для ${\mathrm{}}^{\mathrm{10-2}}$, 10-3 ${\mathrm{}}^{\mathrm{}}$и 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 для сетки $\mathit{x}$ и раствора $\mathit{y}\mathit{(}x$). При продолжении решатель может обрабатывать разрывы при $\mathit{}x\mathrm{}$= 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
%-------------------------------------------