В этом примере показано, как решить задачу о собственных значениях оператора Лапласа на L-образной области.
Рассмотрите мембрану, которая фиксируется за пределами из области в плоскости. Его смещение описан задачей о собственных значениях , где оператор Лапласа и скалярный параметр. Граничное условие \forall .
Оператор Лапласа является самопримыкающим и отрицательный определенный, то есть, только действительные отрицательные собственные значения существовать. Существует максимальное (отрицательное) дискретное собственное значение, соответствующая собственная функция называется стандартным состоянием. В этом примере, L-образная область, и стандартное состояние, сопоставленное с этой областью, является L-образной мембраной, которая является логотипом MATLAB®.
Самый простой подход к задаче о собственных значениях должен аппроксимировать Лапласиан приближением конечной разности (шаблон) на квадратной сетке точек с расстояниями hx
\in направление и расстояния hy
\in направление. В этом примере, аппроксимированном с суммой S_h
из девяти регулярных узлов решетки вокруг средней точки . Неизвестные являются весами .
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
относительно и . Начиная с S_h
представляет , коэффициенты всех других производных u
должен быть нуль. Извлеките коэффициенты, заменив все производные u
, кроме и , 0. Замена и 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 =
parameters
parameters =
Используйте subs
функционируйте, чтобы заменить веса их вычисленными значениями.
C = simplify(subs(C));
Выражения C(7)
, C(6)
, и C(4)
содержа 0th, 1-е, и 3-и производные u
исчезнуть.
[C(7), C(6), C(4)]
ans =
Выражение C(5)
Лапласиан u
.
C(5)
ans =
Таким образом, со значениями весов, вычисленных выше, шаблон 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) =
Выберите различные производные 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 =
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 =
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 =
Поэтому можно выбрать 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) =
Теперь рассмотрите члены третьего порядка в hx
, hy
.
simplify(C(2))
ans =
Поскольку никакие такие условия не существуют в расширении шаблона, термина в (1) имеет на самом деле порядок . Рассмотрите члены четвертого порядка шаблона.
factor(simplify(C(1)))
ans =
Проверяйте, могут ли эти условия быть идентифицированы с еще одной степенью оператора Лапласа. Однако сравнение с
Laplace(Laplace(Laplace(u)))
ans(x, y) =
показывает что выражения порядка не может быть идентифицирован как некоторое кратное третья степень оператора Лапласа, потому что коэффициенты не могут быть соответствующими.
Для квадратной сетки с расстоянием h
между соседними узлами решетки и вышеупомянутым выбором весов, вы добираетесь:
Используйте это расширение в числовом подходе к задаче о собственных значениях . Добавьте некоторое кратное к уравнению собственного значения.
Левая сторона этого уравнения хорошо аппроксимирована шаблоном . Таким образом, с помощью (2), числовое собственное значение из удовлетворения шаблона должно быть приближение собственного значения из оператора Лапласа с
Для данного , решите для получить лучшее приближение собственного значения Лапласа. Обратите внимание на то, что в решении квадратного уравнения в правильный знак квадратного корня дан требованием это для .
Рассмотрите L-образную область состоя из трех модульных квадратов.
Задайте координатные значения углов области.
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(i,j)
представляйте значения u(xmin + (j - 1)*hx,ymin + (i - 1)*hy)
из решения u(x,y)
из задачи о собственных значениях .
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
Область с и не принадлежит . Установите соответствующие матричные записи (i = ny + 1:Ny, j = nx + 1:Nx)
обнулять. Они не играют дальнейшей роли и будут проигнорированы.
u(ny + 1:Ny,nx + 1:Nx) = 0;
Неизвестные проблемы являются следующими матричными записями:
u
u =
Внутренние точки области соответствуйте индексам это содержит неизвестные значения из проблемы. Соберите эти неизвестные в векторном 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 =
С тех пор для этого приближения вы использовали сетку с маленьким числом точек, только ведущие цифры собственных значений правильны.
Третье по высоте собственное значение оператора Лапласа на L-образной области известен точно. Точная собственная функция оператора Лапласа является функцией сопоставленный с (точным) собственным значением . Действительно, с помощью уравнения (3) выше, можно вывести лучшее приближение собственного значения Лапласа от собственного значения шаблона :
mu = D(3,3)
mu =
lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)
lambda =
Постройте собственную функцию, сопоставленную с третьим по высоте собственным значением.
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)
аппроксимирует точное собственное значение . Используя уравнение (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);