Решите BVP с сингулярным термином

В этом примере показано, как решить уравнение Эмдена, которое является краевой задачей с сингулярным термином, который возникает в моделировании сферического тела газа.

После сокращения УЧП модели с помощью симметрии уравнение становится ОДУ второго порядка, заданным на интервале $[0, 1]$ ,

$y^{''} + \frac{2}{x} y^{'} + y^{5} = 0$ .

В $x = 0$ , $(2 / x)$ термин сингулярен, но симметрия подразумевает граничное условие $y^{'} (0) = 0$ . С этим граничным условием, термином $(2 / x) y^{'}$ четко определен как $x \to 0$ . Для граничного условия $y (1) = \sqrt{3} / 2$ , BVP имеет аналитическое решение, которое можно сравнить с числовым решением, вычисленным в MATLAB®,

$y (x) = {[\sqrt{\frac{1 + x^{2}}{3}}]}^{- 1}$ .

Решатель BVP bvp4c может решить сингулярные BVPs, которые имеют форму

$y^{'} = S \frac{y}{x} + f (x, y)$ .

Матрица $S$ должно быть постоянным и граничные условия в $x = 0$ должно быть сопоставимо с необходимым условием $S \cdot y (0) = 0$ . Используйте 'SingularTerm' опция bvpset передать S матрица к решателю.

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

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

${y_{2}}^{'} = - \frac{2}{x} y_{2} - y_{1}^{5}$ .

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

$y_{2} (0) = 0$ ,

$y_{1} (1) = \sqrt{3} / 2$ .

В необходимой матричной форме система

$[\begin{matrix} {y_{1}}^{'} \\ {y_{2}}^{'} \end{matrix}] = \frac{1}{x} [\begin{matrix} 0 & 0 \\ 0 & - 2 \end{matrix}] [\begin{matrix} y_{1} \\ y_{2} \end{matrix}] + [\begin{matrix} y_{2} \\ - y_{1}^{5} \end{matrix}]$ .

В матричной форме это ясно это $S = [\begin{matrix} 0 & 0 \\ 0 & - 2 \end{matrix}]$ и $f (x, y) = [\begin{matrix} y_{2} \\ - y_{1}^{5} \end{matrix}]$ .

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

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

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

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

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

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

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

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

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

ya значение граничного условия в начале интервала.
yb значение граничного условия в конце интервала.

Для этой проблемы одно из граничных условий для $y_{1}$ , и другой для $y_{2}$ . Чтобы вычислить остаточные значения, необходимо поместить граничные условия в форму $g (x, y) = 0$ .

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

$y_{2} (0) = 0$ ,

$y_{1} (1) - \sqrt{3} / 2 = 0$ .

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

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

Создайте исходное предположение

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

$y_{1} = \sqrt{3} / 2$ ,

$y_{2} = 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);

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

Вычислите значение аналитического решения в $[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
%-------------------------------------------

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