回顾50年前的状态之后,我使用经典的有限差分方法,其次是推断,找到的第一特征值区域潜在MathWorks标志。

内容

五十年前

我的博士论文,在1965年提交给斯坦福大学数学系题为“有限差分方法对拉普拉斯算子的特征值”。论文不仅对L型地区,但L是我主要的例子。

我的论文导师是乔治活力四射,L在他自己的一些研究调查。当时,许多人感兴趣的方法获得保证微分算子特征值的上下界。活力四射表明,某些类型的有限差分方法可以为凸区域提供下界。但L不是凸,所以一个地区的调查。

没人知道第一特征值非常精确。活力四射和沃尔夫冈Wasow在1960年出版了一本书,“偏微分方程的有限差分方法”(约翰·威利& Sons),这个特性L在几个不同的部分。24.3节中活力四射的引用一些计算1954年他做了,他估计

$ $ \ lambda_1 = 9.636 $ $

他继续报告一个非正式的猜想在50年代末从他的同事r . De Vogelaere

< \ lambda_1 < 9.63968美元9.63972美元美元

我在这里重复这些数字给你一些想法的准确性我们处理。我们将在另一篇博文中,正确的答案是在De Vogelaere上界。

1964年我的计算

这是我的论文的一个页面。页面的底部的等高线图最终演变成MathWorks标志。情节了数控绘图机绘图机涉及计算机控制电磁拿着圆珠笔或液体墨水笔在纸卷上移动一个旋转的鼓。这是一个非常有效的使用设备和计算机中心全球许多年了。

计算显示在这个页面是第一个特征值的有限差分网格大小的拉普拉斯算子十值h美元。最好的网格,h = 1/100美元,有29601个室内点。我记得,第一特征值的计算,我用一种逆幂法,花了半个小时在斯坦福的校园主框架,一个IBM 7090。

这是一个阴谋我计算五十年前的这些值。我添加了一个黑线在我们现在知道的极限值,第一个连续微分算子的特征值。你可以看到,有限差分值收敛非常缓慢。此外,他们从上面收敛。活力四射是不会得到这个非凸区域的下界。

thesis_graph

收敛速度慢

有限差分特征值的显示收敛速度和减少网格大小是由相应的本征函数是凹角的奇异。在那个转角梯度变得无限。你可以看到轮廓线附近的角落里挤在一起。如果你想让一个鼓头或小手鼓拉伸膜/ l型框架,将撕裂。

更准确地说,使用极坐标,r和θ\美元,美元为中心凹角。注意,覆盖区域的内部\θ将美元从$ 0 $美元\压裂{3}{2}\π美元。事实证明,接近角第一本征函数的渐近行为

$ $ u (r,θ\)\大约r ^ \α\{\α、θ}$ $的罪

美元\α= \压裂{2}{3}$。这意味着k th美元衍生品像的本征函数的奇点

$ $ \部分^ k u \大约r ^{\α_k} $ $

当这个事实结合的泰勒级数分析五点有限差分算子的离散化误差我们发现第一特征值收敛速度

$ $ \ lambda_ {1 h} - \ lambda_1 \大约O (h ^{2 \α})= O (h ^ \压裂{4}{3})$ $

这是明显慢于O (h ^ 2)美元时收敛速度得到了本征函数并不奇异。

在MATLAB差分方法

现在,让我们在MATLAB中使用我的笔记本电脑和稀疏功能。的函数numgrid,delsq,间谍,eigs已经包含在吗sparfun目录以来介绍。我们首先生成一个小l型编号方案。

n = 10;L = numgrid (“L”n + 1)

L = 0 0 0 0 0 0 0 0 0 0 0 0 1 5 9 13 17 21 30 39 48 14 0 0 2 6 10 18 22日31日40 49 0 0 3 7 11 15 19 23 32 41 50 0 0 4 8 12 16 20 24 33 51 42 0 0 0 0 0 0 0 25 34 43 52 0 0 0 0 0 0 0 26 35 44 53 0 0 0 0 0 0 0 27 36 45 54 0 0 0 0 0 0 0 28 37 46 55 0 0 0 0 0 0 0 29 38 47 56 0 0 0 0 0 0 0 0 0 0 0 0

为了说明这个方案如何运作,重叠网格有限差分。

Lgrid(左)

该声明

一个= delsq(左);

五点离散拉普拉斯算子的构造矩阵表示的网格l。

谁一个issym = isequal (,) nnzA = nnz (A)

