Solve ODE with Strongly State-Dependent Mass Matrix

This example shows how to solve Burgers' equation using a moving mesh technique [1]. The problem includes a mass matrix, and options are specified to account for the strong state dependence and sparsity of the mass matrix, making the solution process more efficient.

Problem Setup

Burgers' equation is a convection-diffusion equation given by the PDE

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}=ϵ\frac{{\partial }^2\mathit{u}}{{\partial \mathit{x}}^2}-\frac{\partial }{\partial \mathit{x}}\left(\frac{{\mathit{u}}^2}{2}\right),0<\mathit{x}<1,\mathit{t}>0,ϵ=1{\times \mathrm{10}}^{-4}.$$

Applying a coordinate transformation (Eq. 18 in [1]) leads to an extra term on the left-hand side:

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

Converting the PDE into an ODE of one variable is accomplished by using finite differences to approximate the partial derivatives taken with respect to $\mathit{x}$. If the finite differences are written as $\Delta$, then the PDE can be rewritten as an ODE that only contains derivatives taken with respect to $\mathit{t}$:

$$\frac{\mathit{du}}{\mathit{dt}}-\Delta \mathit{u}\frac{\mathit{dx}}{\mathit{dt}}=ϵ{\Delta }^2\mathit{u}-\Delta \left(\frac{{\mathit{u}}^2}{2}\right).$$

In this form, you can use an ODE solver such as ode15s to solve for $\mathit{u}$ and $\mathit{x}$ over time.

For this example, the problem is formulated on a moving mesh of $\mathit{N}$ points, and the moving mesh technique described in [1] positions the mesh points at each time step so that they are concentrated in areas of change. The boundary and initial conditions are

$$\begin{array}{l}\mathit{u}\left(0,\mathit{t}\right)=\mathit{u}\left(1,\mathit{t}\right)=0,\\ \mathit{u}\left(\mathit{x},0\right)=\sin \left(2\pi \mathit{x}\right)+\frac{1}{2}\sin \left(\pi \mathit{x}\right).\end{array}$$

For a given initial mesh of $\mathit{N}$ points, there are $2\mathit{N}$ equations to solve: $\mathit{N}$ equations corresponding to Burgers' equation, and $\mathit{N}$ equations determining the movement of each mesh point. So, the final system of equations is:

$$\begin{array}{l}\frac{\mathit{d}{\mathit{u}}_1}{\mathit{dt}}-\Delta {\mathit{u}}_1\frac{\mathit{d}{\mathit{x}}_1}{\mathit{dt}}=ϵ{\Delta }^2{\mathit{u}}_1-\Delta \left(\frac{{{\mathit{u}}_1}^2}{2}\right),\\ ⋮\\ \frac{\mathit{d}{\mathit{u}}_{\mathit{N}}}{\mathit{dt}}-\Delta {\mathit{u}}_{\mathit{N}}\frac{\mathit{d}{\mathit{x}}_{\mathit{N}}}{\mathit{dt}}=ϵ{\Delta }^2{\mathit{u}}_{\mathit{N}}-\Delta \left(\frac{{{\mathit{u}}_{\mathit{N}}}^2}{2}\right),\\ \frac{{\mathit{d}}^2\dot{{\mathit{x}}_1}}{\mathit{d}{\mathit{t}}^2}=\frac{1}{\tau }\frac{\mathit{d}}{\mathit{dt}}\left(\mathit{B}\left({\mathit{x}}_1,\mathit{t}\right)\frac{\mathit{d}{\mathit{x}}_1}{\mathit{dt}}\right),\\ ⋮\\ \frac{{\mathit{d}}^2\dot{{\mathit{x}}_{\mathit{N}}}}{\mathit{d}{\mathit{t}}^2}=\frac{1}{\tau }\frac{\mathit{d}}{\mathit{dt}}\left(\mathit{B}\left({\mathit{x}}_{\mathit{N}},\mathit{t}\right)\frac{\mathit{d}{\mathit{x}}_{\mathit{N}}}{\mathit{dt}}\right).\end{array}$$

The terms for the moving mesh correspond to MMPDE6 in [1]. The parameter $\tau$ represents a timescale for forcing the mesh toward equidistribution. The term $\mathit{B}\left(\mathit{x},\mathit{t}\right)$ is a monitor function given by Eq. 21 in [1]:

$$\mathit{B}\left(\mathit{x},\mathit{t}\right)=\sqrt{1+{\left(\frac{\mathit{d}{\mathit{u}}_{\mathit{i}}}{\mathit{d}{\mathit{x}}_{\mathit{i}}}\right)}^2}.$$

The approach used in this example to solve Burgers' equation with moving mesh points demonstrates several techniques:

