Проверьте непротиворечивость BVP Используя продолжение

В этом примере показано, как использовать продолжение, чтобы постепенно расширить решение BVP больших интервалов.

Краевые задачи Falkner-Skan [1] являются результатом решений для подобия вязкого, несжимаемого, ламинарного течения по плоской пластине. Уравнение в качестве примера

${f^{'}}^{'}^{'} + f {f^{'}}^{'} + β (1 - {f^{'}}^{2}) = 0$ .

Проблема создана на бесконечном интервале $[0, \infty]$ с $β = 0.5$ Согласно граничным условиям

$f (0) = 0$ ,

$f^{'} (0) = 0$ ,

$f^{'} (\infty) = 1$ .

BVP не может быть решен на бесконечном интервале, и это невозможно решить BVP в очень большом конечном интервале. Вместо этого этот пример решает последовательность проблем, созданных на меньшем интервале $[0, a]$ проверять, что решение имеет сопоставимое поведение как $a \to \infty$ . Эта практика разбивания проблемы в более простые проблемы, с решением каждой проблемы, возвращающейся в как исходное предположение для следующей проблемы, называется продолжением.

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

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

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

x является независимой переменной.
f является зависимой переменной.

Можно переписать уравнение третьего порядка как систему уравнений первого порядка с помощью замен $f_{1} = f$ , $f_{2} = f^{'}$ , и $f_{3} = {f^{'}}^{'}$ . Уравнения становятся

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

${f_{2}}^{'} = f_{3}$ ,

${f_{3}}^{'} = - f_{1} f_{3} - β (1 - f_{2}^{2})$ .

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

function dfdeta = fsode(x,f)
b = 0.5;
dfdeta = [ f(2)
           f(3)
          -f(1)*f(3) - b*(1 - f(2)^2) ];
end

Примечание: Все функции включены как локальные функции в конце примера.

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

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

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

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

$f (0) = 0$ ,

$f^{'} (0) = 0$ ,

$f^{'} (\infty) - 1 = 0$ .

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

function res = fsbc(f0,finf)
res = [f0(1)
       f0(2)
       finf(2) - 1];
end

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

Наконец, необходимо обеспечить исходное предположение для решения. Грубая сетка пяти точек и постоянного предположения, которое удовлетворяет граничным условиям, достаточно хороша, чтобы получить сходимость на интервале $[0, 3]$ . Переменная infinity обозначает правый предел интервала интегрирования. Как значение infinity увеличения на последующих итерациях от 3 до его максимального значения 6, решение от каждой предыдущей итерации действует как исходное предположение для следующей итерации.

infinity = 3;
maxinfinity = 6;
solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);

Решите уравнение и постройте решения

Решите задачу в начальном интервале $[0, 3]$ . Постройте решение и сравните значение ${f^{'}}^{'} (0)$ к аналитическому значению [1].

sol = bvp4c(@fsode,@fsbc,solinit);
x = sol.x;
f = sol.y;
plot(x,f(2,:),x(end),f(2,end),'o');
axis([0 maxinfinity 0 1.4]);
title('Falkner-Skan Equation, Positive Wall Shear, \beta = 0.5.')
xlabel('x')
ylabel('df/dx')
hold on

Figure contains an axes object. The axes object with title F a l k n e r - S k a n blank E q u a t i o n , blank P o s i t i v e blank W a l l blank S h e a r , blank beta blank = blank 0 . 5 . contains 2 objects of type line.

fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n')

Cebeci & Keller report that f''(0) = 0.92768.

fprintf('Value computed using infinity = %g is %7.5f.\n', ...
         infinity,f(3,1))

Value computed using infinity = 3 is 0.92915.

Теперь решите задачу на прогрессивно больших интервалах путем увеличения значения infinity для каждой итерации. bvpinit функция экстраполирует каждое решение нового интервала, чтобы действовать как исходное предположение для следующего значения infinity. Каждая итерация распечатывает расчетное значение ${f^{'}}^{'} (0)$ и накладывает график решения по предыдущим решениям. Когда infinity = 6, сопоставимое поведение решения становится очевидным и значение ${f^{'}}^{'} (0)$ очень близко к ожидаемому значению.

for Bnew = infinity+1:maxinfinity
  solinit = bvpinit(sol,[0 Bnew]); % Extend solution to new interval
  sol = bvp4c(@fsode,@fsbc,solinit);
  x = sol.x;
  f = sol.y;
  
  fprintf('Value computed using infinity = %g is %7.5f.\n', ...
           Bnew,f(3,1))
  plot(x,f(2,:),x(end),f(2,end),'o');
  drawnow
end

Value computed using infinity = 4 is 0.92774.
Value computed using infinity = 5 is 0.92770.
Value computed using infinity = 6 is 0.92770.

hold off

Figure contains an axes object. The axes object with title F a l k n e r - S k a n blank E q u a t i o n , blank P o s i t i v e blank W a l l blank S h e a r , blank beta blank = blank 0 . 5 . contains 8 objects of type line.

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

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

function dfdeta = fsode(x,f) % equation being solved
dfdeta = [ f(2)
           f(3)
          -f(1)*f(3) - 0.5*(1 - f(2)^2) ];
end 
%-------------------------------------------
function res = fsbc(f0,finf) % boundary conditions
res = [f0(1)
       f0(2)
       finf(2) - 1];
end
%-------------------------------------------

Ссылки

[1] Cebeci, T. и Х. Б. Келлер. "Стрельба и Параллельные Методы Стрельбы для Решения Уравнения пограничного слоя Falkner-Skan". J. Аккомпанируйте Физике, Изданию 7, 1971, стр 289-300.

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

bvp4c | bvp5c | bvpinit

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