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

Создайте функцию для кодирования уравнений для $\mathit{f}\mathit{(}x\mathit{,}$y). Эта функция должна иметь подписьdydx = emdenode(x,y), где:

Эти входные данные автоматически передаются функции решателем, но имена переменных определяют способ кодирования уравнений. В этом случае:

function dydx = emdenode(x,y)
dydx = [y(2) 
       -y(1)^5];
end

Термин, который содержит S передается решателю с помощью опций, поэтому этот член не включается в функцию.

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

Теперь запишите функцию, которая возвращает остаточное значение граничных условий в граничных точках. Эта функция должна иметь подпись res = emdenbc(ya,yb), где:

Для этой задачи одно из граничных условий - для ${\mathit{}}_{\mathrm{y1}}$, а другое - для ${\mathit{}}_{\mathrm{y2}}$. Для вычисления остаточных значений необходимо поместить граничные условия в форму $\mathit{g}\mathit{(}x\mathit{,}y\mathrm{)}$= 0.

В этой форме граничными условиями являются

${\mathit{}}_{\mathrm{y2}}\mathrm{}(0)\mathrm{=0}$,

${\mathit{}}_{\mathrm{y1}}\mathrm{}(1)\sqrt{\mathrm{}}\mathrm{}\mathrm{-3/2=0}$.

Затем выполняется соответствующая функция.

function res = emdenbc(ya,yb)
res = [ya(2)
       yb(1) - sqrt(3)/2];
end

Создать начальное предположение

Наконец, создайте первоначальное предположение о решении. Для этой задачи используйте постоянное начальное предположение, удовлетворяющее граничным условиям, и простую сетку из пяти точек между 0 и 1. Использование многих точек сетки не требуется, поскольку решатель BVP адаптирует эти точки в процессе решения.

${\mathit{}}_{\mathrm{}}\sqrt{\mathrm{}}\mathrm{y1=3/2}$,

${\mathit{}}_{\mathrm{}}\mathrm{y2=0}$.

guess = [sqrt(3)/2; 0];
xmesh = linspace(0,1,5);
solinit = bvpinit(xmesh, guess);

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

Создание матрицы для S и передать его в bvpset как значение 'SingularTerm' вариант. Наконец, позвоните bvp4c с функцией ОДУ, функцией граничного условия, начальным предположением и структурой опций.

S = [0 0; 0 -2];
options = bvpset('SingularTerm',S);
sol = bvp4c(@emdenode, @emdenbc, solinit, options);

Решение для построения графика

Вычислите значение аналитического решения в $[\mathrm{}\mathrm{0,1}$].

x = linspace(0,1);
truy = 1 ./ sqrt(1 + (x.^2)/3);

Постройте график аналитического решения и решения, рассчитанного по bvp4c для сравнения.

plot(x,truy,sol.x,sol.y(1,:),'ro');
title('Emden Problem -- BVP with Singular Term.')
legend('Analytical','Computed');
xlabel('x');
ylabel('solution y');

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

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

function dydx = emdenode(x,y) % equation being solved 
dydx = [y(2) 
       -y(1)^5];
end
%-------------------------------------------
function res = emdenbc(ya,yb) % boundary conditions
res = [ya(2)
       yb(1) - sqrt(3)/2];
end
%-------------------------------------------