Собственные значения оператора Лапласа

В этом примере показано, как решить проблему собственных значений оператора Лапласа в L-образной области.

Проблема с мембраной

Рассмотрим мембрану, которая зафиксирована на граничной $∂Ω$ области $Λ$ в плоскости. Его смещение $u (x,$ y) описывается $задачей$ собственного значения Δu = λ u, $где_{Δu} =_{}$ uxx + uyy - оператор Лапласа, а λ - скалярный параметр. Граничное $условие u$ (x, y) $= 0 для$ всех (x, y) ∈∂Ω.

Оператор Лапласа является самосопряжённым и отрицательным определённым, то есть существуют только действительные отрицательные собственные значения λ. Существует максимальное (отрицательное) дискретное собственное значение, соответствующая собственная функция $u$ называется состоянием основания. В этом примере Λ представляет собой L-образную область, а основное состояние, связанное с этой областью, представляет собой L-образную мембрану, которая является логотипом MATLAB ®.

Аппроксимация конечных разностей в девяти точках

Простейший подход к задаче собственного значения заключается в аппроксимации лапласиана $Δu$ конечной разностной аппроксимацией (трафаретом) на квадратной сетке точек с расстояниями hx в направлении x и расстояниях hy в направлении y. В этом примере аппроксимировать $Δu$ суммой S_h девяти правильных точек сетки вокруг средней точки $(x, y$ ). Неизвестные - веса $_{} а-1-1, . ._{.,}$ а11.

syms u(x,y) Eps a11 a10 a1_1 a01 a00 a0_1 a_11 a_10 a_1_1
syms hx hy positive
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) +... 
      a10  * u(x + Eps*hx,y) +...
      a_1_1* u(x - Eps*hx,y - Eps*hy) +...
      a0_1 * u(x,y - Eps*hy) +...
      a1_1 * u(x + Eps*hx,y - Eps*hy);

Использовать символический параметр Eps сортировать расширение этого выражения по полномочиям hx и hy. Зная веса, можно приблизить лапласиан, установив Eps = 1.

t = taylor(S_h, Eps, 'Order', 7);

Используйте coeffs функция извлечения их коэффициентов членов с одинаковыми степенями Eps. Каждый коэффициент является выражением, содержащим степени hx, hyи производные от u относительно $x$ и $y$ . С тех пор S_h представляет $_{собой uxx} +_{}$ uyy, коэффициенты всех других производных u должно быть равно нулю. Извлеките коэффициенты путем замены всех производных u, кроме $_{uxx}$ и $_{uyy}$ , на 0. Замените $_{uxx}$ и $_{uyy}$ на 1. Это уменьшает расширение Тейлора до коэффициента, который требуется вычислить, и приводит к следующим шести линейным уравнениям.

