Gilbert Strang , LA for engineers, part 1
recitation 01
1 : four Special Matrices
2 : differential Eqns and difference Eqns
3 : solving a Linear System
4 : Delta Function
recitation 02
5 : Eigenvalues (Part 1)
6 : EigenValues (part 2)
7 : Positive Definite
recitation 03
8 : Springs and Masses
9 : Oscillation
recitation 04
10 : Finite Differences in Time
11 : Least Squares
12 : Graphs and Networks
13 : Kirchhoff's Current Law
recitation 05
recitation 06
14 : exam review
15 : trusses and A⁺*C*A
16 : trusses (part 2)
recitation 07
17 : Finite Elements in 1D (part 1)
18 : Finite Elements in 1D (part 2)
19 : Quadratic/Cubic Elements
20 : Element Matrices; 4th Order Bending Equations
recitation 08
21 : Boundary Conditions, Splines, Gradient, Divergence
22 : Gradient and Divergence
23 : Laplace's Equation
recitation 09
24 : Laplace's Equation (part 2)
25 : Fast Poisson Solver (part 1)
26 : Fast Poisson Solver (part 2); Finite Elements in 2D
01 : Four Special Matrices
octave:
N = nnz (A)
return the number of nonzero elements in A
матрицы Тоеплитца - матрицы, где все строки - копия первой, но со сдвигом
octave> toeplitz([2 -1 0 0])
ans =
2 -1 0 0
-1 2 -1 0
0 -1 2 -1
0 0 -1 2
octave> k = toeplitz([2 -1 zeros(1,4)])
k =
2 -1 0 0 0 0
-1 2 -1 0 0 0
0 -1 2 -1 0 0
0 0 -1 2 -1 0
0 0 0 -1 2 -1
0 0 0 0 -1 2
octave> ks = sparse(k)
ks =
Compressed Column Sparse (rows = 6, cols = 6, nnz = 16 [44%])
(1, 1) -> 2
(2, 1) -> -1
(1, 2) -> -1
(2, 2) -> 2
(3, 2) -> -1
(2, 3) -> -1
(3, 3) -> 2
(4, 3) -> -1
(3, 4) -> -1
(4, 4) -> 2
(5, 4) -> -1
(4, 5) -> -1
(5, 5) -> 2
(6, 5) -> -1
(5, 6) -> -1
(6, 6) -> 2
определитель матриц с целыми числами - целое число
диагональные верхне- и нижне- треугольные матрицы обратимы, если не имеют
нулей на диагонали
у обратимых матриц нуль-пространство состоит только из нуль-вектора
double diff matrices:
% fixed-fixed
k = toeplitz([2 -1 zeros(1,4)])
2 -1 0 0 0 0
-1 2 -1 0 0 0
0 -1 2 -1 0 0
0 0 -1 2 -1 0
0 0 0 -1 2 -1
0 0 0 0 -1 2
% circular
c = k ; c(1,6) = c(6,1) = -1 ;
c =
2 -1 0 0 0 -1
-1 2 -1 0 0 0
0 -1 2 -1 0 0
0 0 -1 2 -1 0
0 0 0 -1 2 -1
-1 0 0 0 -1 2
такие матрицы например определяют систему пружин, где все они соединены
последовательно и по кругу - конец последней соединен с началом первой
% fixed-free
t = k ; t(1,1) = 1 ;
t =
1 -1 0 0 0 0
-1 2 -1 0 0 0
0 -1 2 -1 0 0
0 0 -1 2 -1 0
0 0 0 -1 2 -1
0 0 0 0 -1 2
% free-fixed
b = k ; b(6,6) = 1 ;
b =
2 -1 0 0 0 0
-1 2 -1 0 0 0
0 -1 2 -1 0 0
0 0 -1 2 -1 0
0 0 0 -1 2 -1
0 0 0 0 -1 1
% free-free
f = k ; f(1,1) = f(6,6) = 1 ;
f =
1 -1 0 0 0 0
-1 2 -1 0 0 0
0 -1 2 -1 0 0
0 0 -1 2 -1 0
0 0 0 -1 2 -1
0 0 0 0 -1 1
эта матрица - сингулярна
например при преобразованиях Фурье :
cos матрицы "free-free"
sin матрицы "fix-fix"
Eqns
матрица "двойной разности"
...
-1 2 -1
-1 2 -1
-1 2 -1
...
вычисляет вторую производную для константного вектора, для вектора линейного
и для вектора квадратов, а также для векторов [sin(x), 2.sin(x),3.sin(x), ...] и т.д.
recitation 01
умножение прямоугольной матрицы на транспонированную - всегда дает
симметричную
03 : Solving a Linear System
если матрица sparse - перегруппировать строки чтобы максимальный pivot был
в первой строке, тогда при элиминации метод будет более стабильным, потому что :
element/pivot < 1
и еще надо стараться, чтобы матрицы с ведущими нулями оказались сразу внизу
строки с максимальным pivot (и чем больше таких нулей - тем ниже строка), а еще
ниже уже шли строки не с ведущими нулями
если матрица А имеет структуру, то скорее всего после элиминации матрица А⁻
структуры иметь уже не будет
A \ I = A⁻ патамушта A⁻ * A = I
стобцы A⁻ - это ответ на импульное воздействие на систему, определяемую
матрицей A, а само импульсное воздействие - это мартица I
-----------------------------------------------------------------------------
четыре задачи линейки : четыре декомпозиции матрицы :
-----------------------------------------------------------------------------
elimintation : A = I * U = L * D * L⁺ lu(A)
orthogonalization : A = Q * R qr(A)
eigenvalues : A = S * A * S⁻ eig(A)
singular value : A = U * Σ * V⁻ svd(A)
-----------------------------------------------------------------------------
U - верхне-треугольные матрицы
L - нижне-треугольные
D - диагональные
инверсия симметричной матрицы будет тоже симметричной матрицей
04 : Delta Function
-u'' = δ(x-a) u(0)=0 u(1)=0
/\
/ \
octave> m = [2 -1 0; -1 2 -1; 0 -1 2] ;
octave> m \ [0;1;0] % это сигма-импульс воздействия
ans =
0.50000
1.00000
0.50000
-u'' = δ(x-a) u'(0)=0 u(1)=0
----
\
\
octave> m = [1 -1 0; -1 2 -1; 0 -1 2] ;
octave> m \ [0;1;0] % это сигма-импульс воздействия
ans =
2
2
1
recitation 02
-u'' = δ(x-a) - δ(x-b) u'(0)=0 u'(1)=0
---
\
\
-----
octave> m = [1 -1 0; -1 2 -1; 0 -1 1] ;
octave> m \ [0;1;-1] % это сигма-импульс воздействия
warning: matrix singular to machine precision
ans =
0.33333
0.33333
-0.66667
octave> r2 = [3;3;3] + ans
r2 =
2.6667
2.6667
3.6667
octave:80> m * r2
ans =
0.00000
-1.00000
1.00000
octave> m = [1 -1 0; -1 2 -1; 0 -1 2] ;
octave> inv(m)
ans =
3 2 1
2 2 1
1 1 1
octave> m = [2 -1 0; -1 2 -1; 0 -1 2] ;
octave> inv(m)
ans =
0.75000 0.50000 0.25000
0.50000 1.00000 0.50000
0.25000 0.50000 0.75000
beam equation : u'''' = f(x)
u(0)=0, u''(0)=0, u(1)=0, u''(0)=0
fourth difference matrix:
1 -4 6 -4 1 0 0 0 ...
0 1 -4 6 -4 1 0 0 ...
0 0 1 -4 6 -4 1 0 ...
...
05 : : Eigenvalues (Part 1)
ни один вектор не может иметь вещественных собственных значений, если матрица
воздействия - это матрица вращения, но собственные значения в таком случае
будут комплексными числами
симметричные вещественные матрицы будут всегда иметь вещественные
собственные значения и ортогональные собственные векторы
антисимметричные вещественные матрицы будут всегда иметь чисто мнимые
собственные значения
чтобы решить дифур u' = A * u , u(0) = z надо
1) найти собсвенные значения матрицы А и ее собственные векторы
2) разложить z = c₁*v₁ + ... + cₙ*vₙ, где v - собственные векторы матрицы А
3) тогда u' = c₁*exp(l₁*t)*v₁ + ... + cₙ*exp(lₙ*t)*vₙ,
где li - собственные значения матрицы А
06 : Eigenvalues (Part 2)
A = S * Λ * S⁻ , где
S - матрица из столбцов собственных векторов матрицы А
Λ - диагональная матрица собственных чисел матрицы А
такое разложение называется "диагонализацией"
A⁻ = S * Λ⁻ * S⁻
если одно из собственных значений равно нулю, значит матрица необратима
для симметричных матриц А все собственные числа - вещественны, все вектора
ортогональны и могут быть нормированы и
S⁺ * S = I
а это значит, что
S⁺ = S⁻
"fixed-fixed" difference matrix:
octave> k = toeplitz ([2,-1 zeros(1,3)])
k =
2 -1 0 0 0
-1 2 -1 0 0
0 -1 2 -1 0
0 0 -1 2 -1
0 0 0 -1 2
octave> output_precision (2)
octave> [V,L,W] = eig (k)
V =
2.9e-01 -5.0e-01 5.8e-01 5.0e-01 -2.9e-01
5.0e-01 -5.0e-01 1.1e-17 -5.0e-01 5.0e-01
5.8e-01 2.3e-16 -5.8e-01 -1.6e-16 -5.8e-01
5.0e-01 5.0e-01 1.8e-16 5.0e-01 5.0e-01
2.9e-01 5.0e-01 5.8e-01 -5.0e-01 -2.9e-01
L =
Diagonal Matrix
0.27 0 0 0 0
0 1.00 0 0 0
0 0 2.00 0 0
0 0 0 3.00 0
0 0 0 0 3.73
W =
2.9e-01 -5.0e-01 5.8e-01 5.0e-01 -2.9e-01
5.0e-01 -5.0e-01 1.1e-17 -5.0e-01 5.0e-01
5.8e-01 2.3e-16 -5.8e-01 -1.6e-16 -5.8e-01
5.0e-01 5.0e-01 1.8e-16 5.0e-01 5.0e-01
2.9e-01 5.0e-01 5.8e-01 -5.0e-01 -2.9e-01
07 : Positive Definite
у PD-матрицы все "угловые определители" должны быть положительными
например, для матрицы 3х3:
[ a b c
b d e
f g h ]
a > 0
a*d - b*b > 0
значение первого pivot будет a
значение второго : (a*d - b*b)/a
значение третьего : ???
энергия: x⁺ * K * x > 0 if x > 0
где K - nxn positive-definite matrix
графики для I:
[1 3 ; 3 9] :
[1 5 ; 5 9] :
x⁺ * K * x = x⁺ * A⁺ * A * x = (A * x)⁺ * (A * x) > 0
равенство достигается только если A * x = 0
например, free-free и circular матрицы - PSD-матрицы : x⁺ * K * x = 0
матрица К PD-матрица <==> матрица K имеет независимые столбцы
инверсия PD-матрицы - положительно-определенная матрица
сумма двух положительно-определенных - положительно-определенная матрица
recitation 03
если H и K - PD-матрицы, тогда блочная матрица [H zeros; zeros K] будет иметь те же самые
pivots и те же самые собственные числа, а собственные вектора - просто дополнятся нулями,
снизу - для векторов H и сверху - для векторов K
матрица [K K ; K K] - необратимая, поскольку имеет одинаковые столбцы и одинаковые строки
ее нуль-пространство - большое
если мы сложим PD-матрицу с единичной, умноженной на число, тогда
(K + n*I) * x = (l + n) * x
т.е. все собственные числа увеличатся на число n
octave:
[R, P] = chol (A)
compute the upper Cholesky factor, R, of the real symmetric
or complex Hermitian positive definite matrix A
R' * R = A
called with two or more output arguments P flags whether the matrix
A was positive definite and 'chol' does not fail
a zero value of P indicates that matrix A is positive definite and
R gives the factorization. otherwise, P will have a positive value
08 : Springs and Masses
... skipped ...
recitation 04
уравение диффузии с конвекцией:
-u'' + v * u' = δ , u(0) = 0 , u(1) = 0
вторая производная - это процесс диффузии, а первая производная - процесс конвекции
элементарная матрица, это произведение вектор-столбца на вектор-строку, когда в первом
лишь один элемент ненулевой (в позиции i) и во втором лишь один элемент ненулевой
(в позиции j). результат - матрица, где на месте ij стоит элемент из вектор-столбца
сложение такой матрицы с матрицей I дает операцию, которая называется транзвекцией
SVD : A * v = σ * u , σ ∈ ℝ₊
mxn nx1 mx1
все такие v - ортогональны друг другу и все такие u тоже ортогональны друг другу
даже если матрица А - квадратная, но сингулярная - все равно можно найти такие v и u
которые будут отвечать условиям ортогональности
уравнение говорит о том, что для матрицы А найдется базис в пространстве размерности n
такой что при преобразовании А он перейдет в базис другого пространства - размерности u
оказывается что при этом
v - это собственные векторы матрицы A⁺*A и
u - собственные векторы матрицы A*A⁺
таким образом, если матрица А квадратная, положительно-определенная, то v и u совпадают
с собственными векторами А, σ совпадает с λ
в октаве:
S = svd (A)
[U, S, V] = svd (A)
du/dt = f(u,t)
если случай линейный то матрица G в уравнении un+1 = G * un постоянна и решение может быть
найдено как
un = Gⁿ * u0
и тогда собственные числа матрицы G дадут ответ - стабильно решение или нет
но функция f может быть нелинейной
Forward:
Un+1 - Un
--------- = f (Un, n*Δt)
Δt
Backward:
Un+1 - Un
--------- = f (Un+1, (n+1)*Δt)
Δt
для нелинейных функций существует метод Ньютона, который вычисляет Un+1 итерациями
в октаве есть функция, которая вычисляет решение методом Рунге-Кутта (четвертого порядка):
[T, Y] = ode45 (FUN, TRANGE, INIT)
SOLUTION = ode45 (...)
Solve a set of non-stiff Ordinary Differential Equations (non-stiff
ODEs) with the well known explicit Dormand-Prince method of order
4.
FUN is a function handle, inline function, or string containing the
name of the function that defines the ODE: 'y' = f(t,y)'. The
function must accept two inputs where the first is time T and the
second is a column vector of unknowns Y.
TRANGE specifies the time interval over which the ODE will be
evaluated. Typically, it is a two-element vector specifying the
initial and final times ('[tinit, tfinal]'). If there are more
than two elements then the solution will also be evaluated at these
intermediate time instances.
The function typically returns two outputs. Variable T is a column
vector and contains the times where the solution was found. The
output Y is a matrix in which each column refers to a different
unknown of the problem and each row corresponds to a time in T.
The output can also be returned as a structure SOLUTION which has a
field X containing a row vector of times where the solution was
evaluated and a field Y containing the solution matrix such that
each column corresponds to a time in X. Use
'fieldnames (SOLUTION)' to see the other fields and additional
information returned.
09 : Oscillation
-- [V, LAMBDA, W] = eig (A, B)
compute the eigenvalues (LAMBDA) and optionally the right
eigenvectors (V) and the left eigenvectors (W)
of a pair of matrices B⁻*A
10: Finite Differences in Time
Trapezoidal:
Un+1 - Un A
--------- = - (Un+1 + Un)
∆t 2
BDF2 :
Un+1 - Un 1 Un+1 -2Un + Un-1
--------- + - ---------------- = A * Un+1
∆t 2 ∆t
11: Least Squares
LSS: A * u = b
mxn nx1 mx1
m > n
тогда решение LSS будет:
A⁺ * A * u = A⁺ * b
nxm mxn nx1 nxm mx1
как получить эту формулу?
столбцы матрицы А определяют плоскость, в которой лежит проекция p вектора b,
а вектор ошибок
e = b - p = A * ū
перпендикулярен этой плоскости, что означает что любой столбец матрицы A
при умножении на вектор e дает 0
здесь ū - это вектор, лежащий в плоскости пространства столбцов прямоугольной
матрицы A. в матричной форме:
A⁺ * e = 0
A⁺ * (b - A * ū) = 0
A⁺ * A * ū = A⁺ * b
вектор проекции вектора b на пространство столбцов матрицы А может быть выражен как
p = b - e = A * ū
A⁺ * A * ū = A⁺ * b
вставляя в первую формулу выражение для ū получаем:
p = A * (A⁺ * A)⁻ * A⁺ * b
p = P * b
P = A * (A⁺ * A)⁻ * A⁺
если бы матрица А была квадратной, то P = I т.е. проекция p совпала бы с b и
вектор ошибок e был бы нулевым
P² = P
применяя этот метод мы должны избегать ситуации, когда A * u = 0 патамушта в этом
случае матрица A⁺ * A будет сингулярной, а значит - необратимой!
в статистике часто используют "взвешенный метод"
w = C * e
A * ú = w
где С - диагональная матрица ненулевые элементы которой обратны к девиациям датчиков σ²
и, тем самым, слабее учитывают сильные оклонения от среднего и больше учитыват слабые
отклонения
разумеется, в случае диагональной матрицы предполагается независимость показаний датчиков
друг от друга (отсутствие корреляций)
12: Graphs and Networks
полный граф - где все вершины попарно соединены
дерево - граф с вершинами, которые могут иметь только одно соединение (листы) или несколько
соединений (ветви) и если граф направленный, то существует вершина, в которую не входят
ребра, а только выходят из нее - основание дерева (root)
13: Kirchhoff's Current Law (KCL)
закон Киргоффа для токов KCL :
закон Киргоффа для напряжений (KVL) :
A * u = v сумма компонент v для замкнутого контура равна нулю
A - матрица сопряжений
recitation 05
антисимметричные матрицы (A⁺=-A) имеют чисто мнимые собственные числа и ортогоналные
собственные векторы - такие системы сохраняют энергию
две последовательности, для которых нужно знать ряд Тайлора наизусть:
- экспонента x = 1 + x + x²/2 + x³/6 + ...
- геометрическая прогрессия 1/(1-x) = 1 + x² + x³ + ...
поскольку в последнем случае x должен быть меньше 1, то это становится критерием
стабильности трапезоидного метода - наш ∆t*λ/2 тоже должен быть меньше 1,
а это значит, что если матрица имеет condition number большой, то для
сохранения стабильности метода мы должны брать ОЧЕНЬ маленький шаг ∆t
octave:
-- cond (A)
-- cond (A, p)
compute the P-norm condition number of a matrix with respect to inversion
'cond (A)' is defined as 'norm (A, p) * norm (inv (A), p)'
by default, 'p = 2' is used which implies a (relatively slow)
singular value decomposition. other possible selections are
'p = 1, Inf, "fro"' which are generally faster
recitation 06
du/dt = A * u
|u(t)| const если A⁺ = -A
trapezoide : |Un| остается неизменной поскольку |Un+1|² = Un+1⁺ * Un+1 = Un⁺ * Un
собственные векторы образуют ортогональный базис если A⁺ * A = A * A⁺ т.е. когда
эти две матрицы коммутируют
-d/dx (C*du/dx) = f(x) соответствует матричная форма (1/∆t)² * A⁺ * C * A * u = f
где A - матрица разностей, С - диагональная матрица
тогда A⁺ * A - это матрица двойных разностей
----------- z
|
dz/dx = δ(x-a) x---------
0 a
14: exam review
"tall" matrix = "thin" matrix
LSS: A⁺ * A * u = A⁺ * b
A = Q * R Q - orthogonal, R - upper-triangular
R⁺ * Q⁺ * Q * R * u = R⁺ * Q⁺ * b
Q⁺ * Q = I
R * u = Q * b
octave:
-- [Q, R] = qr (A)
-- [Q, R, P] = qr (A) # non-sparse A
-- X = qr (A) # non-sparse A
compute the QR factorization of A, using standard LAPACK subroutines
the QR factorization is Q * R = A
where Q is an orthogonal matrix and R is upper triangular
if the matrix A is full, and a third output P is requested, then
'qr' calculates the permuted QR factorization
Q * R = A * P , where
Q is an orthogonal matrix,
R is upper triangular, and
P is a permutation matrix.
the permuted QR factorization has the additional property that
the diagonal entries of R are ordered by decreasing magnitude
in other words, 'abs (diag (R))' will be ordered from largest to smallest
абсолютная ошибка вычисления уравнения K * u = b оценивается как
∆u = ∆b/λ₁ , где ∆b - максимальная погрешность измерений
такая оценка не учитывает накопленной погрешности при преобразованиях K
относительная ошибка |∆u|/|u| оценивается как
ε * λₙ/λ₁ , где ε - инструментальная ошибка измерений
т.е. такая ошибка пропорциональна cond(K). cond(K) ≈ n², где n - размерность u
если матрица несимметричная, да еще и с большими значениями выше диагонали и с
маленькими - ниже (или наоборот), тогда
cond (A) = cond (A⁺ * A)
например A = [1 10⁶ ; 0 1] имеет собственными числами 1 и 1, но ее cond равен 10¹²
15 : trusses and A⁺*C*A
16 : trusses (part 2)
recitation 06
17: Finite Elements in 1D (part 1)
- d/dx = (d/dx)⁺
(A * u) * w = u⁺ * A⁺ * w
mxn nx1 mx1 1xn nxm mx1
(A * u)⁺ * w = u⁺ * A⁺ * w
mxn nx1 mx1 1xn nxm mx1
граничные условия и матрицы:
u f
A * u ↓ ↑ A⁺ w
e → C → w
free-fixed task:
A = [ -1 1 0 [ u0 A⁺ = [ -1 0 0 [ w1
0 -1 1 u1 1 -1 0 w2
0 0 -1 ] u2 ] 0 1 -1 ] w3 ]
u3 = 0 w0 = 0
Dirichle problem - fixed-fixed
Neumann problem - free-free
staticly determinate случай - это когда матрица А - квадратная, тогда по стрелкам
можно пойти и в обратную сторону: u = (A)⁻ * (C⁻) * (A⁺)⁻ * f
если при диффуре первой степени заданы оба начальных условия, то матица А перестает
быть квадратной и такой случай называют indeterminate
18 : Finite Elements in 1D (part 2)
linear approximation by hat-functions with slopes ±1 and duration 2∆t
piece-wize linear approximation
19 : Quadratic/Cubic Elements
quadratic approximation:
cubic approximation :
TODO
20 : Element Matrices; 4th Order Bending Equations
∫ c(x) * (du/dt)² dx = u⁺ * K * u ,
где К - матрица суммирования для выражения u = u₀*Ф₀ + u₁*Ф₁ + ...
при c(x) = 1
первая производная - склон
вторая производная - кривизна
recitation 08
образцы кода для матлаба можно посмотреть тут https://math.mit.edu/cse
hat-functions don't have the second deriatives
нельзя проитегрировать произведение двух разных дельта-функций
потомушта такой интеграл равен ∞
21 : Boundary Conditions, Splines, Gradient, Divergence
если число изолированных точек становится большим (скажем, больше 6) то
интерполировать с помощью полиномов - не слишком хорошая идея. в таком
случае можно применить cubic_splines-интерполяцию
в октаве:
-- PP = splinefit (X, Y, BREAKS)
-- PP = splinefit (X, Y, P)
-- PP = splinefit (..., "periodic", PERIODIC)
-- PP = splinefit (..., "robust", ROBUST)
-- PP = splinefit (..., "beta", BETA)
fit a piecewise cubic spline with breaks (knots) BREAKS to the noisy data, X and Y
X is a vector, and Y is a vector or N-D array
if Y is an N-D array, then X(j) is matched to Y(:,...,:,j)
P is a positive integer defining the number of intervals along X,
and P+1 is the number of breaks
the number of points in each interval differ by no more than 1
the optional property PERIODIC is a logical value which specifies
whether a periodic boundary condition is applied to the spline
the length of the period is 'max (BREAKS) - min (BREAKS)'
the default value is 'false'.
the optional property ROBUST is a logical value which specifies if
robust fitting is to be applied to reduce the influence of outlying
data points. three iterations of weighted least squares are
performed. weights are computed from previous residuals
the sensitivity of outlier identification is controlled by the property BETA
0 < BETA < 1. the default value is 1/2
Values close to 0 give all data equal weighting
increasing values of BETA reduce the influence of outlying data
Values close to unity may cause instability or rank deficiency
the fitted spline is returned as a piecewise polynomial, PP,
and may be evaluated using 'ppval'
22 : Gradient and Divergence
a*x² + b*x*y + c*y² = 1
[ a b
b c ]
если эта матрица положительно-определена, то уравнение будет эллиптическим
если матрица неопределенная, то уравнение будет гиперболическим
если определитель матрицы равен нулю, то уравнение будет параболическим
чтобы определить, является ли векторное поле v градиентом какого-нибуть
скалярного поля u(x,y) надо проверить, что
∂v/∂x - ∂v/∂y = 0
если это выполняется, то поле v является градиентом
второй тест - это тест по любому замкнутому контуру C:
⨚ v₁dx + v₂dy = 0
C
это аналог KVL - "напряжение на замкнутом контуре равно нулю"
есть и аналог для KCL: ∇ . w = 0
"сколько тока в узел втекает, столько из узла и вытекает"
23 : Laplace's Equation
пусть есть функция u(x,y)
уравение Лапласа (Laplace) :
div (grad u) = 0
или:
∂²u/(∂x)² + ∂²u/(∂y)² = 0
уравенине Пуассона (Poisson) :
div (grad u) = f
или:
∂²u/(∂x)² + ∂²u/(∂y)² = f(x,y)
уравнение Лапласа - это уравнение "без источников"
уравнение Пуассона - "с источниками"
оба - изотропны
уравения Коши-Риманна (Cauchy-Riemann) для
- функции потенциала u(x,y) и
- функции потока s(x,y) :
∂u/∂x = ∂s/∂y
∂u/∂y = - ∂s/∂y
т.е.:
grad s = A * grad u
A = [ 0 -1
1 0 ]
пусть z = x + i * y, тогда:
в декартовых в полярных
------------------------------------- ------------------------------
Re (zⁿ) Im (zⁿ) or Re (zⁿ) Im (zⁿ)
u s
----------------------------------- ------------------------------
z¹ x y r *cos( Θ) r *sin( Θ)
z² x² - y² 2*x*y r²*cos(2Θ) r²*sin(2Θ)
z³ x³ - 3*x*y² 3*x²*y - y³ r³*cos(3Θ) r³*sin(2Θ)
... ........... ...........
-------------------------------------- ------------------------------
все эти функции удовлетворяют Коши-Риманну
recitation 09
градиентные поля не содержат "вращений" - если поместить туда частичку, то она
будет двигаться, но не будет при этом вращаться. таковы все потенциальные поля -
поле точечного заряда, как пример
24 : Laplace's Equation (part 2)
пусть z = x + i * y, тогда f(zⁿ) удовлетворяет уравнению Лапласса:
fxx = n * (n-1) * (1 + i*y)^(n-2)
fyy = n * (n-1) * (1 + i*y)^(n-2) * i² = - (n-1) * (1 + i*y)^(n-2)
fxx + fyy = 0
и вообще, любая функция f(z) будет удовлетворять уравнению Лапласа
в полярной системе координат любая функция f(r*exp(i*Θ)) тоже будет удовлетворять
уравнению Лапласа
функция exp(z) аналитична всюду, а ее вещественная и мнимая части являются
гармоническими функциями
функция ln(z) аналитична всюду, кроме точки (0,0)
ln(z) = ln(x + i*y) = ln(r*exp(i*Θ)) = ln(r) + i*Θ = √ x² + y² + i * arctg(y/x)
ее вещественная и мнимая части - гармонические функции u и s
функция Грина: это решение уравнения Пуассона, когда правая часть есть δ
а левая берется из аналитической функции ln(z)
uxx + uyy = 2πδ
u = ln(r)/2π
25 : Fast Poisson Solver (part 1)
формула Грина:
(grad u, w) = (u, -div w) + ⨚ u * (w , n) ds
S
где скобки означают "скалярное произведение",
а n - нормаль к поверхности интегрирования S
grad⁺ = - div
u(x,y) - скаляр
w = A * u = [∂u/∂x ; ∂u/∂y]
w - вектор
fixed (Dirichlet) : u=u0 на части (или на всей) границе S
free (Noemann) : w.n на части (или на всей) границе S
26 : Fast Poisson Solver (part 2); Finite Elements in 2D
квадратный полигон, граничные условия - u0, K * u = f
матрица двойных разностей:
octave:12> k = toeplitz([4 -1 0 -1 zeros(1,5)])
k =
4 -1 0 -1 0 0 0 0 0
-1 4 -1 0 -1 0 0 0 0
0 -1 4 -1 0 -1 0 0 0
-1 0 -1 4 -1 0 -1 0 0
0 -1 0 -1 4 -1 0 -1 0
0 0 -1 0 -1 4 -1 0 -1
0 0 0 -1 0 -1 4 -1 0
0 0 0 0 -1 0 -1 4 -1
0 0 0 0 0 -1 0 -1 4
octave:13> det (k)
ans = 69388.00000
или
octave> a = toeplitz([2 -1 0])
a =
2 -1 0
-1 2 -1
0 -1 2
octave> b = eye(3)
b =
Diagonal Matrix
1 0 0
0 1 0
0 0 1
octave> kron(b,a) + kron(a,b)
ans =
4 -1 0 -1 0 0 0 0 0
-1 4 -1 0 -1 0 0 0 0
0 -1 4 0 0 -1 0 0 0
-1 0 0 4 -1 0 -1 0 0
0 -1 0 -1 4 -1 0 -1 0
0 0 -1 0 -1 4 0 0 -1
0 0 0 -1 0 0 4 -1 0
0 0 0 0 -1 0 -1 4 -1
0 0 0 0 0 -1 0 -1 4
------- случай размерности 1
уравение Фурье fixed-fixed :
-u'' = λ * u , u(0) = 0 = u(1)
"собственные функции" - sin (πkx)
"собственные числа" - λ = π²k²
yk = [sin 1πkh, sin 2πkh, ..., sin Nπkh]
λk = 2 - 2 * cos πkh > 0
т.о. матрица K будет положительно-определенной
если уравнение free-free, то одно из собственных значений будет равно 0 и
матрица К будет полуопределенной
------- случай размерности 2
уравение Фурье fixed-fixed :
-uxx - uyy = λ*u , u(0) задан по всей границе
u(x,y) = (sin πkx) * (sin πly)
λkl = π²k² + π²l²
что касается К, то команда октава 'kron(b,a) + kron(a,b)' все сделает за нас
b - единичная матрица умноженная на 2
a - матрица двойных разностей для одной "молекулы" (минимальная матрица Тоеплитца)
(см. пример выше и рисунок выше)
K * u = f
1) находим собственные векторы ykl и собственные числа λkl
2) выражаем f = Σ ckl * ykl , как линейную комбинацию собственных векторов с
некими коэффициентами сkl
3) тогда решением будет : u = Σ ckl/λkl * ykl
слабая интегральная форма уравнения Пуассона :
∫∫ (∂u/∂x * ∂v/∂x + ∂u/∂y * ∂v/∂y) dx dy = ∫∫ f(x,y) * v(x,y) dx dy
где v(x) - тестовая функция
u(x) - искомая функция
идея финитных элементов заключается в том, чтобы подобрать тестовые функции Ф(x)
так, чтобы
n
u = Σ uk * Фk(x)
k=1