拉普拉斯算子的特征值

打开生活的脚本

这个例子展示了如何解决L形区域上拉普拉斯算子的特征值问题。

膜问题

考虑一个固定在边界上的膜 $\partial Ω$ 的一个地区 $Ω$ 在平面上。它的位移 $u （ x ， y ）$ 是由特征值问题描述的 $Δ u ＝ λ u$ 哪里 $Δ u ＝ u_{x x} + u_{y y}$ 是拉普拉斯算子和吗 $λ$ 是标量参数。边界条件为 $u （ x ， y ）＝ 0$ 总的来说 $（ x ， y ） \in \partial Ω$ ．

拉普拉斯算子是自伴和负定的，也就是说，只有实负特征值 $λ$ 存在。有一个极大的(负的)离散特征值，对应的特征函数 $u$ 被称为基态．在这个例子中, $Ω$ 为L形区域，与该区域相关的基态为L形膜，即MATLAB®徽标。

九点有限差分近似

特征值问题的最简单方法是近似拉普拉斯方程 $Δ u$ 通过有限差分近似(a模板)在有距离的点组成的正方形网格上hx在 $x$ 方向和距离沪元在 $y$ 方向。在这个例子中，是近似的 $Δ u$ 有一笔钱S_h围绕中点的九个常规栅格点的 $（ x ， y ）$ ．未知量就是权重 ${一个}_{- 1 - 1} ，．.. ， {一个}_{11}$ ．

符号u（x，y）每股收益a11a10a1_1a01a00a0_1a_11a_10a_1_1符号hx沪元积极的S_h = a_11 * u(x - Eps*hx,y + Eps*hy) +．..a01*u（x，y+Eps*hy）+．..a11*u（x+Eps*hx，y+Eps*hy）+．..a_10*u（x-Eps*hx，y）+．..a00*u（x，y）+．..(b + b) b + b (c + b) b．..a_1_1*u（x-Eps*hx，y-Eps*hy）+．..a0_1*u（x，y-每股收益*hy）+．..a1_1*u（x+Eps*hx，y-Eps*hy）；

使用符号参数每股收益按的幂对这个表达式的展开进行排序hx和沪元．知道了权重，你可以通过设置来近似拉普拉斯矩阵每股收益= 1．

t =泰勒(S_h, Eps，“秩序”, 7);

使用多项式系数函数来提取幂相同的项的系数每股收益．每个系数都是包含的幂的表达式hx，沪元的导数u关于 $x$ 和 $y$ ．自S_h代表 $u_{x x} + u_{y y}$ ，其他导数的系数u必须为零。通过替换的所有导数提取系数u,除了 $u_{x x}$ 和 $u_{y y}$ , 0。取代 $u_{x x}$ 和 $u_{y y}$ 1。这将泰勒展开式简化为您想要计算的系数，并引出以下六个线性方程。