类属性名称大小字节56 x56 3156双稀疏issym = 1 nnzA = 244

所以对于这个网格,离散拉普拉斯算子一个是一个56-by-56稀疏双精度与244年对称矩阵非零元素。这是一个平均的

率= nnz (A) /尺寸(A, 1)

率= 4.3571

非零元素每一行。每个内部网格点连接到它的四个最近的邻居。例如,44号连接指向点35,43岁,45岁和53岁。所以非零元素列44

(:44)

ans =(35岁,1)1(43岁,1)1(44岁,1)4(45岁,1)1(53岁,1)1

对角线上有4和1的位置相对应的非对角的邻居之间的连接。他们可以看到间谍情节,展示了非零的位置。

间谍(A)

规模一个通过筛孔尺寸的平方,然后使用eigs找到五个最小的特征值。

h = 2 / n = A / h ^ 2;eigvals = eigs (5“sm”)

h = 0.2000 eigvals = 30.3140 27.7177 19.0983 14.6708 9.6786

2014年我的计算

竞选减少网格大小,我的笔记本电脑的内存容量。

类型L_finite_diffs

% L_finite_diffs %剧本MathWorks标志,第二部分%有限选择的差异。托尔= 1.0 e-12;eigvals = 0(13日5);抽搐n = 50:50:650 L = numgrid (L, n + 1);h = 2 / n;A = delsq(左)/ h ^ 2;e = eigs (5“sm”,选择);:eigvals (n / 50) = fliplr (e”);流(4.0 ' % % 10.6 f % 9.0 f % 10.0 f % 16.12 f \ n ',…n、h、大小(A, 1), nnz (A), e(5))结束toc拯救eigvals eigvals

类型L_finite_diffs_results

大小n h (A) nnz lambda_1, h 1776 0.040000 8684 9.661331286788 100 0.020000 7301 36109 9.649547111146 150 0.013333 16576 9.643932372382 250 0.008000 82284 9.645736333124 200 0.010000 29601 147209 46376 230884 9.642903412412 300 0.006667 66901 333309 9.642247498956 350 0.005714 91176 454484 9.641797238078 400 0.005000 119201 594409 9.641471331746 450 0.004444 150976 753084 9.641225852818 500 0.004000 186501 9.640883203502 930509 9.641035124947 550 0.003636 225776 1126684 600 268801 1341609 9.640759699387 650 0.003077 315576 0.003333 1575284 9.640657572983运行时间是29.729023秒。

外推法

我知道的泰勒级数分析离散误差包括与两美元h ^{4/3} $和$ h ^ 2美元,我们使用这些作为推断的基础。使用反斜杠来做一个最小二乘匹配。

负载eigvalseigvals = eigvals (: 1);n = (50:50:650)”;h = 2. / n;e = 1(大小(h));X = [e h。^ (4/3) h。^ 2];c = X(4:最后,:)\ eigvals(4:结束);λ₁= c (1);n = (50:10:650)”;h = 2. / n;适合c = c (1) + (2) * h。^ (4/3) + c (3) * h。^ 2; plot(n(1:5:end),eigvals,“o”n,健康,“- - -”700年,[0],[λ₁λ₁),“k -”甘氨胆酸)组(,“xtick”100:100:600,“xticklabel”,…num2str (sprintf (“2 / % 3.0 f \ n”,100:100:600)))包含(“h”)

坦率地说,我很惊讶的推断这个结果产生的有限差分特征值。

格式长λ₁

λ₁= 9.639723753239963

我知道从技术,我在后面的文章将描述这个同意连续微分算子的特征值7位小数。之前我写这个博客我从来没有尝试过这个计算在计算机上有足够的内存来处理网格细,也从来没有能够获得这种外推有限差分特征值的准确性。

本征函数

让我们做一个本征函数的等高线图。棘手的部分是使用网格编号l插入的特征向量回电网。

近n = 100;h = 2 / n;L = numgrid (“L”n + 1);A = delsq(左)/ h ^ 2;(V, E) = eigs(一个5“sm”);L (L > 0) = - v (:, 5);contourf (fliplr (L)、12);轴广场轴从

情节展示parula有柄,MATLAB默认colormap图形II和释放2014 b。的colormap从而得名parula,这是一个小莺在北美东部发现展示这些颜色。作为我们的等高线图表明,颜色强度变化平滑和均匀的轮廓水平增加。这是不正确的飞机,旧的违约。史蒂夫·埃丁有更多的关于parula在他的博客。

得到了MATLAB代码

发表与MATLAB®R2014b