Мембранная задача

Рассмотрим мембрану, которая зафиксирована на контуре $$\partial \Omega$$ области $$\Omega$$ в плоскости. Его перемещение $$u(x,y)$$ описывается задачей собственного значения $$\Delta u=\lambda u$$, где $$\Delta u=u_{xx}+u_{yy}$$ является оператором Laplace и $$\lambda$$ является скалярным параметром. Граничное условие $$u(x,y)=0$$ для всех $$(x,y)\in \partial \Omega$$.

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

Девятиточечное приближение конечного различия

Самый простой подход к задаче собственных значений состоит в том, чтобы аппроксимировать лапласиан $$\Delta u$$ конечным различием приближения (трафарет) на квадратной сетке точек с расстояниями hx в $$x$$ направление и расстояния hy в $$y$$ направление. В этом примере аппроксимируйте $$\Delta u$$ с суммой S_h из девяти правильных точек сетки вокруг средней точки $$(x,y)$$. Неизвестные - это веса $$a_{-1-1},\dots ,a_{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 представляет $$u_{xx}+u_{yy}$$, коэффициенты всех других производных u должно быть нулем. Извлечь коэффициенты путем замены всех производных u, кроме $$u_{xx}$$ и $$u_{yy}$$, по 0. Заменить $$u_{xx}$$ и $$u_{yy}$$ по 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 = 
$$\left(\begin{array}{ccc}z& \frac{1}{{\mathrm{hy}}^2}-2\hspace{0.17em}z& z\\ \frac{1}{{\mathrm{hx}}^2}-2\hspace{0.17em}z& 4\hspace{0.17em}z-\frac{2}{{\mathrm{hx}}^2}-\frac{2}{{\mathrm{hy}}^2}& \frac{1}{{\mathrm{hy}}^2}-2\hspace{0.17em}z\\ z& \frac{1}{{\mathrm{hy}}^2}-2\hspace{0.17em}z& z\end{array}\right)$$

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 = $$\left(\begin{array}{ccc}0& 0& 0\end{array}\right)$$

Выражение C(5) - лапласиан u.

C(5)

ans = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)$$

Таким образом, с значениями весов, вычисленными выше, трафарет S_h аппроксимирует лапласиан до порядка hx^2, hy^2 для любых значений произвольного параметра z, при условии, что z выбран таким образом, чтобы иметь порядок <reservedrangesplaceholder0>.

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

Хотя решение содержит свободный параметр z, выражение C(3) содержащие производные четвертого порядка u невозможно превратить в нуль подходящим выбором z. Другая опция - превратить его в произведение квадрата оператора Laplace.

syms d
Laplace = @(u) laplacian(u,[x,y]);
expand(d*Laplace(Laplace(u)))

ans(x, y) = 
$$d\hspace{0.17em}\frac{{\partial }^4}{\partial x^4}\mathrm{ }u\left(x,y\right)+2\hspace{0.17em}d\hspace{0.17em}\frac{{\partial }^2}{\partial y^2}\mathrm{ }\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+d\hspace{0.17em}\frac{{\partial }^4}{\partial y^4}\mathrm{ }u\left(x,y\right)$$

Выберите различные производные 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{{\mathrm{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 = $${\mathrm{hx}}^2\hspace{0.17em}{\mathrm{hy}}^2\hspace{0.17em}z=2\hspace{0.17em}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{{\mathrm{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.

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(h^3)$$ в (1) находится на самом деле в порядке $$O(h^4)$$. Рассмотрим условия четвертого порядка трафарета.

factor(simplify(C(1)))

ans = 
$$\left(\begin{array}{cccccc}\frac{1}{360}& h& h& h& h& \frac{{\partial }^6}{\partial x^6}\mathrm{ }u\left(x,y\right)+5\hspace{0.17em}\frac{{\partial }^2}{\partial y^2}\mathrm{ }\frac{{\partial }^4}{\partial x^4}\mathrm{ }u\left(x,y\right)+5\hspace{0.17em}\frac{{\partial }^4}{\partial y^4}\mathrm{ }\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^6}{\partial y^6}\mathrm{ }u\left(x,y\right)\end{array}\right)$$

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

Laplace(Laplace(Laplace(u)))

ans(x, y) = 
$$\frac{{\partial }^6}{\partial x^6}\mathrm{ }u\left(x,y\right)+3\hspace{0.17em}\frac{{\partial }^2}{\partial y^2}\mathrm{ }\frac{{\partial }^4}{\partial x^4}\mathrm{ }u\left(x,y\right)+3\hspace{0.17em}\frac{{\partial }^4}{\partial y^4}\mathrm{ }\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^6}{\partial y^6}\mathrm{ }u\left(x,y\right)$$

показывает, что выражения порядка $$O(h^4)$$ не может быть идентифицировано как некоторое число, кратное третьей степени оператора Laplace, поскольку коэффициенты не могут совпадать.

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

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

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

Левая сторона этого уравнения хорошо аппроксимируется трафаретом $$S_h$$. Таким образом, использование (2), числового собственного значения $$\mu$$ трафарета удовлетворяющего $$S_{h}u=\mu u$$ должно быть приближением собственного значения $$\lambda$$ оператора Laplace с

Для данного $$\mu$$, решить для $$\lambda$$ получить лучшее приближение собственного значения Лапласа. Обратите внимание, что в решении квадратичного уравнения в $$\lambda$$ правильный знак квадратного корня задается требованием, что $$\lambda \to \mu$$ для $$h\to 0$$.

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

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

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

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-by- Nx символьная матрица $$u$$. Его записи u(i,j) представлять значения u(xmin + (j - 1)*hx,ymin + (i - 1)*hy) решения u(x,y) задачи собственного значения $$\Delta u=\lambda u$$.

u = sym('u',[Ny,Nx]);

Область контуров $$\Omega$$ соответствуют следующим индексам:

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

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

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

u

u = 
$$\left(\begin{array}{ccccccccccc}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}\right)$$

Внутренние точки области $$\Omega \setminus \partial \Omega$$ соответствуют индексам $$(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 = $$\left(\begin{array}{ccc}-9.5214641572625960021345709535953& -14.431096242107969492574666743957& -18.490392088545609858994660377955\end{array}\right)$$

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

Третье по величине собственное значение оператора Laplace на L-образной области $$\Omega$$ Точно известно. Точная собственная функция оператора Laplace является функцией $$u(x,y)=\sin (\pi x)\sin (\pi y)$$ сопоставлен с (точным) собственным значением $$-2{\pi }^2=-19.7392...$$. Действительно, используя уравнение (3) выше, можно вывести лучшее приближение собственного значения Лапласа $$\lambda$$ из собственного значения трафарета $$\mu$$:

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

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

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

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{\pi }^2=-19.7392088...$$. Используя уравнение (3) выше, выведите более точное приближение собственного значения Лапласа $$\lambda$$ из собственного значения трафарета $$\mu$$.

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

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

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