Настройка проблемы

Уравнение Бургера - уравнение конвекции-диффузии, данное PDE

Применение преобразования координат (Eq. 18 в [1]) приводит к дополнительному члену в левой части:

Преобразование PDE в ODE одной переменной осуществляется с использованием конечных разностей для аппроксимации частных производных, взятых по отношению к $\mathit{}$x. Если конечные разности записаны как $\Delta$, то PDE может быть переписан как ODE, который содержит только производные, взятые по отношению к $\mathit{t}$:

В этой форме можно использовать решатель ОДУ, например ode15s для решения для $\mathit{u}$ и $\mathit{x}$ с течением времени.

Для этого примера проблема формулируется на движущейся сетке из $\mathit{N}$ точек, и метод движущейся сетки, описанный в [1], позиционирует точки сетки на каждом временном этапе так, что они концентрируются в областях изменения. Граница и исходные условия:

Для данной начальной сетки из $\mathit{N}$ точек существуют $\mathrm{}\mathit{2N}$ уравнения для решения: $\mathit{N}$ уравнений, соответствующих уравнению Бургера, и $\mathit{N}$ уравнений, определяющих движение каждой точки сетки. Итак, конечная система уравнений:

Термины для движущейся сетки соответствуют MMPDE6 в [1]. Параметр $\mathrm{startпредставляет}$ собой временную шкалу для принудительной передачи сетки в сторону эквидистантного распределения. Термин $\mathit{B}\mathit{(}x\mathit{,}$t) представляет собой функцию контроля, заданную уравнением 21 в [1]:

Подход, используемый в этом примере для решения уравнения Бургера с движущимися точками сетки, демонстрирует несколько методов:

Чтобы решить это уравнение в MATLAB ®, запишите производную функцию, массовую матричную функцию, функцию для шаблона разреженности Jacobian $\mathrm{}\mathrm{}$dF/dy и функцию для шаблона разреженности $\mathit{}d\mathrm{(}Mv\mathrm{)}$/dy. Требуемые функции можно либо включить в качестве локальных функций в конце файла (как здесь сделано), либо сохранить их как отдельные именованные файлы в каталоге по пути MATLAB.

Матрица массы кода

Левая часть системы уравнений включает линейные комбинации первых производных, поэтому для представления всех членов требуется массовая матрица. Установите левую сторону системы уравнений равной $\mathit{M}{\mathit{y}}^{}$′ для извлечения формы матрицы масс. Массовая матрица состоит из четырех блоков, каждый из которых является квадратной матрицей порядка $\mathit{}$N :

$[\begin{array}{c}\frac{{\mathit{}}_{\mathrm{}}}{\mathit{}}\frac{{\mathit{}}_{\mathrm{}}}{{\mathit{}}_{\mathrm{}}}\frac{{\mathit{}}_{\mathrm{}}}{\mathit{}}\\ \\ \frac{{\mathit{}}_{\mathit{}}}{\mathit{}}\frac{{\mathit{}}_{\mathit{}}}{{\mathit{}}_{\mathit{}}}\frac{{\mathit{}}_{\mathit{}}}{\mathit{}}\\ \frac{{}^{\mathrm{}}\stackrel{}{{\mathit{}}_{\mathrm{}}}}{{\mathit{}}^{\mathrm{}}}\\ \\ \frac{{}^{\mathrm{}}\stackrel{}{{\mathit{}}_{\mathit{}}}}{{\mathit{}}^{\mathrm{}}}\end{array}\mathit{\partial u1\partial t-\partial u1\partial x1\partial x1\partial t⋮\partial uN\partial t-\partial uN\partial xN\partial xN\partial t\partial 2x1\dot{}\partial t2⋮\partial 2xN\dot{}\partial t2]=M}{\mathit{}}^{}\begin{array}{cc}{\mathit{}}_{\mathrm{}}& {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}& {\mathit{}}_{\mathrm{y\prime =[M1M2M3M4}}\end{array}]\begin{array}{c}\stackrel{}{{\mathit{[}}_{\mathrm{}}}\\ \\ \stackrel{}{{\mathit{}}_{\mathit{}}}\\ \stackrel{}{{\mathit{}}_{\mathrm{}}}\\ \\ \stackrel{}{{\mathit{}}_{\mathit{}}}\end{array}u1\dot{}⋮uN\dot{}x1\dot{}⋮xN\dot{}$].

