2024-11-30 |欢迎您，第位读者！|

Wave_PDE_In_Matlab求解时变波动方程

波动方程

这都很正经波传播结果，不要多想……不要乱想……这些个波都很正经。

网格也很正经。

这是在一个心脏模型上求解时变波动方程的结果。

$$\frac{\partial ^2}{\partial t^2} u\left(t,x,y\right) =\frac{\partial ^2}{\partial x^2} u\left(t,x,y\right)+\frac{\partial ^2}{\partial y^2} u\left(t,x,y\right)$$

问题的边界按照网格图上分为了四个部分，都能够设置边界条件，当然这里是二阶方程，所以设置边界条件是必须设置$u$和$\frac{du}{dt}$。当然作为一个时变方程，还有初值也需要设置。

该怎么求解呢？

PDE工具箱系数推导

偏微分方程的散度形式

偏微分工具箱（Partial Differential equation Toolbox $^{\text{TM}}$）求解的方程形式为：

$$ \mathbf{m} \frac{\partial^2 \mathbf{u}}{\partial t^2} + \mathbf{d} \frac{\partial \mathbf{u}}{\partial t} - \nabla \cdot (\mathbf{c} \otimes \nabla \mathbf{u}) + \mathbf{a} \mathbf{u} = \mathbf{f} $$

其中，$\mathbf{u}$ 是未知函数，$\mathbf{m}$ 是质量矩阵，$\mathbf{d}$ 是阻尼矩阵，$\mathbf{c}$ 是扩散矩阵，$\mathbf{a}$ 是刚度矩阵，$\mathbf{f}$ 是外力矩阵。

当一个偏微分方程，可以写成上述散度形式（Divergence Form）,即可以采用有限元方法很方便的求解，这里有一个条件就是方程的系数矩阵必须不包含函数的偏微分（可以包括函数、坐标、时间等变量）。

在PDE工具箱中，把方程组描述为散度形式并设定对应的系数矩阵，然后设置边界、初始条件并通过求解器求解。

1model = createpde()
2% ...
3specifyCoefficients(model, m=m, d=d, c=c, a=a, f=f)
4% ...

所以，通常需要进行的一个问题就是，如何把一个偏微分方程、方程组转化为这样的形式。

当然，数学达人很容易通过变换和推导，有些时候需要求解一些代数方程来得到系数矩阵。但是这对于一般的工程师来说，可能就有些困难了。

这个时候就有一个很好用的Matlab工具箱能够派上用场：符号计算工具箱。

符号计算工具箱提供了两个函数：pdeCoefficents和pdeCoefficientsToDouble来完成这个转化。得到对应的系数矩阵后，就能够调用PDE工具箱进行求解。这个功能需要2021a之后的版本，如果版本更低，只好升级。

推导与求解

这个函数pdfCeofficents非常简单。

 1>> help pdeCoefficients
 2--- sym/pdeCoefficients 的帮助 ---
 3
 4  pdeCoefficients Extract the coefficients of a pde.
 5     S = pdeCoefficients(EQ, U) extracts the coefficients
 6     of a symbolic expression EQ that represents a PDE with the dependent
 7     variable U. The result is a struct S that can be used in the function
 8     specifyCoefficients of the PDE toolbox.
 9 
10     S = pdeCoefficients(EQ, U, 'Symbolic', true)
11     extracts the coefficients and returns them as a struct S of
12     symbolic expressions.
13 
14    Example:
15        syms u(x, y)
16        S = pdeCoefficients(laplacian(u) - x, [u]);
17        model = createpde();
18        geometryFromEdges(model,@lshapeg);
19        specifyCoefficients(model, S)
20 
21        creates a pde model for the equation laplacian(u) == x.
22        This example requires the PDE Toolbox.
23 
24    See also pdeCoefficientsToDouble, CREATEPDE, pde.PDEModel.
25
26    sym/pdeCoefficients 的文档
27       doc sym/pdeCoefficients

这里的例子也给得很清楚：

1syms u(x, y)
2S = pdeCoefficients(laplacian(u) - x, [u]);
3model = createpde();
4geometryFromEdges(model,@lshapeg);
5specifyCoefficients(model, S)

首先，要把被积函数的函数关系定义出来，作为符号表达式。接着，就给出被积函数的整体形式作为函数的第一个参数，并把因变量，只有一个可以省略[]，作为第二个参数。

函数返回的结果是一个结构体，这个结构体可以直接作为specifyCoefficients的参数，当然，这个形式我们知道，跟分开写key=value或者，'key', value,...的调用形式是一样的，为了更加清晰，甚至可以故意写成：

1specifyCoefficients(model, 'm', S.m, 'a', S.a, 'c', S.c, 'f', S.f, 'd', S.d)

其他，都可以不用仔细说明。

当然，如果这个推导结果表明偏微分方程无法写成散度形式时，会提示：

1Warning: After extracting m, d, and c, some gradients remain. Writing all remaining terms to f.

这个时候，可以通过下面的方式显示表达式f的具体形式，因为转换后剩余的部分，全部都被归到f这个项中。当f还有偏导项时，PDE工具箱就不能求解这个方程。

1S.f('show')

全部代码

最后，完整给出代码。

wave

代码中唯一个好玩一点点的就是最后的gif生成，当目录下已经有文件的时候，就增加一个编号，这样保证gif文件是新的，当exportgraphics采用Append=true的形式调用时，就会把图形增加到后面，如果是一个已有的gif，就会发生不希望的情况。

最后就是，固定ColorMap增加图像动画的一致性，可以通过hsv(n)中的n来调节图形的精细程度。

 1fprintf("Model and mesh...\n")
 2model = createpde;
 3geometryFromEdges(model, @cardg);
 4generateMesh(model, "Hmax", 0.1);
 5figure(1)
 6clf;
 7pdegplot(model, 'EdgeLabels', 'on')
 8hold on
 9pdemesh(model)
10exportgraphics(gca, 'MeshWithEdgeLabels.png');
11
12fprintf("Setting up the PDE...\n");
13syms u(t, x, y)
14pdeeq = diff(u, t, t) - laplacian(u, [x, y])
15
16
17coeffs = pdeCoefficients(pdeeq, u)
18
19specifyCoefficients(model, 'm', coeffs.m, 'd', coeffs.d, 'c', coeffs.c, 'a', coeffs.a, 'f', coeffs.f);
20
21fprintf("Setting up the boundary conditions...\n");
22setInitialConditions(model, 0, 0);
23% setInitialConditions(model, 1, 0, 'Edge', 2:3)
24applyBoundaryCondition(model, 'dirichlet', 'Edge', [1], 'u', 1);
25
26fprintf("Solving the PDE...\n");
27results = solvepde(model, linspace(0, 10, 50));
28fprintf("Done\n");
29
30%% Visualize the solution
31figure(3);
32
33fnk = @(k)sprintf("./wave-%d.gif", k);
34k = 1;
35fn = fnk(k);
36
37while exist(fn, 'file')
38    k = k + 1;
39    fn = fnk(k);
40    % delete(fn);
41end
42
43fprintf("Visualized to file %s...\n", fn)
44
45for idx = 1:50
46    pdeplot(model, 'XYData', results.NodalSolution(:, idx), 'ColorMap', hsv(40));
47    title(['t = ', num2str(results.SolutionTimes(idx))]);
48    axis equal
49    axis off
50    colorbar off
51    clim([-1, 3])
52    exportgraphics(gca, fn, "append", true);
53end