Setup задачи

Уравнение Бургера является конвекционно-диффузионным уравнением, заданным УЧП

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

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

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

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

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

Условия для движущегося mesh соответствуют MMPDE6 в [1]. Параметр $\tau$ представляет собой шкалу времени для принуждения mesh к равнораспределению. Термин $\mathit{B}\left(\mathit{x},\mathit{t}\right)$ - функция монитора, заданная Eq. 21 в [1]:

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

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

Кодовая Большая матрица

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

$\left[\begin{array}{c}\frac{\partial {\mathit{u}}_1}{\partial \mathit{t}}-\frac{\partial {\mathit{u}}_1}{\partial {\mathit{x}}_1}\frac{\partial {\mathit{x}}_1}{\partial \mathit{t}}\\ ⋮\\ \frac{\partial {\mathit{u}}_{\mathit{N}}}{\partial \mathit{t}}-\frac{\partial {\mathit{u}}_{\mathit{N}}}{\partial {\mathit{x}}_{\mathit{N}}}\frac{\partial {\mathit{x}}_{\mathit{N}}}{\partial \mathit{t}}\\ \frac{{\partial }^2\dot{{\mathit{x}}_1}}{\partial {\mathit{t}}^2}\\ ⋮\\ \frac{{\partial }^2\dot{{\mathit{x}}_{\mathit{N}}}}{\partial {\mathit{t}}^2}\end{array}\right]=\mathit{M}{\mathit{y}}^{\prime }=\left[\begin{array}{cc}{\mathit{M}}_1& {\mathit{M}}_2\\ {\mathit{M}}_3& {\mathit{M}}_4\end{array}\right]\left[\begin{array}{c}\dot{{\mathit{u}}_1}\\ ⋮\\ \dot{{\mathit{u}}_{\mathit{N}}}\\ \dot{{\mathit{x}}_1}\\ ⋮\\ \dot{{\mathit{x}}_{\mathit{N}}}\end{array}\right]$.

Эта формулировка показывает, что ${\mathit{M}}_1$ и ${\mathit{M}}_2$ сформировать левую сторону уравнений Бургерса (первая $\mathit{N}$ уравнения в системе), в то время как ${\mathit{M}}_3$ и ${\mathit{M}}_4$ сформировать левую сторону сетки уравнений (последняя $\mathit{N}$ уравнения в системе). Матрицы блоков:

${\mathit{I}}_{\mathit{N}}$ является $\mathit{N}\times \mathit{N}$ единичная матрица. Частные производные в ${\mathit{M}}_2$ оцениваются с использованием конечных различий, в то время как частная производная в ${\mathit{M}}_4$ использует матрицу Лапласа. Заметьте, что ${\mathit{M}}_3$ содержит только нули, потому что ни одно из уравнений для движения сетки не зависит от $\dot{\mathit{u}}$.

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

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

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

Производная по коду функция

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

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

$\mathit{f}\left(\mathit{u}\right)=ϵ\frac{{\partial }^2\mathit{u}}{{\partial \mathit{x}}^2}-\frac{\partial }{\partial \mathit{x}}\left(\frac{{\mathit{u}}^2}{2}\right)$.

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

По краям mesh (для которой $\mathit{i}=1$ и $\mathit{i}=\mathit{N}$), вместо этого используются только односторонние различия. Этот пример использует $ϵ=1\times {10}^{-4}$.

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

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

После оценки функции монитора применяется пространственное сглаживание (уравнения 14 и 15 в [1]). Этот пример использует $\gamma =2$ и $\mathit{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

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

Якобиан $\mathrm{dF}/\mathrm{dy}$ для производной функции является a $2\mathit{N}\times 2\mathit{N}$ матрица, содержащая все частные производные производной функции, 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{dF}/\mathrm{dy}$ для $\mathit{N}=80$ использование spy.

spy(JPat(80))

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

Функция для шаблона разреженности $\mathit{d}\left(\mathrm{Mv}\right)/\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}\left(\mathrm{Mv}\right)/\mathrm{dy}$ для $\mathit{N}=80$ использование spy.

spy(MvPat(80))

Решающая система уравнений

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

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}$, точки mesh $\mathit{x}\left(\mathit{t}\right)$, и решение $\mathit{u}\left(\mathit{x},\mathit{t}\right)$. Извлеките эти значения из структуры.

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

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

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

Теперь, выборка $\mathit{u}\left(\mathit{x},\mathit{t}\right)$ при нескольких значениях $\mathit{t}$ и постройте график эволюции решения с течением времени. Mesh в концах интервала фиксированы, так что 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}\left(\mathit{x},0\right)$ - это гладкая волна, которая развивает крутой градиент с течением времени, когда она движется к $\mathit{x}=1$. Mesh отслеживают перемещение разрыва, так что дополнительные точки оценки находятся в соответствующем положении на каждом временном шаге.

Ссылки

[1] Huang, Weizhang, et al. Методы скользящей сетки, основанные на движущихся дифференциальных уравнениях с частными производными сетки. Журнал вычислительной физики, том 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
% -------------------------------------------------------------------------