Эта формулировка показывает, что ${\mathit{}}_{\mathrm{M1}}$ и ${\mathit{}}_{\mathrm{M2}}$ образуют левую сторону уравнений Бургера (первые $\mathit{N}$ уравнений в системе), в то время как ${\mathit{}}_{\mathrm{M3}}$ и ${\mathit{}}_{\mathrm{M4}}$ образуют левую сторону уравнений сетки (последние $\mathit{N}$ уравнений в системе). Блочные матрицы:

${\mathit{}}_{\mathit{IN}}$ - единичная $\mathit{}\mathit{}$матрица N × N. Частные производные ${\mathit{в}}_{\mathrm{}}$M2 оцениваются с использованием конечных разностей, в то время как частная производная ${\mathit{в}}_{\mathrm{}}$M4 использует матрицу Лапласа. Обратите внимание, ${\mathit{что}}_{\mathrm{}}$M3 содержит только нули, поскольку ни одно из уравнений для перемещения сетки не зависит от $\stackrel{}{\mathit{}}$u˙.

Теперь можно записать функцию, которая вычисляет массовую матрицу. Функция должна принимать два входа для времени $\mathit{t}$ и вектора решения $\mathit{y}$. Поскольку вектор решения $\mathit{y}$ содержит половину $\stackrel{u\dot{}}{\mathit{}}$ компонентов и половину $\stackrel{x\dot{}}{\mathit{}}$ компонентов, функция извлекает их первой. Затем функция формирует все блочные матрицы (принимая во внимание граничные значения задачи) и собирает массовую матрицу, используя четыре блока.

function M = mass(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2; 
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% M1 and M2 are the portions of the mass matrix for Burgers' equation. 
% The derivative du/dx is approximated with finite differences, using 
% single-sided differences on the edges and centered differences in between.
M1 = speye(N);
M2 = sparse(N,N);
M2(1,1) = - (u(2) - u0)/(x(2) - x0);
for i = 2:N-1
    M2(i,i) = - (u(i+1) - u(i-1))/(x(i+1) - x(i-1));
end
M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));

% M3 and M4 define the equations for mesh point evolution, corresponding to 
% MMPDE6 in the reference paper. Since the mesh functions only involve d/dt(dx/dt), 
% the M3 portion of the mass matrix is all zeros. The second derivative in M4 is 
% approximated using a finite difference Laplacian matrix. 
M3 = sparse(N,N);
e = ones(N,1);
M4 = spdiags([e -2*e e],-1:1,N,N);

% Assemble mass matrix
M = [M1 M2
     M3 M4];
end

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

Функция кодовой производной

Производная функция для этой задачи возвращает вектор с $\mathrm{}\mathit{2N}$ элементами. Первые $\mathit{N}$ элементов соответствуют уравнениям Бургера, в то время как последние $\mathit{N}$ элементов являются уравнениями движущейся сетки. Функция movingMeshODE для оценки правых сторон всех уравнений в системе:

Первые $\mathit{N}$ уравнений в производной функции кодируют правую сторону уравнений Бургера. Уравнения бургеров можно рассматривать как дифференциальный оператор, включающий пространственные производные вида:

$\mathit{f}\mathit{(}u)\frac{{}^{\mathrm{}}\mathit{}}{{\mathit{}}^{\mathrm{}}}\frac{}{\mathit{}}\frac{{\mathit{=ϵ\partial 2u\partial x2-\partial \partial x}}^{\mathrm{}}}{\mathrm{(}}$u22 ).

