Этот пример показывает, как решить уравнение Эмдена, которое является краевой задачей с сингулярным термином, который возникает в моделировании сферического тела газа.
После сокращения УЧП модели с помощью симметрии уравнение становится ОДУ второго порядка, заданным на интервале ,
.
В , термин сингулярен, но симметрия подразумевает граничное условие . С этим граничным условием, термином четко определен как . Для граничного условия , BVP имеет аналитическое решение, которое можно сравнить с числовым решением, вычисленным в MATLAB®,
.
bvp4c
решателя BVP может решить сингулярные BVPs, которые имеют форму
.
Матрица должно быть постоянным и граничные условия в должно быть сопоставимо с необходимым условием . Используйте опцию 'SingularTerm'
bvpset
, чтобы передать матрицу S
решателю.
Можно переписать уравнение Эмдена как систему использования уравнений первого порядка и как
,
.
Граничные условия становятся
,
.
В необходимой матричной форме система
.
В матричной форме это ясно это и .
Чтобы решить эту систему уравнений в MATLAB, необходимо закодировать уравнения, граничные условия и опции прежде, чем вызвать решатель для краевой задачи bvp4c
. Вы любой может включать необходимые функции как локальные функции в конце файла (как сделано здесь) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.
Создайте функцию, чтобы закодировать уравнения для . Эта функция должна иметь подпись dydx = emdenode(x,y)
, где:
x
является независимой переменной.
y
является зависимой переменной.
dydx(1)
дает уравнение для , и dydx(2)
уравнение для .
Эти входные параметры автоматически передаются функции решателем, но имена переменных определяют, как вы кодируете уравнения. В этом случае:
function dydx = emdenode(x,y) dydx = [y(2) -y(1)^5]; end
Термин, который содержит S
, передается решателю с помощью опций, так, чтобы термин не был включен в функцию.
Теперь, запишите функцию, которая возвращает остаточное значение граничных условий в граничных точках. Эта функция должна иметь подпись res = emdenbc(ya,yb)
, где:
ya
является значением граничного условия в начале интервала.
yb
является значением граничного условия в конце интервала.
Для этой проблемы одно из граничных условий для , и другой для . Чтобы вычислить остаточные значения, необходимо поместить граничные условия в форму .
В этой форме граничные условия
,
.
Соответствующая функция затем
function res = emdenbc(ya,yb) res = [ya(2) yb(1) - sqrt(3)/2]; end
Наконец, создайте исходное предположение решения. Для этой проблемы используйте постоянное исходное предположение, которое удовлетворяет граничные условия и простую сетку пяти точек между 0 и 1. Используя многие mesh точки являются ненужными, поскольку решатель BVP адаптирует эти точки во время процесса решения.
,
.
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);
Вычислите значение аналитического решения в .
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 %-------------------------------------------
bvp4c
| bvp5c
| bvpinit
| bvpset