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

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

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

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

Самый простой подход к задаче о собственных значениях должен аппроксимировать Лапласиан $$\Delta u$$ приближением конечной разности (шаблон) на квадратной сетке точек с расстояниями hx \in $$x$$ направление и расстояния hy \in $$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) содержа 0th, 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 выбран, чтобы быть порядка 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\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(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)$$ не может быть идентифицирован как некоторое кратное третья степень оператора Лапласа, потому что коэффициенты не могут быть соответствующими.

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

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

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

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

Для данного $$\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- 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)$$

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

Третье по высоте собственное значение оператора Лапласа на L-образной области $$\Omega$$ известен точно. Точная собственная функция оператора Лапласа является функцией $$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);

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

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