В справочном документе [1] описан процесс аппроксимации дифференциального оператора $\mathit{f}$ с использованием центрированных конечных разностей на

На кромках сетки (для которой $\mathit{i}\mathrm{=}$ 1 $\mathit{и}\mathit{i}$= N) используются только односторонние различия. В этом примере $\mathrm{используется}{\mathrm{}}^{\mathrm{}}$ϵ=1×10-4.

Уравнения, определяющие сетку (состоящую из последних $\mathit{N}$ уравнений в производной функции):

Подобно уравнениям Бургера, можно использовать конечные разности для аппроксимации функции монитора $\mathit{B}\mathit{(}x\mathit{,}$t):

После оценки функции контроля применяется пространственное сглаживание (уравнения 14 и 15 в [1]). В этом примере для $\mathrm{}$$\mathit{}\mathrm{}$параметров пространственного сглаживания используются γ = 2 и p = 2.

Функция, кодирующая систему уравнений,

function g = movingMeshODE(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2;
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% Preallocate g vector of derivative values.
g = zeros(2*N,1);

% Use centered finite differences to approximate the RHS of Burgers'
% equations (with single-sided differences on the edges). The first N
% elements in g correspond to Burgers' equations.
for i = 2:N-1
    delx = x(i+1) - x(i-1);
    g(i) = 1e-4*((u(i+1) - u(i))/(x(i+1) - x(i)) - ...
        (u(i) - u(i-1))/(x(i) - x(i-1)))/(0.5*delx) ...
        - 0.5*(u(i+1)^2 - u(i-1)^2)/delx;
end
delx = x(2) - x0;
g(1) = 1e-4*((u(2) - u(1))/(x(2) - x(1)) - (u(1) - u0)/(x(1) - x0))/(0.5*delx) ...
    - 0.5*(u(2)^2 - u0^2)/delx;
delx = xNP1 - x(N-1);
g(N) = 1e-4*((uNP1 - u(N))/(xNP1 - x(N)) - ...
    (u(N) - u(N-1))/(x(N) - x(N-1)))/delx - ...
    0.5*(uNP1^2 - u(N-1)^2)/delx;

% Evaluate the monitor function values (Eq. 21 in reference paper), used in
% RHS of mesh equations. Centered finite differences are used for interior
% points, and single-sided differences are used on the edges.
M = zeros(N,1);
for i = 2:N-1
    M(i) = sqrt(1 + ((u(i+1) - u(i-1))/(x(i+1) - x(i-1)))^2);
end
M0 = sqrt(1 + ((u(1) - u0)/(x(1) - x0))^2);
M(1) = sqrt(1 + ((u(2) - u0)/(x(2) - x0))^2);
M(N) = sqrt(1 + ((uNP1 - u(N-1))/(xNP1 - x(N-1)))^2);
MNP1 = sqrt(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);

% Apply spatial smoothing (Eqns. 14 and 15) with gamma = 2, p = 2.
SM = zeros(N,1);
for i = 3:N-2
    SM(i) = sqrt((4*M(i-2)^2 + 6*M(i-1)^2 + 9*M(i)^2 + ...
        6*M(i+1)^2 + 4*M(i+2)^2)/29);
end
SM0 = sqrt((9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);
SM(1) = sqrt((6*M0^2 + 9*M(1)^2 + 6*M(2)^2 + 4*M(3)^2)/25);
SM(2) = sqrt((4*M0^2 + 6*M(1)^2 + 9*M(2)^2 + 6*M(3)^2 + 4*M(4)^2)/29);
SM(N-1) = sqrt((4*M(N-3)^2 + 6*M(N-2)^2 + 9*M(N-1)^2 + 6*M(N)^2 + 4*MNP1^2)/29);
SM(N) = sqrt((4*M(N-2)^2 + 6*M(N-1)^2 + 9*M(N)^2 + 6*MNP1^2)/25);
SMNP1 = sqrt((4*M(N-1)^2 + 6*M(N)^2 + 9*MNP1^2)/19);
for i = 2:N-1
    g(i+N) = (SM(i+1) + SM(i))*(x(i+1) - x(i)) - ...
        (SM(i) + SM(i-1))*(x(i) - x(i-1));
end
g(1+N) = (SM(2) + SM(1))*(x(2) - x(1)) - (SM(1) + SM0)*(x(1) - x0);
g(N+N) = (SMNP1 + SM(N))*(xNP1 - x(N)) - (SM(N) + SM(N-1))*(x(N) - x(N-1));

% Form final discrete approximation for Eq. 12 in reference paper, the equation governing
% the mesh points.
tau = 1e-3;
g(1+N:end) = - g(1+N:end)/(2*tau);
end

Функции кода для шаблонов разреженности

Jacobian $\mathrm{}\mathrm{dF/dy}$ для производной функции является $\mathrm{}\mathit{}\mathrm{}\mathit{}$матрицей 2N × 2N, содержащей все частные производные производной функции,movingMeshODE. ode15s оценивает якобиан, используя конечные различия, когда матрица не представлена в структуре опций. Вы можете предоставить разреженность модели якобиана, чтобы помочь ode15s рассчитать его быстрее.

Функция для разреженности узора якобиана

function out = JPat(N)
S1 = spdiags(ones(N,3),-1:1,N,N);
S2 = spdiags(ones(N,9),-4:4,N,N);
out = [S1 S1
       S2 S2];
end

Постройте график разреженности $\mathrm{}\mathrm{dF/dy}$ для $\mathit{N}\mathrm{=}$ 80 с помощьюspy.

spy(JPat(80))

Другой способ сделать расчет более эффективным заключается в обеспечении схемы разреженности $\mathit{d}\mathrm{(}Mv\mathrm{)/}$dy. Этот узор разреженности можно найти, изучив, какие члены ${\mathit{}}_{\mathit{}}$ui ${\mathit{и}}_{\mathit{}}$xi присутствуют в конечных разностях, вычисленных в функции матрицы масс.

Функция для шаблона разреженности $\mathit{d}\mathrm{(}Mv\mathrm{)/}$dy

function S = MvPat(N)
S = sparse(2*N,2*N);
S(1,2) = 1;
S(1,2+N) = 1;
for i = 2:N-1
    S(i,i-1) = 1;
    S(i,i+1) = 1;
    S(i,i-1+N) = 1;
    S(i,i+1+N) = 1;
end
S(N,N-1) = 1;
S(N,N-1+N) = 1;
end

Постройте график разреженности $\mathit{d}\mathrm{(}Mv\mathrm{)/}$dy $\mathit{для}\mathrm{N}$ = 80 с использованиемspy.

spy(MvPat(80))

Система решений уравнений

Вычислите систему со значением $\mathit{N}\mathrm{=}$80. Для начальных условий $\mathit{инициализируйте}$ x с помощью однородной сетки и $\mathit{вычислите}\mathit{u}\mathrm{(}x$, 0) на сетке.

N = 80;
h = 1/(N+1);
xinit = h*(1:N);
uinit = sin(2*pi*xinit) + 0.5*sin(pi*xinit);
y0 = [uinit xinit];

Использовать odeset для создания структуры опций, которая задает несколько значений:

opts = odeset('Mass',@mass,'MStateDependence','strong','JPattern',JPat(N),...
   'MvPattern',MvPat(N),'RelTol',1e-5,'AbsTol',1e-4);

Наконец, позвоните ode15s для решения системы на интервале $\left[0,1\right]$ с помощью movingMeshODE производная функция, временной интервал, начальные условия и структура опций.

tspan = [0 1];
sol = ode15s(@movingMeshODE,tspan,y0,opts);

Результаты графика

Результатом интеграции является структура sol содержит временные шаги $\mathit{t}$, точки сетки $\mathit{x}\mathit{(}t$) и решение $\mathit{}u\mathit{}(\mathit{x},$t). Извлеките эти значения из структуры.

t = sol.x;
x = sol.y(N+1:end,:);
u = sol.y(1:N,:);

Постройте график перемещения точек сетки по времени. График показывает, что точки сетки сохраняют достаточно ровный интервал во времени (благодаря функции монитора), но они способны кластеризоваться вблизи разрыва в решении по мере его перемещения.

plot(x,t)
xlabel('t')
ylabel('x(t)')
title('Burgers'' equation: Trajectories of grid points')

Теперь, образец $\mathit{u}\mathit{(}x\mathit{,}$t) при нескольких значениях $\mathit{}$t и график эволюции раствора во времени. Точки сетки на концах интервала фиксированы, поэтомуx(0) = 0 и x(N+1) = 1. Граничные значения: u(t,0) = 0 и u(t,1) = 0, которые необходимо добавить к известным значениям, вычисленным для фигуры.

tint = 0:0.2:1;
yint = deval(sol,tint);
figure
labels = {};
for j = 1:length(tint)
   solution = [0; yint(1:N,j); 0];
   location = [0; yint(N+1:end,j); 1];
   labels{j} = ['t = ' num2str(tint(j))];
   plot(location,solution,'-o')
   hold on
end
xlabel('x')
ylabel('solution u(x,t)')
legend(labels{:},'Location','SouthWest')
title('Burgers equation on moving mesh')
hold off

График показывает, что $\mathit{u}\mathit{(}x\mathrm{,}$0) является гладкой волной, которая развивает крутой градиент с течением времени по мере движения $\mathit{к}\mathrm{x}$= 1. Точки сетки отслеживают движение разрыва, так что дополнительные точки оценки находятся в соответствующем положении на каждом временном шаге.

Ссылки

[1] Huang, Weizhang и др. «Методы перемещения сетки на основе дифференциальных уравнений в частных производных сетки». Журнал вычислительной физики, том 113, № 2, август 1994, стр. 279-90. https://doi.org/10.1006/jcph.1994.1135.

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

Здесь перечислены локальные вспомогательные функции, которые решатель ode15s вызывает для вычисления решения. Кроме того, эти функции можно сохранить в виде собственных файлов в каталоге по пути MATLAB.

function g = movingMeshODE(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2;
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% Preallocate g vector of derivative values.
g = zeros(2*N,1);

% Use centered finite differences to approximate the RHS of Burgers'
% equations (with single-sided differences on the edges). The first N
% elements in g correspond to Burgers' equations.
for i = 2:N-1
    delx = x(i+1) - x(i-1);
    g(i) = 1e-4*((u(i+1) - u(i))/(x(i+1) - x(i)) - ...
        (u(i) - u(i-1))/(x(i) - x(i-1)))/(0.5*delx) ...
        - 0.5*(u(i+1)^2 - u(i-1)^2)/delx;
end
delx = x(2) - x0;
g(1) = 1e-4*((u(2) - u(1))/(x(2) - x(1)) - (u(1) - u0)/(x(1) - x0))/(0.5*delx) ...
    - 0.5*(u(2)^2 - u0^2)/delx;
delx = xNP1 - x(N-1);
g(N) = 1e-4*((uNP1 - u(N))/(xNP1 - x(N)) - ...
    (u(N) - u(N-1))/(x(N) - x(N-1)))/delx - ...
    0.5*(uNP1^2 - u(N-1)^2)/delx;

% Evaluate the monitor function values (Eq. 21 in reference paper), used in
% RHS of mesh equations. Centered finite differences are used for interior
% points, and single-sided differences are used on the edges.
M = zeros(N,1);
for i = 2:N-1
    M(i) = sqrt(1 + ((u(i+1) - u(i-1))/(x(i+1) - x(i-1)))^2);
end
M0 = sqrt(1 + ((u(1) - u0)/(x(1) - x0))^2);
M(1) = sqrt(1 + ((u(2) - u0)/(x(2) - x0))^2);
M(N) = sqrt(1 + ((uNP1 - u(N-1))/(xNP1 - x(N-1)))^2);
MNP1 = sqrt(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);

% Apply spatial smoothing (Eqns. 14 and 15) with gamma = 2, p = 2.
SM = zeros(N,1);
for i = 3:N-2
    SM(i) = sqrt((4*M(i-2)^2 + 6*M(i-1)^2 + 9*M(i)^2 + ...
        6*M(i+1)^2 + 4*M(i+2)^2)/29);
end
SM0 = sqrt((9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);
SM(1) = sqrt((6*M0^2 + 9*M(1)^2 + 6*M(2)^2 + 4*M(3)^2)/25);
SM(2) = sqrt((4*M0^2 + 6*M(1)^2 + 9*M(2)^2 + 6*M(3)^2 + 4*M(4)^2)/29);
SM(N-1) = sqrt((4*M(N-3)^2 + 6*M(N-2)^2 + 9*M(N-1)^2 + 6*M(N)^2 + 4*MNP1^2)/29);
SM(N) = sqrt((4*M(N-2)^2 + 6*M(N-1)^2 + 9*M(N)^2 + 6*MNP1^2)/25);
SMNP1 = sqrt((4*M(N-1)^2 + 6*M(N)^2 + 9*MNP1^2)/19);
for i = 2:N-1
    g(i+N) = (SM(i+1) + SM(i))*(x(i+1) - x(i)) - ...
        (SM(i) + SM(i-1))*(x(i) - x(i-1));
end
g(1+N) = (SM(2) + SM(1))*(x(2) - x(1)) - (SM(1) + SM0)*(x(1) - x0);
g(N+N) = (SMNP1 + SM(N))*(xNP1 - x(N)) - (SM(N) + SM(N-1))*(x(N) - x(N-1));

% Form final discrete approximation for Eq. 12 in reference paper, the equation governing
% the mesh points.
tau = 1e-3;
g(1+N:end) = - g(1+N:end)/(2*tau);
end

% -----------------------------------------------------------------------
function M = mass(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2; 
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% M1 and M2 are the portions of the mass matrix for Burgers' equation. 
% The derivative du/dx is approximated with finite differences, using 
% single-sided differences on the edges and centered differences in between.
M1 = speye(N);
M2 = sparse(N,N);
M2(1,1) = - (u(2) - u0)/(x(2) - x0);
for i = 2:N-1
    M2(i,i) = - (u(i+1) - u(i-1))/(x(i+1) - x(i-1));
end
M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));

% M3 and M4 define the equations for mesh point evolution, corresponding to 
% MMPDE6 in the reference paper. Since the mesh functions only involve d/dt(dx/dt), 
% the M3 portion of the mass matrix is all zeros. The second derivative in M4 is 
% approximated using a finite difference Laplacian matrix. 
M3 = sparse(N,N);
e = ones(N,1);
M4 = spdiags([e -2*e e],-1:1,N,N);

% Assemble mass matrix
M = [M1 M2
     M3 M4];
end

% -------------------------------------------------------------------------
function out = JPat(N)  % Jacobian sparsity pattern
S1 = spdiags(ones(N,3),-1:1,N,N);
S2 = spdiags(ones(N,9),-4:4,N,N);
out = [S1 S1
    S2 S2];
end

% -------------------------------------------------------------------------
function S = MvPat(N)  % Sparsity pattern for the derivative of the Mass matrix times a vector
S = sparse(2*N,2*N);
S(1,2) = 1;
S(1,2+N) = 1;
for i = 2:N-1
    S(i,i-1) = 1;
    S(i,i+1) = 1;
    S(i,i-1+N) = 1;
    S(i,i+1+N) = 1;
end
S(N,N-1) = 1;
S(N,N-1+N) = 1;
end
% -------------------------------------------------------------------------