Setup задач

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

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

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

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

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

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

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

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

Чтобы решить это уравнение в 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$ сформируйте левую сторону уравнений mesh (последнее $\mathit{N}$ уравнения в системе). Матрицы блока:

${\mathit{I}}_{\mathit{N}}$ $\mathit{N}\times \mathit{N}$ единичная матрица. Частные производные в ${\mathit{M}}_2$ оцениваются с помощью конечных разностей, в то время как частная производная в ${\mathit{M}}_4$ использует Матрицу Лапласа. Заметьте, что ${\mathit{M}}_3$ содержит только нули, потому что ни одно из уравнений для перемещения mesh не зависит от $\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}$ элементы для движущихся уравнений mesh. Функциональный 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}$ матрица, содержащая все частные производные производной функции, movingMeshODEode15s оценивает якобиевские конечные разности использования, когда матрица не предоставляется в структуре опций. Можно предоставить шаблон разреженности якобиана, чтобы помочь 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] Хуан, Вэйчжан, и др. “Перемещая Методы Mesh На основе Движущихся Дифференциальных уравнений с частными производными Mesh”. Журнал Вычислительной Физики, издания 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
% -------------------------------------------------------------------------