Проверьте согласованность BVP с помощью продолжения

Этот пример показывает, как использовать продолжение для постепенного расширения решения BVP на большие интервалы.

Задачи контура значения Фалькнера-Скана [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

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

Наконец, вы должны предоставить начальное предположение для решения. Неочищенный mesh из пяти точек и постоянное предположение, удовлетворяющее граничным условиям, достаточно хороши, чтобы получить сходимость по интервалу $[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. The axes with title Falkner-Skan Equation, Positive Wall Shear, \beta = 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 для каждой итерации. The 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. The axes with title Falkner-Skan Equation, Positive Wall Shear, \beta = 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. and H. B. Keller. Методы съёмки и параллельной съёмки для решения уравнения контурного слоя Фалькнера-Скана. J. Comp. Phys., Vol. 7, 1971, pp. 289-300.

См. также

bvp4c | bvp5c | bvpinit

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