不久前，我在https://mathematica.stackexchange.com上写了一些这样的例子。你可以在这里、这里和这里找到它们。最后一个链接可能是最好理解的——即使这不适用于耦合 PDE；但没有太大区别。在 Mathematicas 的函数NDEigensystem的参考页面上也有示例，它可以满足您的需求。

让我在这里复制一个小例子，它更适合可能不熟悉该语言的观众。

为了模拟梁的模式，我们定义了平面应力算子（耦合 PDE）：

planeStress = {
 Div[{{0, -((Y*ν)/(1 -ν^2))}, {-(Y*(1 -ν))/(2*(1 - ν^2)), 0}}.Grad[v[t, x, y], {x, y}], {x, y}] + Div[{{-(Y/(1 - ν^2)),  0}, {0, -(Y*(1 - ν))/(2*(1 - ν^2))}}.Grad[u[t, x, y], {x, y}], {x, y}], 
 Div[{{0, -(Y*(1 - ν))/(2*(1 - ν^2))}, {-((Y*ν)/(1 - ν^2)), 0}}.Grad[u[t, x, y], {x, y}], {x, y}] + Div[{{-(Y*(1 -ν))/(2*(1 -ν^2)), 0}, {0, -(Y/(1 - ν^2))}}.Grad[v[t, x, y], {x, y}], {x, y}]
} /. {Y -> 210 10^9, ν -> 33/100};

我们把左手边的横梁固定住。

(*held fixed at left*)
bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x == 0];

并调用数值微分特征求解器：

\[Rho] = 7850; length = 10; height = 0.7;
{vals, funs} = 
  NDEigensystem[{Thread[\[Rho]*D[{u[t, x, y], v[t, x, y]}, {t, 1}] == 
      planeStress], 
    DirichletCondition[{u[t, x, y] == 0., v[t, x, y] == 0.}, 
     x == 0]}, {u[t, x, y], v[t, x, y]}, 
   t, {x, y} \[Element] Rectangle[{0, 0}, {length, height}], 9];

这将返回特征值和特征函数的列表。绘制它们：

Show[{Graphics3D[
     {Gray, 
      GraphicsComplex[{{0, 0, 0}, {length, 0, 0}, {length, height, 
         0}, {0, height, 0}}, Line[{{1, 2, 3, 4, 1}}]]}],
    Plot3D[
     Sqrt[10000*Total[#^2]], {x, y} \[Element] 
      Rectangle[{0, 0}, {length, height}], 
     ColorFunction -> "TemperatureMap", Axes -> False, 
     Mesh -> False]}, Boxed -> False] & /@ funs

有关其工作原理的更多详细信息，请查看我在开头提到的帖子。