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

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

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

y+2xy+y5=0.

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

y(x)=[1+x23]-1.

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

y=Syx+f(x,y).

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

Можно переписать уравнение Эмдена как систему использования уравнений первого порядка y1=y и y2=y как

y1=y2,

y2=-2xy2-y15.

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

y2(0)=0,

y1(1)=3/2.

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

[y1y2]=1x[000-2][y1y2]+[y2-y15].

В матричной форме это ясно это S=[000-2] и f(x,y)=[y2-y15].

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

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

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

  • x независимая переменная.

  • y зависимая переменная.

  • dydx(1) дает уравнение для y1, и dydx(2) уравнение для y2.

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

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

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

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

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

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

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

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

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

y2(0)=0,

y1(1)-3/2=0.

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

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

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

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

y1=3/2,

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);

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

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

Смотрите также

| | |

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте