Кодовое уравнение

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

Можно переписать уравнение третьего порядка как систему уравнений первого порядка, используя замены ${\mathit{}}_{\mathrm{f1}}\mathit{=}$f${\mathit{,}}_{\mathrm{}}{\mathit{f2}}^{}$= f ′ ${\mathit{}}_{\mathrm{}}{{\mathit{и}}^{}}^{}$f3 = f ′ ′. Уравнения становятся

${{\mathit{}}_{\mathrm{}}}^{f1\; }={\mathit{}}_{\mathrm{f2}}$,

${{\mathit{}}_{\mathrm{}}}^{f2\; }={\mathit{}}_{\mathrm{f3}}$,

${{\mathit{}}_{\mathrm{}}}^{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\mathrm{f3\prime =-f1f3-\beta }\mathrm{(}{\mathit{}}_{\mathrm{}}^{\mathrm{}}1-f22$).

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

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), где:

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

$\mathit{f}\mathrm{(}0)\mathrm{}$= 0,

${\mathit{f}}^{}\prime \mathrm{}(0\mathrm{)}$= 0,

${\mathit{f}}^{}\prime (\infty \mathrm{)}\mathrm{-}$1 = 0.

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

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

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

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

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

Решения для уравнений и графиков

Решите проблему в начальном интервале $[0,3$]. Постройте график решения и сравните значение ${{\mathit{f}}^{}}^{\prime }\mathrm{\prime }$(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

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. Каждая итерация печатает вычисленное значение ${{\mathit{f}}^{}}^{\prime }\mathrm{\prime }$(0) и накладывает график решения на предыдущие решения. Когдаinfinity = 6непротиворечивое поведение решения становится очевидным, и значение ${{\mathit{f}}^{}}^{\prime }\mathrm{\prime }$(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

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

Здесь перечислены локальные вспомогательные функции, которые решатель 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. и Х. Б. Келлер. «Методы съемки и параллельной съемки для решения уравнения граничного слоя Фалькнера-Скана». J. Comp. Phys., том 7, 1971, стр. 289-300.