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

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

Мембранная проблема

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

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

Приближение конечной разности с девятью точками

Самый простой подход к задаче о собственных значениях должен аппроксимировать Лапласиан Δu приближением конечной разности (шаблон) на квадратной сетке точек с расстояниями hx в x направление и расстояния hy в y направление. В этом примере, аппроксимированном Δu с суммой S_h девяти регулярных узлов решетки вокруг средней точки (x,y). Неизвестные являются весами a-1-1,,a11.

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 = 

(z1hy2-2zz1hx2-2z4z-2hx2-2hy21hy2-2zz1hy2-2zz)

parameters
parameters = z

Используйте функцию subs, чтобы заменить веса их вычисленными значениями.

C = simplify(subs(C));

Выражения C(7), C(6) и C(4), содержащий 0th, 1-е, и 3-и производные u исчезают.

[C(7), C(6), C(4)]
ans = (000)

Выражение C(5) является Лапласианом u.

C(5)
ans = 

2x2 u(x,y)+2y2 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) = 

d4x4 u(x,y)+2d2y2 2x2 u(x,y)+d4y4 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 = 0.0833hx2=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 = hx2hy2z=2d
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 = 0.0833hy2=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=Δ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 = 

(0.0028hhhh6x6 u(x,y)+52y2 4x4 u(x,y)+54y4 2x2 u(x,y)+6y6 u(x,y))

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

Laplace(Laplace(Laplace(u)))
ans(x, y) = 

6x6 u(x,y)+32y2 4x4 u(x,y)+34y4 2x2 u(x,y)+6y6 u(x,y)

показывает что выражения порядка O(h4) не может быть идентифицирован как некоторое кратное третья степень оператора Лапласа, потому что коэффициенты не могут быть соответствующими.

Сводные данные

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

Sh=Δu+h212Δ2u+O(h4).(2)

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

Δu+h212Δ2u=(λ+h212λ2)u.

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

μ=λ+h212λ2+O(h4).

Для данного μ, решите для λ получить лучшее приближение собственного значения Лапласа. Обратите внимание на то, что в решении квадратного уравнения в λ правильный знак квадратного корня дан требованием это λμ для h0.

λ=6h2(1+μh23-1)=2μ1+μh23+1.(3)

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

Рассмотрите L-образную область Ω состоя из трех модульных квадратов.

Ω={(x,y);-1x0,-1y0}{(x,y);0x1,-1y0}{(x,y);-1x0,0y1}.

Задайте координатные значения углов области.

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);

Создайте Nx Ny символьной матрицей 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<x1 и 0<y1 не принадлежит Ω. Установите соответствующие матричные записи (i = ny + 1:Ny, j = nx + 1:Nx) обнулять. Они не играют дальнейшей роли и будут проигнорированы.

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

Неизвестные проблемы являются следующими матричными записями:

u
u = 

(000000000000u2,2u2,3u2,4u2,5u2,6u2,7u2,8u2,9u2,1000u3,2u3,3u3,4u3,5u3,6u3,7u3,8u3,9u3,1000u4,2u4,3u4,4u4,5u4,6u4,7u4,8u4,9u4,1000u5,2u5,3u5,4u5,5u5,6u5,7u5,8u5,9u5,1000u6,2u6,3u6,4u6,50000000u7,2u7,3u7,4u7,50000000u8,2u8,3u8,4u8,50000000u9,2u9,3u9,4u9,50000000u10,2u10,3u10,4u10,500000000000000000)

Внутренние точки области ΩΩ соответствуйте индексам (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 = (-9.5215-14.4311-18.4904)

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

Третье по высоте собственное значение оператора Лапласа на L-образной области Ω известен точно. Точная собственная функция оператора Лапласа является функцией u(x,y)=sin(πx)sin(πy) сопоставленный с (точным) собственным значением -2π2=-19.7392.... Действительно, с помощью уравнения (3) выше, можно вывести лучшее приближение собственного значения Лапласа λ от собственного значения шаблона μ:

mu = D(3,3)
mu = -18.4904
lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)
lambda = -19.7968

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

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);

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

Когда вы используете символьные матрицы, увеличивание числа узлов решетки решительно не рекомендуется, потому что символьные вычисления значительно медленнее, чем численные расчеты с 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);

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

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

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