C=公式（系数（t，Eps，“所有”))（7）u（x，y，y，y），1）==0；eq11=sub（C（6），[dif（u，x，x，x，（u，y，y）及[1，0，（1，0）0，[1，0，（0）eq0=0；eq0=0；EQQ1=0；EQ1=0；EQ1=0=0；eq11=0=0（C（C（C（6），[C（C（6），[C（6），[dif（u，x，（u，x，（u，x，x，x，x，x，x，x，x，x，x，x，x，（u，x，（u，x，x，（u，x，x，x，x，x，x），u，u，x），u，y，y）（u，（u，y，y，u，y）及y）]，，，[u，y，（u，y，（u，（u，y，x，（u，y，y 1,0]）==0；eq23=subs（C（5），[diff（u，x，x），diff（u，x，y），diff（u，y，y）]，[0,0,1]）==1；

因为里面有九个未知的重量S_h，通过要求所有的三阶导数u都是0。

eq31 =潜艇(C (4), (diff (u, x, x, x), diff (u, x, x, y), diff (u, x, y, y) diff (u, y, y, y)], [1, 0, 0, 0) = = 0;eq32 =潜艇(C (4), (diff (u, x, x, x), diff (u, x, x, y), diff (u, x, y, y) diff (u, y, y, y)], [0 1 0,0]) = = 0;eq33 =潜艇(C (4), (diff (u, x, x, x), diff (u, x, x, y), diff (u, x, y, y) diff (u, y, y, y)], [0, 0, 1, 0]) = = 0;eq34 =潜艇(C (4), (diff (u, x, x, x), diff (u, x, x, y), diff (u, x, y, y) diff (u, y, y, y)], [0, 0, 0, 1]) = = 0;

解得到的10个方程，得到9个未知权重。使用ReturnConditions求包含任意参数的所有解。万博尤文图斯

(a11 a10, a1_1、a01 a00、a0_1, a_11, a_10, a_1_1,参数、条件)=．..解决([eq0, eq11、eq12 eq21, eq22, eq23, eq31, eq32, eq33, eq34),．..(a11 a10, a1_1、a01 a00、a0_1, a_11, a_10, a_1_1),．..“返回条件”,真正的);扩大([a_11 a01, a11;．..a_10，a00，a01；．..a1_1，a0_1，a_1_1]）

ans =
                 
                   $（ \begin{array}{c} z & \frac{1}{{沪元}^{2}} - 2 z & z \\ \frac{1}{{hx}^{2}} - 2 z & 4 z - \frac{2}{{hx}^{2}} - \frac{2}{{沪元}^{2}} & \frac{1}{{沪元}^{2}} - 2 z \\ z & \frac{1}{{沪元}^{2}} - 2 z & z \end{array} ）$

参数

参数=
                  $z$

使用潜艇函数将权重替换为其计算值。

C=简化（subs（C））；

表情C (7)，C (6),C (4)包含的0阶，1阶和3阶导数u消失。

[C(7)、C(6)、C (4))

ans =
                  $（ \begin{array}{c} 0 & 0 & 0 \end{array} ）$

表达式C (5)是拉普拉斯的吗u．

C (5)

ans =
                 
                   $\frac{\partial^{2}}{\partial x^{2}} u （ x ， y ） + \frac{\partial^{2}}{\partial y^{2}} u （ x ， y ）$

因此，根据上面计算的权重值，模板S_h近似拉普拉斯阶hx ^ 2，hy ^ 2对于任意参数的任意值z,前提是z是有秩序的吗O（1/hx^2,1/hy^2）．

包含四阶和高阶导数的项

虽然解决方案包含一个自由参数z,该表达式C (3)包含的四阶导数u不能通过选择合适的z．另一个选择是把它变成拉普拉斯算子的平方的倍数。

符号d拉普拉斯= @(u)拉普拉斯(u，[x,y])扩大(d *拉普拉斯(拉普拉斯(u)))

ans（x，y）=
                 
                   $d \frac{\partial^{4}}{\partial x^{4}} u （ x ， y ） + 2 d \frac{\partial^{2}}{\partial y^{2}} \frac{\partial^{2}}{\partial x^{2}} u （ x ， y ） + d \frac{\partial^{4}}{\partial y^{4}} u （ x ， y ）$

取不同的导数u在C (3)，并使它们的系数与相应的项相等。

潜艇(C (3), [diff (u, x, x, x, x), diff (u, x, x, y, y), diff (u, y, y, y, y)], [1, 0, 0)) = = d

ans =
                 
                   $\frac{{hx}^{2}}{12} ＝ d$

潜艇(C (3), [diff (u, x, x, x, x), diff (u, x, x, y, y), diff (u, y, y, y, y)], [0, 1,0]) = = 2 * d

ans =
                  ${hx}^{2} {沪元}^{2} z ＝ 2 d$

subs（C（3），[diff（u，x，x，x，x），diff（u，x，x，y，y），diff（u，y，y，y）]，[0,0,1]）==d

ans =
                 
                   $\frac{{沪元}^{2}}{12} ＝ d$

因此，你可以选择d=hx^2/12=hy^2/12和z = 2 * d / (hx hy ^ ^ 2 * 2),这意味着hx =沪元和z = 1 / (6 * hx *衔接)．因此,模板S_h在正方形网格上近似一个改进的拉普拉斯算子Hx = hy = h．

${年代}_{h} ＝ Δ u + \frac{h^{2}}{12} Δ^{2} u + O （ h^{3.} ）．（ 1 ）$

符号hhx = h;hy = h;d = h ^ 2/12;

取代hx和沪元通过h．

C=subs（C）；

取代z的价值,1/（6*h^2）．因为z在MATLAB®工作空间中不存在，您只能将其作为存储在参数数组中。

C=subs（C，参数，1/（6*h^2））；

验证公式（1）。

简化（C（3）-d*Laplace（Laplace（u）））

ans（x，y）=
                  $0$

现在考虑三阶项hx，沪元．

简化(C (2))

ans =
                  $0$

由于在模板的扩展中不存在这样的术语，该术语 $O （ h^{3.} ）$ 在(1)是事实上的秩序 $O （ h^{4} ）$ ．考虑模具的四阶项。

因子(简化(C (1)))

ans =
                 
                   $（ \begin{array}{c} \frac{1}{360} & h & h & h & h & \frac{\partial^{6}}{\partial x^{6}} u （ x ， y ） + 5 \frac{\partial^{2}}{\partial y^{2}} \frac{\partial^{4}}{\partial x^{4}} u （ x ， y ） + 5 \frac{\partial^{4}}{\partial y^{4}} \frac{\partial^{2}}{\partial x^{2}} u （ x ， y ） + \frac{\partial^{6}}{\partial y^{6}} u （ x ， y ） \end{array} ）$

检查这些项是否可以和拉普拉斯算子的另一次幂联系起来。然而,与

拉普拉斯(拉普拉斯(拉普拉斯(u)))

ans（x，y）=
                 
                   $\frac{\partial^{6}}{\partial x^{6}} u （ x ， y ） + 3. \frac{\partial^{2}}{\partial y^{2}} \frac{\partial^{4}}{\partial x^{4}} u （ x ， y ） + 3. \frac{\partial^{4}}{\partial y^{4}} \frac{\partial^{2}}{\partial x^{2}} u （ x ， y ） + \frac{\partial^{6}}{\partial y^{6}} u （ x ， y ）$

显示顺序的表达式 $O （ h^{4} ）$ 无法识别为拉普拉斯算子三次幂的某些倍数，因为系数无法匹配。

总结

对于具有距离的正方形网格h在相邻的网格点和上面选择的权重之间，你得到:

${年代}_{h} ＝ Δ u + \frac{h^{2}}{12} Δ^{2} u + O （ h^{4} ）．（ 2 ）$

将此展开式用于特征值问题的数值方法 $Δ u ＝ λ u$ ．加上的倍数 $Δ^{2} u ＝ λ^{2} u$ 到特征值方程。

$Δ u + \frac{h^{2}}{12} Δ^{2} u ＝（ λ + \frac{h^{2}}{12} λ^{2} ） u ．$

这个等式的左边可以很好地用模板来近似 ${年代}_{h}$ ．因此，使用(2)，一个数值特征值 $μ$ 模板的效果令人满意 ${年代}_{h} u ＝ μ u$ 一定是特征值的近似值 $λ$ 的拉普拉斯算子

$μ ＝ λ + \frac{h^{2}}{12} λ^{2} + O （ h^{4} ）．$

对于给定 $μ$ ，解决 $λ$ 为了获得拉普拉斯特征值的更好近似值，请注意，在求解二次方程时 $λ$ 平方根的正确符号由以下要求给出： $λ \to μ$ 为 $h \to 0$ ．

$λ ＝ \frac{6}{h^{2}} （ \sqrt{1 + \frac{μ h^{2}}{3.}} - 1 ）＝ \frac{2 μ}{\sqrt{1 + \frac{μ h^{2}}{3.}} + 1} ．（ 3. ）$

用符号矩阵求解特征值问题

考虑一个l型区域 $Ω$ 由三个单元正方形组成。

$Ω ＝｛（ x ， y ）； - 1 \leq x \leq 0 ， - 1 \leq y \leq 0 ｝ \cup ｛（ x ， y ）； 0 \leq x \leq 1 ， - 1 \leq y \leq 0 ｝ \cup ｛（ x ， y ）； - 1 \leq x \leq 0 ， 0 \leq y \leq 1 ｝．$

定义区域角落的坐标值。

xmin=-1；xmax=1；ymin=-1；ymax=1；

考虑一个由奇数组成的正方形网格Nx = 2 * nx-1网格点的x方向与奇数纽约= 2 * ny-1网格点的y方向。

nx = 6;Nx = 2 * nx-1;hx = (xmax-xmin) / (Nx-1);纽约= 6;纽约= 2 * ny-1;hy = (ymax-ymin) / (Ny-1);

创建一个纽约-借-Nx符号矩阵 $u$ ．它的条目u（i，j）表示值u（xmin+（j-1）*hx，ymin+（i-1）*hy）解决方案u（x，y）关于特征值问题 $Δ u ＝ λ u$ ．

u=sym(“u”(纽约Nx]);

的边界 $Ω$ 对应以下指数:

左边界对应于(i = 1:Ny, j = 1)．
下边界对应于(i = 1, j = 1:Nx)．
右边界对应于(i = 1:ny, j = Nx)和(i = ny: ny, j = nx)．
上界对应于(i = Ny, j =1:nx)和(i = ny, j = nx: nx)．

u (: 1) = 0;%左边界u (1) = 0;%下边界u(1:纽约,Nx) = 0;%右边界，上部u(纽约:纽约,nx) = 0;%右边界，下部u（Ny，1:nx）=0；%左上边界u（ny，nx:nx）=0；%右上边界

该地区与 $0 < x \leq 1$ 和 $0 < y \leq 1$ 不属于 $Ω$ ．设置相应的矩阵项（i=ny+1:ny，j=nx+1:nx）归零。它们不再扮演进一步的角色，将被忽略。

u(ny + 1: ny,nx + 1: nx) = 0;

问题的未知量为以下矩阵项:

u=
                 
                   $（ \begin{array}{c} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & u_{2 ， 2} & u_{2 ， 3.} & u_{2 ， 4} & u_{2 ， 5} & u_{2 ， 6} & u_{2 ， 7} & u_{2 ， 8} & u_{2 ， 9} & u_{2 ， 10} & 0 \\ 0 & u_{3. ， 2} & u_{3. ， 3.} & u_{3. ， 4} & u_{3. ， 5} & u_{3. ， 6} & u_{3. ， 7} & u_{3. ， 8} & u_{3. ， 9} & u_{3. ， 10} & 0 \\ 0 & u_{4 ， 2} & u_{4 ， 3.} & u_{4 ， 4} & u_{4 ， 5} & u_{4 ， 6} & u_{4 ， 7} & u_{4 ， 8} & u_{4 ， 9} & u_{4 ， 10} & 0 \\ 0 & u_{5 ， 2} & u_{5 ， 3.} & u_{5 ， 4} & u_{5 ， 5} & u_{5 ， 6} & u_{5 ， 7} & u_{5 ， 8} & u_{5 ， 9} & u_{5 ， 10} & 0 \\ 0 & u_{6 ， 2} & u_{6 ， 3.} & u_{6 ， 4} & u_{6 ， 5} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & u_{7 ， 2} & u_{7 ， 3.} & u_{7 ， 4} & u_{7 ， 5} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & u_{8 ， 2} & u_{8 ， 3.} & u_{8 ， 4} & u_{8 ， 5} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & u_{9 ， 2} & u_{9 ， 3.} & u_{9 ， 4} & u_{9 ， 5} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & u_{10 ， 2} & u_{10 ， 3.} & u_{10 ， 4} & u_{10 ， 5} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array} ）$

该区域的内部点

Ω ∖ \partial Ω

对应于指数

（ 我 ， j ）

包含未知值

u （ 我 ， j ）

这个问题。将这些未知数收集到一个向量中var．

将符号表达式(由本例第一部分中派生的模板给出)与每个索引(即每个未知数)关联。

因为这个表达式在未知元素中是线性的u（储存于var)，可以将其视为作用于向量的矩阵var．

你可以请客S_h作为拉普拉斯算子的矩阵近似。计算它的特征向量和特征值。

因为对于这种近似，您使用了一个具有少量点的网格，所以只有特征值的前导数字是正确的。

l形区域上拉普拉斯算子的第三大特征值

Ω

是清楚的。拉普拉斯算子的准确本征函数就是这个函数

u （ x ， y ） ＝ 罪 （ π x ） 罪 （ π y ）

与（精确的）特征值有关

- 2 π^{2} ＝ - 19 ． 7 3. 92 ． ． ．

．实际上，使用上面的方程(3)，你可以得到拉普拉斯特征值的一个更好的近似

λ

从模板特征值

μ

：

用双精度矩阵求解特征值问题

当您使用符号矩阵时，不建议大幅增加网格点的数量，因为符号计算比MATLAB双精度矩阵的数值计算要慢得多。这个例子的这一部分演示了如何使用稀疏双精度算法来细化数值网格。l型区域

Ω

设置的方式和以前一样。不要用符号未知数来表示内部点，而是初始化网格值u一个一个地定义

Ω

通过将边界点和外部点的值设为零。不是为每个内部点定义一个符号表达式并计算模具作为雅可比矩阵，而是直接将模具矩阵设置为一个稀疏矩阵。

在这里,S_h为(稀疏)模板矩阵。使用eigs它处理稀疏矩阵来计算三个最大的特征值。

D（3,3）近似精确的特征值

- 2 π^{2} ＝ - 19 ． 7 3. 92088 ． ． ．

．利用上面的方程(3)，导出拉普拉斯特征值的更精确的近似

λ

从模板特征值

μ

．

注意MATLAB膜函数通过不同的方法计算拉普拉斯算子的特征函数。

拉普拉斯算子的特征值

膜问题

九点有限差分近似

包含四阶和高阶导数的项

总结

用符号矩阵求解特征值问题

用双精度矩阵求解特征值问题

符号数学工具箱文档

万博1manbetx

数学建模与符号数学工具箱