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

Открыть Live Script

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

После сокращения PDE модели с помощью симметрии уравнение становится ОДУ второго порядка, заданной на интервале $[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{1 + \frac{x^{2}}{3}}]}^{- 1}$ .

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

$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

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

Наконец, создайте начальное предположение о решении. Для этой задачи используйте постоянное начальное предположение, которое удовлетворяет граничным условиям, и простой mesh из пяти точек между 0 и 1. Использование многих сетчатый не требуется, поскольку решатель 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');

Figure contains an axes. The axes with title Emden Problem -- BVP with Singular Term. contains 2 objects of type line. These objects represent Analytical, Computed.

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

Здесь перечислены локальные вспомогательные функции, которые решатель 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
%-------------------------------------------

См. также

bvp4c | bvp5c | bvpinit | bvpset

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