这个例子展示了如何使用移动网格技术[1]求解Burgers方程。该问题包括一个质量矩阵,并指定选项,以说明质量矩阵的强状态依赖性和稀疏性,使求解过程更有效。
Burgers方程是偏微分方程给出的对流扩散方程
应用坐标变换([1]中的Eq. 18)会导致左边有一个额外的项:
将PDE转换为一个变量的ODE是通过使用有限差分来近似于的偏导数来完成的 .如果有限差分可以写成 ,则PDE可以重写为只包含对求导数的ODE :
在这种形式中,您可以使用ODE求解器,例如ode15s
来解
而且
随着时间的推移。
对于这个例子,这个问题是在移动网的 点,[1]中描述的移动网格技术在每个时间步骤中定位网格点,以便它们集中在变化的区域。边界条件和初始条件为
对于给定的初始网格 点,有 需要求解的方程: 与Burgers方程相对应的方程,和 确定每个网格点运动的方程。所以,最终的方程组是
移动网格的术语对应于[1]中的MMPDE6。的参数 表示迫使网格向均等分布的时间刻度。这个词 为[1]中Eq. 21给出的监控函数:
本例中使用的方法来求解带有移动网格点的Burgers方程,演示了几种技术:
方程组用质量矩阵公式表示,
.质量矩阵提供给ode15s
解算器是一个函数。
导数函数不仅包括Burgers方程的方程,还包括一组控制运动网格选择的方程。
雅可比矩阵的稀疏模式 质量矩阵的导数乘以一个向量 作为函数提供给求解器。提供这些稀疏模式有助于求解器更有效地运行。
有限的差异用来近似几个偏导数。
要在MATLAB®中求解这个方程,请编写一个导数函数,一个质量矩阵函数,一个雅可比矩阵的稀疏模式函数 的稀疏模式函数 .您可以将所需的函数作为本地函数包含在文件的末尾(如这里所做的),或者将它们保存为单独的命名文件,保存在MATLAB路径的目录中。
方程组的左边涉及一阶导数的线性组合,因此需要一个质量矩阵来表示所有的项。令方程组左边等于 求出质量矩阵的形式。质量矩阵由四个方块组成,每个方块都是一个有序的方阵 :
这个公式表明 而且 形成Burgers方程的左边(第一个 方程组中的方程),而 而且 形成左边的网格方程(最后 方程组中的方程)。分块矩阵为:
是 单位矩阵。的偏导数 都是用有限差分估计的,而偏导数在 使用拉普拉斯矩阵。请注意, 只包含零,因为网格运动的方程都不依赖于 .
现在你可以写一个函数来计算质量矩阵。该函数必须接受两个时间输入 解向量 .因为解向量 包含了一半 成分和一半 组件,函数首先提取这些组件。然后,该函数形成所有的块矩阵(考虑问题的边界值),并使用四个块组装质量矩阵。
函数M =质量(t,y)提取解u的y分量,并对x进行网格划分N =长度(y)/2;u = y(1:N);x = y(N+1:结束);%解u与网格x的边界值U0 = 0;uNP1 = 0;X0 = 0;xNP1 = 1;% M1和M2是Burgers方程的质量矩阵的部分。du/dx的导数近似于有限差分,使用%单面差异的边缘和中心之间的差异。M1 = speye(N);M2 =稀疏(N,N);M2(1,1) = - (u(2) - u0)/(x(2) - x0);为我= 2:n - 1平方米(我)= - (u (i + 1) - u(张))/ (x (i + 1) - x(张));结束M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));% M3和M4定义了网格点演化方程,对应于% MMPDE6在参考论文。由于网格函数只涉及d/dt(dx/dt),质量矩阵的M3部分全为零。M4的二阶导数是%近似使用有限差分拉普拉斯矩阵。M3 =稀疏(N,N);e = ones(N,1);M4 = spdigs ([e -2*e],-1:1,N,N);组装质量矩阵M = [m1 m2 m3 m4];结束
注意:所有函数都作为局部函数包含在示例的末尾。
这个问题的导数函数返回一个向量
元素。第一个
元素对应于Burgers方程,而最后一个元素对应于Burgers方程
元素是用于移动网格方程的。这个函数movingMeshODE
通过这些步骤来计算系统中所有方程的右边:
用有限差分来评估Burgers方程(首先 元素)。
评估监视器功能(最后 元素)。
应用空间平滑来监测函数和评估移动网格方程。
第一个 导数函数中的方程编码Burgers方程的右侧。Burgers方程可以被认为是包含如下形式的空间导数的微分算子:
.
参考文献[1]描述了微分算子的逼近过程 用中心有限差分法
在网格的边缘(为其中 而且 ),只使用单面差异。这个例子使用了 .
控制网格的方程(包括最后一个 导数函数中的方程)是
与Burgers方程一样,可以使用有限差分来近似监测函数 :
一旦评估了监测函数,就应用空间平滑([1]中的公式14和15)。这个例子使用了 而且 为空间平滑参数。
编码方程组的函数是
函数g = movingMeshODE(t,y)提取解u的y分量,并对x进行网格划分N =长度(y)/2;u = y(1:N);x = y(N+1:结束);%解u与网格x的边界值U0 = 0;uNP1 = 0;X0 = 0;xNP1 = 1;预分配g矢量的导数值。g = 0 (2*N,1);用中心有限差分近似Burgers的RHS%方程(边有单面差异)。第一个Ng中的%元素对应于Burgers方程。为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)²- u(i-1)²)/delx;结束Delx = x(2) - x0;g(1) = 1的军医* ((u (2) - u (1)) / (x (2) - (1) - (u(1) -情况)/ (x (1) - x0)) / (0.5 * delx)...- 0.5*(u(2)²- u0²)/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²- u(N-1)²)/delx;计算监测函数值(参考文献中公式21),用于网格方程的% RHS。内部采用中心有限差分%点,边沿使用单面差异。M = 0 (N,1);为我= 2:n - 1 M (i) =√6 (1 + ((u (i + 1) - u(张))/ (x (i + 1) - x(张)))^ 2);结束M0 =√(1 + (u(1) - u0)/(x(1) - x0) ^2);M(1) =√(1 + ((u(2) - u0)/(x(2) - x0) ^2);M (N =√(1 + (uNP1 - u (N - 1) / (xNP1 - x (N - 1))) ^ 2);MNP1 =√(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);应用空间平滑(公式14和15),gamma = 2, p = 2。SM = 0 (N,1);为我= 3:n - SM (i) =√(4 * M(我2)^ 2 + 6 * M(张)^ 2 + 9 * M (i) ^ 2 +...6*M(i+1)²+ 4*M(i+2)²)/29;结束SM0 =√(9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);SM(1) =√(M0 ^ 2 + 6 * 9 * M (1) ^ 2 + 6 * M (2) ^ 2 + 4 * M (3) ^ 2) / 25);SM(2) =√(4 * M0 ^ 2 + 6 * M (1) ^ 2 + 9 * M (2) ^ 2 + 6 * M (3) ^ 2 + 4 * M (4) ^ 2) / 29);SM (N - 1) =√(4 * M (N) ^ 2 + 6 * M (N - 2) ^ 2 + 9 * M (N - 1) ^ 2 + 6 * M (N) ^ 2 + 4 * MNP1 ^ 2) / 29);SM (N =√(4 * M (N - 2) ^ 2 + 6 * M (N - 1) ^ 2 + 9 * M (N) ^ 2 + 6 * MNP1 ^ 2) / 25);SMNP1 =√(4 * M (N - 1) ^ 2 + 6 * M (N) ^ 2 + 9 * MNP1 ^ 2) / 19);为我= 2:N - 1 g (i + N) = (SM (i + 1) + SM(我))* (x (i + 1) - x(我))-...(SM(i) + SM(i-1))*(x(i) - x(i-1));结束g (1 + N) = (SM (2) + SM (1)) * (x (2) - (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));形成参考论文中方程12的最终离散近似,即控制方程%网格点。Tau = 1e-3;g(1+N:end) = - g(1+N:end)/(2*tau);结束
雅可比矩阵
导数函数是a
矩阵包含了导数函数的所有偏导数,movingMeshODE
.ode15s
当矩阵没有在选项结构中提供时,使用有限差分估计雅可比矩阵。你可以提供雅可比矩阵的稀疏模式ode15s
更快地计算。
雅可比矩阵的稀疏图的函数是
函数out = JPat(N) S1 = spdigs (ones(N,3),-1:1,N,N);S2 = spdigs (ones(N,9),-4:4,N,N);out = [S1 S1 S2 S2];结束
的稀疏图
为
使用间谍
.
间谍(JPat (80))
另一种使计算更有效的方法是提供的稀疏模式 .您可以通过检查的项来找到这种稀疏模式 而且 存在于质量矩阵函数计算的有限差分中。
的稀疏模式的函数 是
函数S = MvPat(N) S =稀疏(2*N,2*N);S(1,2) = 1;S(1,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;结束S(n, n -1) = 1;S(n, n -1+ n) = 1;结束
的稀疏图
为
使用间谍
.
间谍(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];
使用odeset
创建一个选项结构,设置几个值:
质量矩阵的函数句柄
质量矩阵的状态依赖关系,对于这个问题“强”
因为质量矩阵是两者的函数
而且
计算雅可比稀疏模式的函数句柄
一个函数句柄,用于计算质量矩阵的导数乘以矢量的稀疏性模式
绝对误差和相对误差容差
Opts = odeset(“质量”@mass,“MStateDependence”,“强”,“JPattern”JPat (N),...“MvPattern”MvPat (N),“RelTol”1 e-5“AbsTol”1的军医);
最后,调用ode15s
在区间上解方程组
使用movingMeshODE
导数函数,时间跨度,初始条件和期权结构。
Tspan = [0 1];sol = ode15s(@movingMeshODE,tspan,y0,opts);
整合的结果是一个结构索尔
它包含了时间步长
,网格点
,以及解决方案
.从结构中提取这些值。
T = sol.x;x = sol.y(N+1:end,:);u = sol.y(1:N,:);
绘制网格点随时间的运动。该图显示网格点随着时间的推移保持合理的均匀间距(由于监测功能),但它们能够在解决方案移动时聚集在不连续的附近。
情节(x, t)包含(“t”) ylabel (“x (t)”)标题(汉堡方程:网格点轨迹)
现在,示例
的一些值
并绘制出解决方案随时间的演变。区间两端的网格点是固定的,所以X (0) = 0
而且x(N+1) = 1
.边界值为U (t,0) = 0
而且U (t,1) = 0
,必须将其添加到为该图计算的已知值中。
色调= 0:0.2:1;Yint = deval(sol,色调);图标签= {};为J = 1:长度(着色)解= [0;yint (1: N, j);0);位置= [0;yint (N + 1:结束,j);1);标签{j} = ['t = 'num2str(色(j)));情节(位置、解决方案“o”)举行在结束包含(“x”) ylabel (“解决方案u (x, t)”{:})传说(标签,“位置”,“西南”)标题(“运动网格上的Burgers方程”)举行从
图表显示 平滑波在移动的过程中会随着时间的推移产生陡峭的梯度吗 .网格点跟踪不连续性的运动,使额外的评价点在每个时间步中处于适当的位置。
[1]黄伟章等“基于动网格偏微分方程的动网格方法”。计算物理杂志,第113卷,no。2, 1994年8月,第279-90页。https://doi.org/10.1006/jcph.1994.1135.
这里列出的是局部帮助函数ode15s
调用来计算解决方案。或者,您可以将这些函数保存为它们自己的文件,保存在MATLAB路径的目录中。
函数g = movingMeshODE(t,y)提取解u的y分量,并对x进行网格划分N =长度(y)/2;u = y(1:N);x = y(N+1:结束);%解u与网格x的边界值U0 = 0;uNP1 = 0;X0 = 0;xNP1 = 1;预分配g矢量的导数值。g = 0 (2*N,1);用中心有限差分近似Burgers的RHS%方程(边有单面差异)。第一个Ng中的%元素对应于Burgers方程。为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)²- u(i-1)²)/delx;结束Delx = x(2) - x0;g(1) = 1的军医* ((u (2) - u (1)) / (x (2) - (1) - (u(1) -情况)/ (x (1) - x0)) / (0.5 * delx)...- 0.5*(u(2)²- u0²)/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²- u(N-1)²)/delx;计算监测函数值(参考文献中公式21),用于网格方程的% RHS。内部采用中心有限差分%点,边沿使用单面差异。M = 0 (N,1);为我= 2:n - 1 M (i) =√6 (1 + ((u (i + 1) - u(张))/ (x (i + 1) - x(张)))^ 2);结束M0 =√(1 + (u(1) - u0)/(x(1) - x0) ^2);M(1) =√(1 + ((u(2) - u0)/(x(2) - x0) ^2);M (N =√(1 + (uNP1 - u (N - 1) / (xNP1 - x (N - 1))) ^ 2);MNP1 =√(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);应用空间平滑(公式14和15),gamma = 2, p = 2。SM = 0 (N,1);为我= 3:n - SM (i) =√(4 * M(我2)^ 2 + 6 * M(张)^ 2 + 9 * M (i) ^ 2 +...6*M(i+1)²+ 4*M(i+2)²)/29;结束SM0 =√(9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);SM(1) =√(M0 ^ 2 + 6 * 9 * M (1) ^ 2 + 6 * M (2) ^ 2 + 4 * M (3) ^ 2) / 25);SM(2) =√(4 * M0 ^ 2 + 6 * M (1) ^ 2 + 9 * M (2) ^ 2 + 6 * M (3) ^ 2 + 4 * M (4) ^ 2) / 29);SM (N - 1) =√(4 * M (N) ^ 2 + 6 * M (N - 2) ^ 2 + 9 * M (N - 1) ^ 2 + 6 * M (N) ^ 2 + 4 * MNP1 ^ 2) / 29);SM (N =√(4 * M (N - 2) ^ 2 + 6 * M (N - 1) ^ 2 + 9 * M (N) ^ 2 + 6 * MNP1 ^ 2) / 25);SMNP1 =√(4 * M (N - 1) ^ 2 + 6 * M (N) ^ 2 + 9 * MNP1 ^ 2) / 19);为我= 2:N - 1 g (i + N) = (SM (i + 1) + SM(我))* (x (i + 1) - x(我))-...(SM(i) + SM(i-1))*(x(i) - x(i-1));结束g (1 + N) = (SM (2) + SM (1)) * (x (2) - (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));形成参考论文中方程12的最终离散近似,即控制方程%网格点。Tau = 1e-3;g(1+N:end) = - g(1+N:end)/(2*tau);结束% -----------------------------------------------------------------------函数M =质量(t,y)提取解u的y分量,并对x进行网格划分N =长度(y)/2;u = y(1:N);x = y(N+1:结束);%解u与网格x的边界值U0 = 0;uNP1 = 0;X0 = 0;xNP1 = 1;% M1和M2是Burgers方程的质量矩阵的部分。du/dx的导数近似于有限差分,使用%单面差异的边缘和中心之间的差异。M1 = speye(N);M2 =稀疏(N,N);M2(1,1) = - (u(2) - u0)/(x(2) - x0);为我= 2:n - 1平方米(我)= - (u (i + 1) - u(张))/ (x (i + 1) - x(张));结束M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));% M3和M4定义了网格点演化方程,对应于% MMPDE6在参考论文。由于网格函数只涉及d/dt(dx/dt),质量矩阵的M3部分全为零。M4的二阶导数是%近似使用有限差分拉普拉斯矩阵。M3 =稀疏(N,N);e = ones(N,1);M4 = spdigs ([e -2*e],-1:1,N,N);组装质量矩阵M = [m1 m2 m3 m4];结束% -------------------------------------------------------------------------函数out = JPat(N)%雅可比稀疏图S1 = spdigs (ones(N,3),-1:1,N,N);S2 = spdigs (ones(N,9),-4:4,N,N);out = [S1 S1 S2 S2];结束% -------------------------------------------------------------------------函数S = MvPat(N)质量矩阵乘以矢量的导数的%稀疏模式S =稀疏(2*N,2*N);S(1,2) = 1;S(1,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;结束S(n, n -1) = 1;S(n, n -1+ n) = 1;结束% -------------------------------------------------------------------------