• The system of equations is expressed using a mass matrix formulation, $\mathit{M}{\mathit{y}}^{\prime }=\mathit{f}\left(\mathit{t},\mathit{y}\right)$. The mass matrix is provided to the ode15s solver as a function.

• The derivative function not only includes the equations for Burgers' equation, but also a set of equations governing the moving mesh selection.

• The sparsity patterns of the Jacobian $\mathrm{dF}/\mathrm{dy}$ and the derivative of the mass matrix multiplied with a vector $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$ are supplied to the solver as functions. Supplying these sparsity patterns helps the solver operate more efficiently.

• Finite differences are used to approximate several partial derivatives.

To solve this equation in MATLAB®, write a derivative function, a mass matrix function, a function for the sparsity pattern of the Jacobian $\mathrm{dF}/\mathrm{dy}$, and a function for the sparsity pattern of $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$. You can either include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

Code Mass Matrix

The left side of the system of equations involves linear combinations of first derivatives, so a mass matrix is required to represent all of the terms. Set the left side of the system of equations equal to $\mathit{M}{\mathit{y}}^{\prime }$ to extract the form of the mass matrix. The mass matrix is composed of four blocks, each of which is a square matrix of order $\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]$.

This formulation shows that ${\mathit{M}}_1$ and ${\mathit{M}}_2$ form the left side of Burgers' equations (the first $\mathit{N}$ equations in the system), while ${\mathit{M}}_3$ and ${\mathit{M}}_4$ form the left side of the mesh equations (the last $\mathit{N}$ equations in the system). The block matrices are:

$$\begin{array}{l}{\mathit{M}}_1={\mathit{I}}_{\mathit{N}},\\ {\mathit{M}}_2=-\frac{\partial {\mathit{u}}_{\mathit{i}}}{\partial {\mathit{x}}_{\mathit{i}}}{\mathit{I}}_{\mathit{N}},\\ {\mathit{M}}_3=0_{\mathit{N}},\\ {\mathit{M}}_4=\frac{{\partial }^2}{\partial {\mathit{t}}^2}{\mathit{I}}_{\mathit{N}}.\end{array}$$

${\mathit{I}}_{\mathit{N}}$ is the $\mathit{N}\times \mathit{N}$ identity matrix. The partial derivatives in ${\mathit{M}}_2$ are estimated using finite differences, while the partial derivative in ${\mathit{M}}_4$ uses a Laplacian matrix. Notice that ${\mathit{M}}_3$ contains only zeros because none of the equations for the mesh movement depend on $\dot{\mathit{u}}$.

Now you can write a function that computes the mass matrix. The function must accept two inputs for time $\mathit{t}$ and the solution vector $\mathit{y}$. Since the solution vector $\mathit{y}$ contains half $\dot{\mathit{u}}$ components and half $\dot{\mathit{x}}$ components, the function extracts these first. Then, the function forms all of the block matrices (taking the boundary values of the problem into account) and assembles the mass matrix using the four blocks.

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

Note: All functions are included as local functions at the end of the example.

Code Derivative Function

The derivative function for this problem returns a vector with $2\mathit{N}$ elements. The first $\mathit{N}$ elements correspond to Burgers' equations, while the last $\mathit{N}$ elements are for the moving mesh equations. The function movingMeshODE goes through these steps to evaluate the right-hand sides of all the equations in the system:

1. Evaluate Burgers' equations using finite differences (first $\mathit{N}$ elements).

2. Evaluate monitor function (last $\mathit{N}$ elements).

3. Apply spatial smoothing to monitor function and evaluate moving mesh equations.

The first $\mathit{N}$ equations in the derivative function encode the right side of Burgers' equations. Burgers' equations can be considered as a differential operator involving spatial derivatives of the form:

$\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)$.

The reference paper [1] describes the process of approximating the differential operator $\mathit{f}$ using centered finite differences by

$${\mathit{f}}_{\mathit{i}}=ϵ\left[\frac{\left(\frac{{\mathit{u}}_{\mathit{i}+1}-{\mathit{u}}_{\mathit{i}}}{{\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}}}\right)-\left(\frac{{\mathit{u}}_{\mathit{i}}-{\mathit{u}}_{\mathit{i}-1}}{{\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{i}-1}}\right)}{\frac{1}{2}\left({\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}-1}\right)}\right]-\frac{1}{2}\left(\frac{{\mathit{u}}_{\mathit{i}+1}^2-{\mathit{u}}_{\mathit{i}-1}^2}{{\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}-1}}\right).$$

On the edges of the mesh (for which $\mathit{i}=1$ and $\mathit{i}=\mathit{N}$), only single-sided differences are used instead. This example uses $ϵ=1\times {10}^{-4}$.

