Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

求解具有强状态相关质量矩阵的 ODE

此示例说明如何使用移动网格方法求解 Burgers 方程 [1]。该问题包括质量矩阵,并指定了选项来描述质量矩阵的强状态依赖性和稀疏性,使求解过程更加高效。

问题设置

Burgers 方程是一种对流扩散方程,由以下 PDE 给出

ut=ϵ2ux2-x(u22),0<x<1,t>0,ϵ=1×10-4.

应用坐标变换([1] 中的方程 18)会使左侧出现额外的项:

ut-uxxt=ϵ2ux2-x(u22).

通过使用有限差分逼近关于 x 的偏导数可将 PDE 转换为仅有一个变量的 ODE。如果有限差分写为 Δ,则 PDE 可以重写为只包含关于 t 的导数的 ODE:

dudt-Δudxdt=ϵΔ2u-Δ(u22).

在此形式中,您可以使用 ODE 求解器(如 ode15s)来求解 ux 随时间的变化。

对于此示例,问题由包含 N 个点的移动网格表示,[1] 中所述的移动网格方法在每个时间步定位网格点,使它们始终集中在变化区域。边界和初始条件是

u(0,t)=u(1,t)=0,u(x,0)=sin(2πx)+12sin(πx).

对于由 N 个点组成的给定初始网格,有 2N 个方程需要求解:N 个方程对应于 Burgers 方程,另外 N 个方程确定每个网格点的运动。因此,最终的方程组是:

du1dt-Δu1dx1dt=ϵΔ2u1-Δ(u122),duNdt-ΔuNdxNdt=ϵΔ2uN-Δ(uN22),d2x1˙dt2=1τddt(B(x1,t)dx1dt),d2xN˙dt2=1τddt(B(xN,t)dxNdt).

移动网格的项对应于 [1] 中的 MMPDE6。参数 τ 表示强制网格趋向均匀分布的时间刻度。项 B(x,t) 是由[1] 中的方程 21 给出的监控函数:

B(x,t)=1+(duidxi)2.

此示例中使用的通过移动网格点求解 Burgers 方程的方式展示了几种方法:

  • 方程组用质量矩阵公式表示,M y=f(t,y)。质量矩阵作为函数提供给 ode15s 求解器。

  • 导函数不仅包括表示 Burgers 方程的方程,还包括一组控制移动网格选择的方程。

  • Jacobian dF/dy 的稀疏模式和质量矩阵的导数乘以向量 d(Mv)/dy 作为函数提供给求解器。提供这些稀疏模式有助于求解器更高效地运行。

  • 有限差分用于逼近几个偏导数。

要在 MATLAB® 中求解此方程,请编写导函数、质量矩阵函数、Jacobian dF/dy 的稀疏模式的函数以及 d(Mv)/dy 的稀疏模式的函数。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写质量矩阵代码

方程组的左侧涉及一阶导数的线性组合,因此需要质量矩阵来表示所有项。将方程组的左侧设置为等于 M y 以提取质量矩阵的形式。质量矩阵由四个分块组成,每个分块均为 N 阶方阵:

[u1t-u1x1x1tuNt-uNxNxNt2x1˙t22xN˙t2]=M y=[M1M2M3M4][u1˙uN˙x1˙xN˙]

此公式表明,M1M2 构成 Burgers 方程的左侧(方程组中的前 N 个方程),而 M3M4 构成网格方程的左侧(方程组中的最后 N 个方程)。分块矩阵是:

M1=IN,M2=-uixiIN,M3=0N,M4=2t2IN.

INN×N 单位矩阵。M2 中的偏导数是使用有限差分估计的,而 M4 中的偏导数使用拉普拉斯矩阵。请注意,M3 只包含零,因为网格移动的所有方程都不依赖于 u˙

现在您可以编写一个用于计算质量矩阵的函数。该函数必须接受时间 t 和解向量 y 两个输入。由于解向量 y 包含一半 u˙ 分量和一半 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

注意:所有函数都作为局部函数包含在示例的末尾。

编写导函数代码

