# 求解具有强状态依赖质量矩阵的 ODE

### 问题设置

`$\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),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}0<\mathit{x}<1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathit{t}>0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}ϵ=1{\text{\hspace{0.17em}}×\mathrm{10}}^{-4}.$`

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

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

`$\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)=\mathrm{sin}\left(2\pi \mathit{x}\right)+\frac{1}{2}\mathrm{sin}\left(\pi \mathit{x}\right).\end{array}$`

`$\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}\stackrel{˙}{{\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}\stackrel{˙}{{\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}$`

`$\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}}.$`

• 方程组用质量矩阵公式表示，$\mathit{M}\text{}{\mathit{y}}^{\prime }=\mathit{f}\left(\mathit{t},\mathit{y}\right)$。质量矩阵作为函数提供给 `ode15s` 求解器。

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

• 雅可比 $\mathrm{dF}/\mathrm{dy}$ 的稀疏模式和质量矩阵的导数乘以向量 $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$ 作为函数提供给求解器。提供这些稀疏模式有助于求解器更高效地运行。

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

### 编写质量矩阵代码

$\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}\stackrel{˙}{{\mathit{x}}_{1}}}{\partial {\mathit{t}}^{2}}\\ ⋮\\ \frac{{\partial }^{2}\stackrel{˙}{{\mathit{x}}_{\mathit{N}}}}{\partial {\mathit{t}}^{2}}\end{array}\right]=\mathit{M}\text{}{\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}\stackrel{˙}{{\mathit{u}}_{1}}\\ ⋮\\ \stackrel{˙}{{\mathit{u}}_{\mathit{N}}}\\ \stackrel{˙}{{\mathit{x}}_{1}}\\ ⋮\\ \stackrel{˙}{{\mathit{x}}_{\mathit{N}}}\end{array}\right]$.

`$\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}}$$\mathit{N}×\mathit{N}$ 单位矩阵。${\mathit{M}}_{2}$ 中的偏导数是使用有限差分估计的，而 ${\mathit{M}}_{4}$ 中的偏导数使用拉普拉斯矩阵。请注意，${\mathit{M}}_{3}$ 只包含零，因为网格移动的所有方程都不依赖于 $\stackrel{˙}{\mathit{u}}$

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

### 编写导函数代码

1. 使用有限差分（前 $\mathit{N}$ 个元素）计算伯格斯方程。

2. 计算监控函数（最后 $\mathit{N}$ 个元素）。

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

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

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

`$\frac{{\partial }^{2}\stackrel{˙}{\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).$`

`$\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}}.$`

```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 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(JPat(80))`

```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(MvPat(80))`

### 求解方程组

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

• 质量矩阵的函数句柄

• 质量矩阵的状态依赖性，对于此问题是 `'strong'`，因为质量矩阵是 $\mathit{t}$$\mathit{y}$ 两者的函数

• 用于计算雅可比稀疏模式的函数句柄

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

• 绝对和相对误差容限

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

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

### 绘制结果

```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')```

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

### 参考

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

### 局部函数

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