The equations governing the mesh (comprising the last $\mathit{N}$ equations in the derivative function) are

$$\frac{{\partial }^2\dot{\mathit{x}}}{\partial {\mathit{t}}^2}=\frac{1}{\tau }\frac{\partial }{\partial \mathit{t}}\left(\mathit{B}\left(\mathit{x},\mathit{t}\right)\frac{\partial \mathit{x}}{\partial \mathit{t}}\right).$$

Just as with Burgers' equations, you can use finite differences to approximate the monitor function $\mathit{B}\left(\mathit{x},\mathit{t}\right)$:

$$\mathit{B}\left(\mathit{x},\mathit{t}\right)=\sqrt{1+{\left(\frac{\partial {\mathit{u}}_{\mathit{i}}}{\partial {\mathit{x}}_{\mathit{i}}}\right)}^2}=\sqrt{1+{\left(\frac{{\mathit{u}}_{\mathit{i}+1}-{\mathit{u}}_{\mathit{i}-1}}{{\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}-1}}\right)}^2}.$$

Once the monitor function is evaluated, spatial smoothing is applied (Equations 14 and 15 in [1]). This example uses $\gamma =2$ and $\mathit{p}=2$ for the spatial smoothing parameters.

The function encoding the system of equations is

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

Code Functions for Sparsity Patterns

The Jacobian $\mathrm{dF}/\mathrm{dy}$ for the derivative function is a $2\mathit{N}\times 2\mathit{N}$ matrix containing all of the partial derivatives of the derivative function, movingMeshODE. ode15s estimates the Jacobian using finite differences when the matrix is not supplied in the options structure. You can supply the sparsity pattern of the Jacobian to help ode15s calculate it more quickly.

The function for the sparsity pattern of the Jacobian is

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

Plot the sparsity pattern of $\mathrm{dF}/\mathrm{dy}$ for $\mathit{N}=80$ using spy.

spy(JPat(80))

Another way to make the calculation more efficient is to provide the sparsity pattern of $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$. You can find this sparsity pattern by examining which terms of ${\mathit{u}}_{\mathit{i}}$ and ${\mathit{x}}_{\mathit{i}}$ are present in the finite differences calculated in the mass matrix function.

The function for the sparsity pattern of $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$ is

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

Plot the sparsity pattern of $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$ for $\mathit{N}=80$ using spy.

spy(MvPat(80))

Solve System of Equations

Solve the system with the value $\mathit{N}=80$. For the initial conditions, initialize $\mathit{x}$ with a uniform grid and evaluate $\mathit{u}\left(\mathit{x},0\right)$ on the grid.

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

Use odeset to create an options structure that sets several values:

• A function handle for the mass matrix

• The state-dependence of the mass matrix, which for this problem is 'strong' since the mass matrix is a function of both $\mathit{t}$ and $\mathit{y}$

• A function handle that calculates the Jacobian sparsity pattern

• A function handle that calculates the sparsity pattern of the derivative of the mass matrix multiplied by a vector

• The absolute and relative error tolerances

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

Finally, call ode15s to solve the system on the interval $\left[0,1\right]$ using the movingMeshODE derivative function, the time span, the initial conditions, and the options structure.

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

Plot Results

The result of the integration is a structure sol that contains the time steps $\mathit{t}$, the mesh points $\mathit{x}\left(\mathit{t}\right)$, and the solution $\mathit{u}\left(\mathit{x},\mathit{t}\right)$. Extract these values from the structure.

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

Plot the movement of the mesh points over time. The plot shows that the mesh points retain a reasonably even spacing over time (due to the monitor function), but they are able to cluster near the discontinuity in the solution as it moves.

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

Now, sample $\mathit{u}\left(\mathit{x},\mathit{t}\right)$ at a few values of $\mathit{t}$ and plot the evolution of the solution over time. The mesh points at the ends of the interval are fixed, so x(0) = 0 and x(N+1) = 1. The boundary values are u(t,0) = 0 and u(t,1) = 0, which you must add to the known values computed for the figure.

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

The plot shows that $\mathit{u}\left(\mathit{x},0\right)$ is a smooth wave that develops a steep gradient over time as it moves towards $\mathit{x}=1$. The mesh points track the movement of the discontinuity so that extra evaluation points are in the appropriate position in each time step.

References

[1] Huang, Weizhang, et al. “Moving Mesh Methods Based on Moving Mesh Partial Differential Equations.” Journal of Computational Physics, vol. 113, no. 2, Aug. 1994, pp. 279–90. https://doi.org/10.1006/jcph.1994.1135.

Local Functions

Listed here are the local helper functions that the solver ode15s calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

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
% -------------------------------------------------------------------------