C = formula(coeffs(t, Eps, 'All'));
eq0  = subs(C(7),u(x,y),1) == 0;
eq11 = subs(C(6),[diff(u,x),diff(u,y)],[1,0]) == 0;
eq12 = subs(C(6),[diff(u,x),diff(u,y)],[0,1]) == 0;
eq21 = subs(C(5),[diff(u,x,x),diff(u,x,y),diff(u,y,y)],[1,0,0]) == 1;
eq22 = subs(C(5),[diff(u,x,x),diff(u,x,y),diff(u,y,y)],[0,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 = subs(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 = subs(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 = subs(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 = subs(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;

Решите результирующие десять уравнений для девяти неизвестных весов. Использовать ReturnConditions для поиска всех решений, включая произвольные параметры.

[a11,a10,a1_1,a01,a00,a0_1,a_11,a_10,a_1_1,parameters,conditions] = ...
    solve([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], ...
          'ReturnConditions',true);

expand([a_11,a01,a11;...
        a_10,a00,a01;...
        a1_1,a0_1,a_1_1])

ans = 
 $(\begin{array}{c} z & \frac{1}{{hy}^{2}} - 2 z & z \\ \frac{1}{{hx}^{2}} - 2 z & 4 z - \frac{2}{{hx}^{2}} - \frac{2}{{hy}^{2}} & \frac{1}{{hy}^{2}} - 2 z \\ z & \frac{1}{{hy}^{2}} - 2 z & z \end{array})$

parameters

parameters =  $z$

Используйте subs для замены весов их вычисленными значениями.

C = simplify(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. Другой вариант - превратить его в кратное квадрату оператора Лапласа.

syms d
Laplace = @(u) laplacian(u,[x,y]);
expand(d*Laplace(Laplace(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)и приравнять их коэффициенты к соответствующим терминам.

subs(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$

subs(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} {hy}^{2} z = 2 d$

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

ans = 
 $\frac{{hy}^{2}}{12} = d$

Поэтому можно выбрать d = hx^2/12 = hy^2/12 и z = 2*d/(hx^2*hy^2), подразумевая, что hx = hy и z = 1/(6*hx*hy). Следовательно, трафарета S_h аппроксимирует модифицированный лапласиан на квадратной сетке с hx = hy = h.

$_{Sh} = \frac{{Δu}^{}}{+}^{} {h212Δ2u}^{} + O ($ h3). (1)

syms h
hx = h;
hy = h;
d = h^2/12;

Заменить hx и hy около h.

C = subs(C);

Заменить z по своей стоимости, 1/(6*h^2). Поскольку z не существует в рабочей области MATLAB ®, доступ к ней возможен только в виде значения, сохраненного в parameters массив.

C = subs(C,parameters,1/(6*h^2));

Проверьте формулу (1).

simplify(C(3) - d*Laplace(Laplace(u)))

ans(x, y) =  $0$

Теперь рассмотрим термины третьего порядка в hx, hy.

simplify(C(2))

ans =  $0$

Поскольку таких терминов в расширении шаблона не существует, термин $O (^{} h3$ ) в (1) фактически имеет порядок $O^{(}$ h4). Рассмотрим термины четвертого порядка шаблона.

factor(simplify(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})$

Проверьте, можно ли идентифицировать эти термины с помощью еще одной мощности оператора Лапласа. Однако сравнение с

Laplace(Laplace(Laplace(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 (^{} h4$ ) не могут быть идентифицированы как несколько кратные третьей степени оператора Лапласа, поскольку коэффициенты не могут быть сопоставлены .

Резюме

Для квадратной сетки с расстоянием h между соседними точками сетки и приведенными выше вариантами весов получается:

$_{Sh} = \frac{{Δu}^{}}{+}^{} {h212Δ2u}^{} + O ($ h4). (2)

Используйте это расширение для числового подхода к задаче собственного значения $Δu =$ λ u. Добавьте к $^{}^{}$ уравнению собственного значения кратное Δ2u = λ 2u.

$Δu \frac{+^{}}{}^{} h212Δ2u = \frac{^{(}}{λ} +^{}$ h212λ 2) u.

Левая часть этого уравнения хорошо аппроксимируется шаблоном $_{Sh}$ . Таким образом, используя (2), числовое собственное значение $λ$ шаблона, удовлетворяющего $_{Shu} =$ micu, должно быть аппроксимацией собственного значения $λ$ оператора Лапласа с

$λ = \frac{λ^{}}{+}^{} {h212λ 2}^{} +$ O (h4).

Для данного $λ$ решите для $λ$ , чтобы получить лучшее приближение собственного значения Лапласа. Заметим, что в решении квадратичного уравнения в $λ$ правильный знак квадратного корня задаётся требованием, которое $λ\toμ$ для $h\to0 .$

$λ \frac{=}{^{}} \sqrt{6h2 \frac{(1^{}}{+}} \frac{}{\sqrt{мкх23-1) \frac{=^{}}{}} 2 мк1} +$ мкх23 + 1. (3)

Использование символьных матриц для решения проблемы собственных значений

Рассмотрим L-образную область $Λ$ , состоящую из трёх единичных квадратов.

$Λ = {(x, y); - 1≤x≤0, - 1\leqy\leq0}\cup {(x, y); 0\leqx\leq1, - 1\leqy\leq0}\cup {(x, y);$ - 1≤x≤0,0≤y≤1}.

Определите значения координат углов области.

xmin =-1;
xmax = 1; 
ymin =-1; 
ymax = 1;

Рассмотрим квадратную сетку, состоящую из нечетного числа Nx=2*nx-1 точек сетки в x направление и нечетное число Ny=2*ny-1 точек сетки в y направление.

nx = 6; 
Nx = 2*nx-1; 
hx = (xmax-xmin)/(Nx-1);

ny = 6; 
Ny = 2*ny-1; 
hy = (ymax-ymin)/(Ny-1);

Создание Nyоколо- Nx символьная матрица $u$ . Его записи u(i,j) представляют значения u(xmin + (j - 1)*hx,ymin + (i - 1)*hy) решения u(x,y) задачи собственного значения $Δu =$ λ u.

u = sym('u',[Ny,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;      % left boundary
u(1,:) = 0;      % lower boundary
u(1:ny,Nx) = 0;  % right boundary, upper part
u(ny:Ny,nx) = 0; % right boundary, lower part
u(Ny,1:nx) = 0;  % upper boundary, left part
u(ny,nx:Nx) = 0; % upper boundary, right part

Область с $0<x≤1$ и $0<y≤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})$

Внутренние точки области $Ω∖∂Ω$ соответствовать индексам $(i, j$ ), которые содержат неизвестные значения $u (i,$ j) проблемы. Собрать неизвестных в векторvars.

[I,J] = find(u~=0);
vars = u(u~=0);

Ассоциируйте символическое выражение (данное шаблоном, полученным в первой части этого примера) с каждым индексом (то есть с каждым неизвестным).

n = length(vars);
Lu = sym(zeros(n,1));
for k=1:n
   i = I(k);
   j = J(k);
   Lu(k) =   1/6*u(i+1,j-1) +  2/3*u(i+1,j) + 1/6*u(i+1,j+1) ...
           + 2/3*u( i ,j-1) - 10/3*u( i ,j) + 2/3*u( i ,j+1) ...
           + 1/6*u(i-1,j-1) +  2/3*u(i-1,j) + 1/6*u(i-1,j+1); 
end
Lu = Lu/hx^2;

Поскольку это выражение является линейным в неизвестных элементах u (хранится в vars), вы можете рассматривать его как матрицу, действующую на вектор vars.

S_h = jacobian(Lu, vars);

Вы можете лечить S_h как матричное приближение оператора Лапласа. Вычислите собственные векторы и собственные значения.

[V,D] = eig(vpa(S_h));

Три максимальных собственных значения задаются первыми диагональными элементами D.

[D(1,1), D(2,2), D(3,3)]

ans =  $(\begin{array}{c} - 9.5214641572625960021345709535953 & - 14.431096242107969492574666743957 & - 18.490392088545609858994660377955 \end{array})$

Поскольку для этого приближения использовалась сетка с небольшим количеством точек, верны только первые цифры собственных значений.

Точно известно третье по величине собственное значение оператора Лапласа в L-образной области Λ. Точная собственная функция оператора Лапласа - это функция $u (x, y) = sin (securityx)$ sin (δy), связанная с (точным) $^{} собственным значением - 2ā2 = -$ 19.7392.... Действительно, используя приведенное выше уравнение (3), можно получить лучшее приближение собственного значения Лапласа λ от собственного значения по трафарету λ:

mu = D(3,3)

mu =  $- 18.490392088545609858994660377955$

lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)

lambda =  $- 19.796765119155672176257649532142$

Постройте график собственной функции, связанной с третьим по величине собственным значением.

v = V(:,3);
for k=1:n
  u(I(k),J(k)) = v(k);
end
u = double(u);
surf(xmin:hx:xmax,ymin:hy:ymax,u');
view(125, 30);

Figure contains an axes. The axes contains an object of type surface.

Использование матриц с двойной точностью для решения проблемы собственных значений

При использовании символьных матриц резко увеличивать число точек сетки не рекомендуется, поскольку символьные вычисления значительно медленнее числовых вычислений с матрицами двойной точности MATLAB. В этой части примера показано, как использовать разреженную двойную арифметику, которая позволяет уточнить числовую сетку. L-образная область $Λ$ устанавливается таким же образом, как и ранее. Вместо обозначения внутренних точек символическими неизвестными, инициализируйте значения сетки u по единицам и определить, задав нулевые значения граничных и внешних точек. Вместо определения символьного выражения для каждой внутренней точки и вычисления шаблона в виде матрицы якобиана настройте шаблон непосредственно в виде разреженной матрицы.

xmin =-1; 
xmax = 1; 
ymin =-1; 
ymax = 1;

nx = 30; 
Nx = 2*nx-1; 
hx = (xmax-xmin)/(Nx-1);

ny = 30; 
Ny = 2*ny-1; 
hy = (ymax-ymin)/(Ny-1);

u = ones(Ny,Nx); 
u(:,1) = 0;      % left boundary
u(1:ny,Nx) = 0;  % right boundary, upper part
u(ny:Ny,nx) = 0; % right boundary, lower part
u(1,:) = 0;      % lower boundary
u(Ny,1:nx) = 0;  % upper boundary, left part
u(ny,nx:Nx) = 0; % upper boundary, right part
u(ny + 1:Ny,nx + 1:Nx) = 0;

[I,J] = find(u ~= 0);
n = length(I);
S_h = sparse(n,n);
for k=1:n 
  i = I(k);
  j = J(k);
  S_h(k,I==i+1 & J==j+1)=  1/6;
  S_h(k,I==i+1 & J== j )=  2/3;
  S_h(k,I==i+1 & J==j-1)=  1/6;
  S_h(k,I== i  & J==j+1)=  2/3;
  S_h(k,I== i  & J== j )=-10/3;
  S_h(k,I== i  & J==j-1)=  2/3;
  S_h(k,I==i-1 & J==j+1)=  1/6;
  S_h(k,I==i-1 & J== j )=  2/3;
  S_h(k,I==i-1 & J==j-1)=  1/6;
end
S_h = S_h./hx^2;

Здесь, S_h - (разреженная) матрица шаблона. Использовать eigs который обрабатывает разреженные матрицы для вычисления трех самых больших собственных значений.

[V,D] = eigs(S_h,3,'la');

Три максимальных собственных значения являются первыми диагональными элементами D.

[D(1,1), D(2,2), D(3,3)]

ans = 1×3

   -9.6493  -15.1742  -19.7006

D(3,3) аппроксимирует точное собственное значение $-^{2ā2} = - 19.7392088$ .... Используя приведенное выше уравнение (3), выведите более точную аппроксимацию $собственного значения$ Лапласа λ из $собственного значения$ λ шаблона.

mu = D(3,3)

mu = -19.7006

lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)

lambda = -19.7393

Постройте график собственной функции, связанной с третьим по величине собственным значением.

v = V(:,3);
for k=1:n
  u(I(k),J(k)) = v(k);
end
surf(xmin:hx:xmax, ymin:hy:ymax,u');
view(125, 30);

Figure contains an axes. The axes contains an object of type surface.

Обратите внимание, что MATLAB membrane функция вычисляет собственные функции оператора Лапласа различными методами.

membrane(3, nx - 1, 8, 8);

Figure contains an axes. The axes contains an object of type surface.

Документация

Собственные значения оператора Лапласа

Проблема с мембраной

Аппроксимация конечных разностей в девяти точках

Термины, содержащие производные четвертого и высшего порядка

Резюме

Использование символьных матриц для решения проблемы собственных значений

Использование матриц с двойной точностью для решения проблемы собственных значений

Документация по инструментам символьной математики

Поддержка