此问题的导函数返回一个包含 2N 个元素的向量。前 N 个元素对应于 Burgers 方程,而后 N 个元素对应于移动网格方程。函数 movingMeshODE 通过这些步骤来计算方程组中所有方程的右侧:

  1. 使用有限差分(前 N 个元素)计算 Burgers 方程。

  2. 计算监控函数(最后 N 个元素)。

  3. 应用空间平滑处理来监控函数和计算移动网格方程。

导函数中的前 N 个方程对 Burgers 方程的右侧进行编码。Burgers 方程可视为微分运算子,涉及以下形式的空间导数:

f(u)=ϵ2ux2-x(u22)

参考文献 [1] 通过以下公式说明使用中心有限差分逼近微分运算子 f 的过程

fi=ϵ[(ui+1-uixi+1-xi)-(ui-ui-1xi-xi-1)12(xi+1-xi-1)]-12(ui+12-ui-12xi+1-xi-1).

在网格的边缘(此时 i=1i=N),仅使用单侧差分。此示例使用 ϵ=1×10-4

控制网格的方程(包括导函数中的最后 N 个方程)是

2x˙t2=1τt(B(x,t)xt).

就像对 Burgers 方程一样,您可以使用有限差分来逼近监控函数 B(x,t)

B(x,t)=1+(uixi)2=1+(ui+1-ui-1xi+1-xi-1)2.

计算监控函数后,应用空间平滑处理([1] 中的方程 14 和 15)。此示例使用 γ=2p=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 dF/dy2N×2N 矩阵,其中包含导函数 movingMeshODE 的所有偏导数。当 options 结构体中没有提供矩阵时,ode15s 使用有限差分估计 Jacobian 矩阵。您可以提供 Jacobian 的稀疏模式来帮助 ode15s 更快地计算它。

Jacobian 矩阵稀疏模式的函数是

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

使用 spy 绘制 N=80dF/dy 稀疏模式。

spy(JPat(80))

Figure contains an axes. The axes contains an object of type line.

另一种提高计算效率的方法是提供 d(Mv)/dy 的稀疏模式。通过检查在质量矩阵函数中计算的有限差分中存在 uixi 的哪些项,可以找到这种稀疏模式。

用于计算 d(Mv)/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

使用 spy 绘制 N=80d(Mv)/dy 稀疏模式。

spy(MvPat(80))

Figure contains an axes. The axes contains an object of type line.

求解方程组

用值 N=80 求解方程组。对于初始条件,用均匀网格初始化 x,并在网格上计算 u(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 创建一个 options 结构体,该结构体设置多个值:

  • 质量矩阵的函数句柄

  • 质量矩阵的状态依赖性,对于此问题是 'strong',因为质量矩阵是 ty 两者的函数

  • 用于计算 Jacobian 稀疏模式的函数句柄

  • 用于计算质量矩阵的导数并乘以向量的稀疏模式的函数句柄

  • 绝对和相对误差容限

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

最后,调用 ode15s 以使用 movingMeshODE 导函数、时间跨度、初始条件和 options 结构体在区间 [0, 1] 上求解方程组。

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

绘制结果

积分的结果是结构体 sol,其中包含时间步 t、网格点 x(t) 和解 u(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')

Figure contains an axes. The axes with title Burgers' equation: Trajectories of grid points contains 80 objects of type line.

现在,在 t 的几个值上对 u(x,t) 进行采样,并绘制解随时间的演变图。区间末端的网格点是固定的,即 x(0) = 0x(N+1) = 1。边界值为 u(t,0) = 0u(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

Figure contains an axes. The axes with title Burgers equation on moving mesh contains 6 objects of type line. These objects represent t = 0, t = 0.2, t = 0.4, t = 0.6, t = 0.8, t = 1.

该图显示 u(x,0) 是平滑的波形,随着时间的推移,它会向 x=1 移动,形成陡峭的梯度。网格点跟踪不连续点的移动,以便在每个时间步中额外的计算点位于适当的位置。

参考

[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.

局部函数

此处列出了求解器 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
% -------------------------------------------------------------------------

另请参阅

|

相关主题