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

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

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

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

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

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

Самый простой подход к задаче о собственных значениях должен аппроксимировать Лапласиан Δu приближением конечной разности (шаблон) на квадратной сетке точек с расстояниями hx \in x направление и расстояния hy \in 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 = 

hx212=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 = 

hy212=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 = 

(1360hhhh6x6 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);

Создайте 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<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.5214641572625960021345709535953-14.431096242107969492574666743957-18.490392088545609858994660377955)

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

Третье по высоте собственное значение оператора Лапласа на L-образной области Ω известен точно. Точная собственная функция оператора Лапласа является функцией u(x,y)=sin(πx)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 object. The axes object 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 object. The axes object contains an object of type surface.

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

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

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

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