Cel: rozwinięcie idei aproksymacji \( f \) przez \( u \), czyli rozwiązywania problemu:
$$ u = f $$
na rozwiązywanie \( \PDE \), np:
$$ -u'' + bu = f,\quad u(0)=C,\ u'(L)=D $$
Szczególny nacisk położony na metodę Galerkina.
$$
\begin{equation*}
\mathcal{L}(u) = 0,\quad x\in\Omega \end{equation*}
$$
Przykłady (dla zagadnień 1D):
$$
\begin{align*}
\mathcal{L}(u) &= \frac{d^2u}{dx^2} - f(x),\\
\mathcal{L}(u) &= \frac{d}{dx}\left(\dfc(x)\frac{du}{dx}\right) + f(x),\\
\mathcal{L}(u) &= \frac{d}{dx}\left(\dfc(u)\frac{du}{dx}\right) - au + f(x),\\
\mathcal{L}(u) &= \frac{d}{dx}\left(\dfc(u)\frac{du}{dx}\right) + f(u,x)
\end{align*}
$$
$$
\begin{equation*}
\mathcal{B}_0(u)=0,\ x=0,\quad \mathcal{B}_1(u)=0,\ x=L
\end{equation*}
$$
Przykłady:
$$
\begin{align*}
\mathcal{B}_i(u) &= u - g,\quad &\hbox{warunek Dirichleta}\\
\mathcal{B}_i(u) &= -\dfc \frac{du}{dx} - g,\quad &\hbox{warunek Neumanna}\\
\mathcal{B}_i(u) &= -\dfc \frac{du}{dx} - h(u-g),\quad &\hbox{warunek Robina}
\end{align*}
$$
Nowe zagadnienia:
Wstawiając \( u=\sum_jc_j\baspsi_j \) w \( \mathcal{L}=f \) otrzymuje się residuum \( R \).
(\( \mathcal{L}(\uex) - f = 0 \), ale \( \mathcal{L}(u) - f \neq 0 = R \))
$$
\begin{equation*}
\mathcal{L}(u) = \mathcal{L}(\sum_j c_j \baspsi_j) = R \neq 0
\end{equation*}
$$
Cel: minimalizacja \( R \) w funkcji \( \sequencei{c} \) (mamy nadzieję, że w ten sposób \( e \) również będzie mały)
$$ R=R(c_0,\ldots,c_N; x)$$
Pomysł: minimalizacja
$$
\begin{equation*}
E = ||R||^2 = (R,R) = \int_{\Omega} R^2 dx
\end{equation*}
$$
Jeśli minimalizujemy względem \( \sequencei{c} \) to:
$$
\frac{\partial E}{\partial c_i} =
\int_{\Omega} 2R\frac{\partial R}{\partial c_i} dx = 0\quad
\Leftrightarrow\quad (R,\frac{\partial R}{\partial c_i})=0,\quad
i\in\If
$$
\( N+1 \) równań i \( N+1 \) niewiadomych \( \sequencei{c} \)
Pomysł: niech \( R \) będzie ortogonalne do \( V \),
$$
(R,v)=0,\quad \forall v\in V
$$
co jest równoważne
$$
(R,\baspsi_i)=0,\quad i\in\If
$$
\( N+1 \) równań i \( N+1 \) niewiadomych \( \sequencei{c} \)
$$
(R,v)=0,\quad \forall v \in W
$$
Jeśli \( \{w_0,\ldots, w_N \} \) to baza dla \( W \):
$$
(R, w_i )=0,\quad i \in \If
$$
Pomysł: Niech \( R=0 \) w \( N+1 \) punktach obszaru
$$ R(\xno{i}; c_0,\ldots,c_N)=0,\quad i\in\If$$
Metoda kolokacji to metoda residuów ważonych gdzie wagi to delty Diraca
$$ 0 = \int_\Omega R(x;c_0,\ldots,c_N)
\delta(x-\xno{i})\dx = R(\xno{i}; c_0,\ldots,c_N)$$
$$
\hbox{własność } \delta(x):\quad
\int_{\Omega} f(x)\delta (x-\xno{i}) dx = f(\xno{i}),\quad \xno{i}\in\Omega
$$
Podanie przykładów zastosowania metod: najmniejszych kwadratów, Galerkina, kolokacji; do rozwiązywania zagadnień 1D z globalnymi funkcjami bazowymi.
$$ -u''(x) = f(x),\quad x\in\Omega=[0,L],\quad u(0)=0,\ u(L)=0$$
Funkcje bazowe:
$$ \baspsi_i(x) = \sinL{i},\quad i\in\If$$
Residuum:
$$
\begin{align*}
R(x;c_0,\ldots,c_N) &= u''(x) + f(x),\nonumber\\
&= \frac{d^2}{dx^2}\left(\sum_{j\in\If} c_j\baspsi_j(x)\right)
+ f(x),\nonumber\\
&= -\sum_{j\in\If} c_j\baspsi_j''(x) + f(x)
\end{align*}
$$
Ponieważ \( u(0)=u(L)=0 \), należy zapewnić, aby dla wszystkich f-cji bazowych \( \baspsi_i(0)=\baspsi_i(L)=0 \). Jeśli tak to:
$$ u(0) = \sum_jc_j{\color{red}\baspsi_j(0)} = 0,\quad
u(L) = \sum_jc_j{\color{red}\baspsi_j(L)} =0 $$
$$
(R,\frac{\partial R}{\partial c_i}) = 0,\quad i\in\If
$$
$$
\begin{equation*}
\frac{\partial R}{\partial c_i} =
\frac{\partial}{\partial c_i}
\left(\sum_{j\in\If} c_j\baspsi_j''(x) + f(x)\right)
= \baspsi_i''(x)
\end{equation*}
$$
Ponieważ:
\small
$$
\frac{\partial}{\partial c_i}\left(c_0\baspsi_0'' + c_1\baspsi_1'' + \cdots +
c_{i-1}\baspsi_{i-1}'' + {\color{red}c_i\baspsi_{i}''} + c_{i+1}\baspsi_{i+1}''
+ \cdots + c_N\baspsi_N'' \right) = \baspsi_{i}''
$$
$$
\begin{equation*}
(\sum_j c_j \baspsi_j'' + f,\baspsi_i'')=0,\quad i\in\If
\end{equation*}
$$
Niewiadome na lewo, dane na prawo:
$$
\begin{equation*}
\sum_{j\in\If}(\baspsi_i'',\baspsi_j'')c_j = -(f,\baspsi_i''),\quad i\in\If \end{equation*}
$$
Co jest równoważne URL:
$$
\begin{equation*} \sum_{j\in\If}A_{i,j}c_j = b_i,\quad i\in\If
\end{equation*}
$$
$$
\begin{align*}
A_{i,j} &= (\baspsi_i'',\baspsi_j'')\nonumber\\
& = \pi^4(i+1)^2(j+1)^2L^{-4}\int_0^L \sinL{i}\sinL{j}\, dx\nonumber\\
&= \left\lbrace
\begin{array}{ll} {1\over2}L^{-3}\pi^4(i+1)^4 & i=j \\ 0, & i\neq j
\end{array}\right.
\\
b_i &= -(f,\baspsi_i'') = (i+1)^2\pi^2L^{-2}\int_0^Lf(x)\sinL{i}\, dx
\end{align*}
$$
Ortogonalność – użyteczna własność funkcji bazowych
$$
\begin{equation*}
\int\limits_0^L \sinL{i}\sinL{j}\, dx = \delta_{ij},\quad
\quad\delta_{ij} = \left\lbrace
\begin{array}{ll} \half L & i=j \\ 0, & i\neq j
\end{array}\right.
\end{equation*}
$$
\( \Rightarrow\ (\baspsi_i'',\baspsi_j'') = \delta_{ij} \),
a więc wyłącznie elementy na przekątnej \( \neq = 0 \), dzięki czemu z łatwością można znaleźć rozwiązanie:
$$
\begin{equation*}
c_i = \frac{2L}{\pi^2(i+1)^2}\int_0^Lf(x)\sinL{i}\, dx
\end{equation*}
$$
Rozwiązanie przy pomocy sympy dla \( f(x)=2 \):
from sympy import *
import sys
i, j = symbols('i j', integer=True)
x, L = symbols('x L')
f = 2
a = 2*L/(pi**2*(i+1)**2)
c_i = a*integrate(f*sin((i+1)*pi*x/L), (x, 0, L))
c_i = simplify(c_i)
print c_i
\small
$$
\begin{equation*}
c_i = 4 \frac{L^{2} \left(\left(-1\right)^{i} + 1\right)}{\pi^{3}
\left(i^{3} + 3 i^{2} + 3 i + 1\right)},\quad
u(x) = \sum_{k=0}^{N/2} \frac{8L^2}{\pi^3(2k+1)^3}\sinL{2k}
\end{equation*}
$$
\( c_i \) szybko zanikają: \( c_2 = c_0/27 \), \( c_4=c_0/125 \) - pierwszy wyraz może być całkiem niezłym przybliżeniem:
$$
\begin{equation*} u(x) \approx \frac{8L^2}{\pi^3}\sin\left(\pi\frac{x}{L}\right) \end{equation*}
$$
$$
\begin{equation*}
(u''+f,v)=0,\quad \forall v\in V,
\end{equation*}
$$
po reorganizacji:
$$
\begin{equation*}
(u'',v) = -(f,v),\quad\forall v\in V \end{equation*}
$$
co jest sformułowaniem wariacyjnym zagadnienia opisanego \( \PDE \)
\( \forall v\in V \) jest równoważne \( \forall v\in\baspsi_i \), \( i\in\If \), i ostatecznie
$$
\begin{equation*}
(\sum_{j\in\If} c_j\baspsi_j'', \baspsi_i)=-(f,\baspsi_i),\quad i\in\If \end{equation*}
$$
$$
\begin{equation*}
\sum_{j\in\If}(\baspsi_j'', \baspsi_i) c_j=-(f,\baspsi_i),\quad i\in\If \end{equation*}
$$
Poniważ \( \baspsi_i''\propto -\baspsi_i \),
metoda Galerkinga daje ten sam URL i to samo rozwiązanie (w tym konkretnym przykładzie) co metoda najmniejszych kwadratów.
\( R=0 \) (czyli równanie różniczkowe) musi być spełnione w \( N+1 \) punktach:
$$
\begin{equation*}
-\sum_{j\in\If} c_j\baspsi_j''(\xno{i}) = f(\xno{i}),\quad i\in\If
\end{equation*}
$$
To daje URL \( \sum_j A_{i,j}=b_i \) o współczynnikach:
$$
\begin{equation*} A_{i,j}=-\baspsi_j''(\xno{i})=
(j+1)^2\pi^2L^{-2}\sin\left((j+1)\pi \frac{x_i}{L}\right),
\quad b_i=2
\end{equation*}
$$
Niech: \( N=0 \), \( x_0=L/2 \), wtedy
$$ c_0=2L^2/\pi^2 $$
>>> import sympy as sym
>>> # Computing with Dirichlet conditions: -u''=2 and sines
>>> x, L = sym.symbols('x L')
>>> e_Galerkin = x*(L-x) - 8*L**2*sym.pi**(-3)*sym.sin(sym.pi*x/L)
>>> e_colloc = x*(L-x) - 2*L**2*sym.pi**(-2)*sym.sin(sym.pi*x/L)
>>> # Verify max error for x=L/2
>>> dedx_Galerkin = sym.diff(e_Galerkin, x)
>>> dedx_Galerkin.subs(x, L/2)
0
>>> dedx_colloc = sym.diff(e_colloc, x)
>>> dedx_colloc.subs(x, L/2)
0
# Compute max error: x=L/2, evaluate numerical, and simplify
>>> sym.simplify(e_Galerkin.subs(x, L/2).evalf(n=3))
-0.00812*L**2
>>> sym.simplify(e_colloc.subs(x, L/2).evalf(n=3))
0.0473*L**2
Od tej pory rozwiązanie równania różniczkowego będzie uzyskiwane poprzez sformułowanie słabe otrzymywane w wyniku całkowania przez części (zastosowanie tw. Greena)
$$
\begin{align*}
\int_0^L u''(x)v(x) dx &= - \int_0^Lu'(x)v'(x)dx
+ [vu']_0^L\nonumber\\
&= - \int_0^Lu'(x)v'(x) dx
+ u'(L)v(L) - u'(0)v(0)
\end{align*}
$$
Dlaczego?
Warunki Dirichleta: \( u(0)=C \) and \( u(L)=D \). Niech \( B(x) \) będzie np.
$$ B(x) = \frac{1}{L}(C(L-x) + Dx):\qquad B(0)=C,\ B(L)=D $$
$$
\begin{equation*}
u(x) = B(x) + \sum_{j\in\If} c_j\baspsi_j(x),
\end{equation*}
$$
$$ u(0) = B(0)= C,\quad u(L) = B(L) = D $$
Warunek Dirichleta: \( u(L)=D \). Niech \( B(x) \) będzie równe:
$$ B(x) = D:\qquad B(L)=D $$
$$
\begin{equation*}
u(x) = B(x) + \sum_{j\in\If} c_j\baspsi_j(x),
\end{equation*}
$$
$$ u(L) = B(L) = D $$
Znaczna część literatury dot. \( \FEM{} \) stosuje specjalną notację jeśli chodzi o sformułowania wariacyjne
Znajdź takie \( (u-B)\in V \), że
$$ a(u,v) = L(v)\quad \forall v\in V $$
$$ -u''=f, \quad u'(0)=C,\ u(L)=D,\quad u=D + \sum_jc_j\baspsi_j$$
Sformułowanie wariacyjne (słabe):
$$
\int_{\Omega} u' v'dx = \int_{\Omega} fvdx - v(0)C
\quad\hbox{or}\quad (u',v') = (f,v) - v(0)C
\quad\forall v\in V
$$
W zaproponowanej abstrakcyjnej notacji: znajdź \( (u-B)\in V \) takie, że
$$ a(u,v) = L(v)\quad \forall v\in V$$
a więc
$$ a(u,v) = (u',v'),\quad L(v) = (f,v) -v(0)C $$
Dla form liniowych
$$ L(\alpha_1 v_1 + \alpha_2 v_2)
=\alpha_1 L(v_1) + \alpha_2 L(v_2),
$$
Dla form dwuliniowych
$$
\begin{align*}
a(\alpha_1 u_1 + \alpha_2 u_2, v) &= \alpha_1 a(u_1,v) + \alpha_2 a(u_2, v),
\\
a(u, \alpha_1 v_1 + \alpha_2 v_2) &= \alpha_1 a(u,v_1) + \alpha_2 a(u, v_2)
\end{align*}
$$
W zagadnieniach nieliniowych: Znależć \( (u-B)\in V \) takie, że \( F(u;v)=0\ \forall v\in V \)
$$ a(u,v) = L(v)\quad \forall v\in V\quad\Leftrightarrow\quad
a(u,\baspsi_i) = L(\baspsi_i)\quad i\in\If$$
URL odpowiadający powyższemu równaniu można otrzymać wstawiając do niego \( u = B + \sum_jc_j\baspsi_j \):
$$ a(B + \sum_{j\in\If} c_j \baspsi_j,\baspsi_i)c_j = L(\baspsi_i)\quad i\in\If$$
Ze wzgledu na liniowość,
$$ \sum_{j\in\If} \underbrace{a(\baspsi_j,\baspsi_i)}_{A_{i,j}}c_j =
\underbrace{L(\baspsi_i) - a(B,\baspsi_i)}_{b_i}\quad i\in\If$$
Jeśli \( a \) jest symetryczne: \( a(u,v)=a(v,u) \), wtedy
$$ a(u,v)=L(v)\quad\forall v\in V$$
jest równoważne minimalizacji funkcjonału
$$ F(v) = {\half}a(v,v) - L(v) $$
dla wszystkich \( v\in V \). Czyli
$$ F(u)\leq F(v)\quad \forall v\in V $$
Wyprowadzić sformułowania wariacyjne dla pewnych typów równań różniczkowych w 1D uwzględniających
$$
\begin{equation*}
-\frac{d}{dx}\left( \dfc(x)\frac{du}{dx}\right) = f(x),\quad x\in\Omega =[0,L],\
u(0)=C,\ u(L)=D
\end{equation*}
$$
$$
u(x) = B(x) + \sum_{j\in\If} c_j\baspsi_i(x),\quad
$$
$$ R = -\frac{d}{dx}\left( a\frac{du}{dx}\right) -f $$
co zapisane w postaci wariacyjnej:
$$
(R, v) = 0,\quad \forall v\in V
$$
lub jawnie:
$$
\int_{\Omega} \left(-\frac{d}{dx}\left( \dfc\frac{du}{dx}\right) -f\right)v \dx = 0,\quad \forall v\in V
$$
$$ -\int_{\Omega} \frac{d}{dx}\left( \dfc(x)\frac{du}{dx}\right) v \dx
= \int_{\Omega} \dfc(x)\frac{du}{dx}\frac{dv}{dx}\dx -
\left[\dfc\frac{du}{dx}v\right]_0^L
$$
Ostatni wyraz prawej strony znika ponieważ \( v(0)=v(L)=0 \)
Znaleźć \( (u-B)\in V \) takie, że
$$
\int_{\Omega} \dfc(x)\frac{du}{dx}\frac{dv}{dx}dx = \int_{\Omega} f(x)vdx,\quad
\forall v\in V
$$
Po zastosowaniu zwięzłej notacji:
$$ \underbrace{(\dfc u',v')}_{a(u,v)} = \underbrace{(f,v)}_{L(v)},
\quad \forall v\in V $$
Uwzględniając
$$ a(u,v) = (\dfc u', v'),\quad L(v) = (f,v) $$
można wygenerować URL gdzie:
$$
\begin{align*}
A_{i,j} &= a(\baspsi_j,\baspsi_i) = (\dfc \baspsi_j', \baspsi_i')
= \int_\Omega \dfc \baspsi_j' \baspsi_i'\dx =
\int_\Omega \baspsi_i' \dfc \baspsi_j'\dx = \\
& = a(\baspsi_i,\baspsi_j) = A_{j,i}\\
b_i &= (f,\baspsi_i) - (\dfc B',\baspsi_i) = \int_\Omega (f\baspsi_i -
\dfc L^{-1}(D-C)\baspsi_i')\dx
\end{align*}
$$
\( v=\baspsi_i \) oraz \( u=B + \sum_jc_j\baspsi_j \):
$$
(\dfc B' + \dfc \sum_{j\in\If} c_j \baspsi_j', \baspsi_i') =
(f,\baspsi_i), \quad i\in\If
$$
Porządkując składniki otrzymuje się
$$ \sum_{j\in\If} (\dfc\baspsi_j', \baspsi_i')c_j =
(f,\baspsi_i) + (aL^{-1}(D-C), \baspsi_i'), \quad i\in\If
$$
czyli URL dany przez \( \sum_j A_{i,j}c_j=b_i \) gdzie
$$
\begin{align*}
A_{i,j} &= (a\baspsi_j', \baspsi_i') = \int_{\Omega} \dfc(x)\baspsi_j'(x)
\baspsi_i'(x)\dx\\
b_i &= (f,\baspsi_i) + (aL^{-1}(D-C),\baspsi_i')=
\int_{\Omega} \left(f\baspsi_i + \dfc\frac{D-C}{L}\baspsi_i'\right) \dx
\end{align*}
$$
$$
-u''(x) + bu'(x) = f(x),\quad x\in\Omega =[0,L],\
u(0)=C,\ u'(L)=E
$$
Nowy problem:
Jak postępować?:
$$ u = C + \sum_{j\in\If} c_j \baspsi_i(x)$$
Metoda Galerkina: mnożymy przez funkcję wagową \( v \), całkujemy nad \( \Omega \), całkujemy przez części:
$$ (-u'' + bu' - f, v) = 0,\quad\forall v\in V$$
$$ (u',v') + (bu',v) = (f,v) + [u' v]_0^L, \quad\forall v\in V$$
\( [u' v]_0^L = u'(L)v(L) - u'(0)v(0)= E v(L) \) ponieważ \( v(0)=0 \) a \( u'(L)=E \)
$$ (u',v') + (bu',v) = (f,v) + Ev(L), \quad\forall v\in V$$
$$ (u',v') + (bu',v) = (f,v) + Ev(L), \quad\forall v\in V$$
Ważne spotrzeżenia:
Uogólniony zapis problemu:
$$ a(u,v)=L(v)\quad\forall v\in V$$
gdzie dla
$$ (u',v') + (bu',v) = (f,v) + Ev(L), \quad\forall v\in V$$
mamy
$$
\begin{align*}
a(u,v)&=(u',v') + (bu',v)\\
L(v)&= (f,v) + E v(L)
\end{align*}
$$
Wstawiając \( u=C+\sum_jc_j\baspsi_j \) oraz \( v=\baspsi_i \) do
$$ (u',v') + (bu',v) = (f,v) + Ev(L), \quad\forall v\in V$$
można otrzymać
$$
\sum_{j\in\If}
\underbrace{((\baspsi_j',\baspsi_i') + (b\baspsi_j',\baspsi_i))}_{A_{i,j}}
c_j =
\underbrace{(f,\baspsi_i) + E \baspsi_i(L)}_{b_i},\quad i\in\If
$$
Spostrzeżenie: \( A_{i,j} \) tym razem nie jest symetryczna, ze względu na istnienie wyrazu
$$
(b\baspsi_j',\baspsi_i)=\int_{\Omega} b\baspsi_j'\baspsi_i dx
\neq \int_{\Omega} b \baspsi_i' \baspsi_jdx = (\baspsi_i',b\baspsi_j)
$$
$$ (u',v') + (bu',v) = (f,v) + u'(L)v(L) - u'(0)v(0)$$
Łatwo zapomnieć o uwzględnieniu warunku brzegowego całkując przez części. W ten sposób pomyłkowo przypisujemy warunek \( u'=0 \) na danej części brzegu!
Problem:
$$
\begin{equation*}
-(\dfc(u)u')' = f(u),\quad x\in [0,L],\ u(0)=0,\ u'(L)=E
\end{equation*}
$$
m. Gal.: przemnóż przez \( v \), scałkuj, scałkuj przez części
$$ \int_0^L \dfc(u)\frac{du}{dx}\frac{dv}{dx}\dx =
\int_0^L f(u)v\dx + [\dfc(u)vu']_0^L\quad\forall v\in V
$$
$$ \int_0^L \dfc(u)\frac{du}{dx}\frac{dv}{dx}v\dx =
\int_0^L f(u)v\dx + \dfc(u(L))v(L)E\quad\forall v\in V
$$
lub
$$ (\dfc(u)u', v') = (f(u),v) + \dfc(u(L))v(L)E\quad\forall v\in V
$$
$$
\begin{equation*}
-u''(x)=f(x),\quad x\in \Omega=[0,1],\quad u'(0)=C,\ u(1)=D
\end{equation*}
$$
Sformułowanie wariacyjne: znaleźć \( (u-B)\in V \) takie, aby
$$
(u',\baspsi_i') = (f,\baspsi_i) - C\baspsi_i(0),\ i\in\If
$$
Uwzględniając \( u(x) = B(x) + \sum_{j\in\If}c_j\baspsi_j \) otrzymuje się
$$ \sum_{j\in\If} A_{i,j}c_j = b_i,\quad i\in\If$$
gdzie
$$ A_{i,j} = (\baspsi_j',\baspsi_i')
$$
$$ b_i = (f,\baspsi_i) - (D,\baspsi_i') -C\baspsi_i(0) $$
$$ A_{i,j} = (\baspsi_j',\baspsi_i') = \int_{0}^1 \baspsi_i'(x)\baspsi_j'(x)dx
= \int_0^1 (i+1)(j+1)(1-x)^{i+j} dx
$$
Wybierzmy \( f(x)=2 \):
$$
\begin{align*}
b_i &= (2,\baspsi_i) - (D,\baspsi_i') -C\baspsi_i(0)\\
&= \int_0^1 \left( 2(1-x)^{i+1} - D(i+1)(1-x)^i\right)dx -C\baspsi_i(0)
\end{align*}
$$
Łatwość obliczeń z wykorzystaniem sympy. \( N=1 \) oraz \( \If = \{0,1\} \):
$$
\begin{equation*}
\left(\begin{array}{cc}
1 & 1\\
1 & 4/3
\end{array}\right)
\left(\begin{array}{c}
c_0\\
c_1
\end{array}\right)
=
\left(\begin{array}{c}
-C+D+1\\
2/3 -C + D
\end{array}\right)
\end{equation*}
$$
$$ c_0=-C+D+2, \quad c_1=-1,$$
$$ u(x) = 1 -x^2 + D + C(x-1)\quad\hbox{(ścisłe rozwiązanie)} $$
Niech, poza warunkami brzegowymi Dirichleta, \( \uex \) leży w przestrzeni \( V \), w której poszukiwane jest \( u \) (\( u_e - B \in V \)). Wtedy:
$$
\begin{align*}
u &= B + {\color{red}F},\quad F\in V\\
a(B+F, v) &= L(v),\quad\forall v\in V\\
\uex & = B + {\color{red}E},\quad E\in V\\
a(B+E, v) &= L(v),\quad\forall v\in V
\end{align*}
$$
Odejmując od siebie równania otrzymuje się: \( a(F-E,v)=0 \)
z czego wynika, że:
\( E=F \) oraz \( u = \uex \)
Problem :
Problem w postaci klasycznej:
$$ -u''(x) = 2,\quad x\in (0,L),\ u(0)=u(L)=0,$$
Sformułowanie wariacyjne:
$$ (u',v') = (2,v)\quad\forall v\in V $$
Ponieważ \( u(0)=0 \) oraz \( u(L)=0 \) należy wymusić
$$ v(0)=v(L)=0,\quad \baspsi_i(0)=\baspsi_i(L)=0$$
Wybieramy jako funkcje bazowe funkcje trójkątne: \( \baspsi_i=\basphi_i \), \( i=0,\ldots,N_n-1 \).
Problem: funkcje skrajne nie spełniają warunków brzegowych \( \basphi_0(0)\neq 0 \) oraz \( \basphi_{N_n-1}(L)\neq 0 \)
Rozwiązanie: wykluczamy z bazy \( \basphi_0 \) oraz \( \basphi_{N_n-1} \) i pracujemy na tak okrojonej bazie:
$$ \baspsi_i=\basphi_{i+1},\quad i=0,\ldots,N=N_n-3$$
Wprowadzając stosowną indeksacje \( \nu(i) \): \( \baspsi_i = \basphi_{\nu(i)} \) otrzymuje się:
$$ u = \sum_{j\in\If}c_j\basphi_{\nu(j)},\quad i=0,\ldots,N,\quad \nu(j) = j+1$$
W przypadku numeracji nieregularnej: tablica \( \nu(j) \) będzie bardziej skomplikowana
$$
\begin{equation*}
A_{i,j}=\int_0^L\basphi_{i+1}'(x)\basphi_{j+1}'(x) dx,\quad
b_i=\int_0^L2\basphi_{i+1}(x) dx
\end{equation*}
$$
Wygodnie jest przeprowadzić reindeksację:
\( i+1\rightarrow i \), \( j+1\rightarrow j \)aby otrzymać wzory: \( \basphi_i'\basphi_j' \), zamiast \( \basphi'_{i+1} \basphi'_{j+1} \)
$$
\begin{equation*}
A_{i-1,j-1}=\int_0^L\basphi_{i}'(x)\basphi_{j}'(x) \dx,\quad
b_{i-1}=\int_0^L2\basphi_{i}(x) \dx
\end{equation*}
$$
$$ \basphi_i' \sim \pm h^{-1} $$
$$ A_{i-1,i-1} = h^{-2}2h = 2h^{-1},\quad
A_{i-1,i-2} = h^{-1}(-h^{-1})h = -h^{-1}$$
oraz \( A_{i-1,i}=A_{i-1,i-2} \)
$$ b_{i-1} = 2({\half}h + {\half}h) = 2h$$
$$
\begin{equation*}
\frac{1}{h}\left(
\begin{array}{ccccccccc}
2 & -1 & 0
&\cdots &
\cdots & \cdots & \cdots &
\cdots & 0 \\
-1 & 2 & -1 & \ddots & & & & & \vdots \\
0 & -1 & 2 & -1 &
\ddots & & & & \vdots \\
\vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\
\vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\
\vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\
\vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\
\vdots & & & & &\ddots & \ddots &\ddots & -1 \\
0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 2
\end{array}
\right)
\left(
\begin{array}{c}
c_0 \\
\vdots\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
\vdots\\
c_{N}
\end{array}
\right)
=
\left(
\begin{array}{c}
2h \\
\vdots\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
\vdots\\
2h
\end{array}
\right)
\end{equation*}
$$
Równanie dla \( i \). węzła:
$$
-\frac{1}{h}c_{i-1} + \frac{2}{h}c_{i} - \frac{1}{h}c_{i+1} = 2h
$$
Ponieważ, \( c_i = u(\xno{i+1})\equiv u_{i+1} \), równanie dla \( i-1 \). węzła ma postać:
$$
-\frac{1}{h}c_{i-2} + \frac{2}{h}c_{i-1} - \frac{1}{h}c_{i} = 2h
$$
i jest równoważne
$$
-\frac{1}{h}u_{i-1} + \frac{2}{h}u_{i} - \frac{1}{h}u_{i+1} = 2h
$$
Równanie \( \FDM{} \) dla problemu \( -u''=2 \) ma postać
$$ -\frac{1}{h^2}u_{i-1} + \frac{2}{h^2}u_{i} - \frac{1}{h^2}u_{i+1} = 2 $$
Po pomnożeniu przez \( h \) okazuje się, że:
\( \FEM \) i \( \FDM \) dają równoważne URL w tym zagadnieniu.
(Równania dla węzłów brzegowych również są w tym przypadku identyczne)
$$
\begin{equation*}
A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} \basphi_i'(x)\basphi_j'(x) \dx
= \int_{-1}^1 \frac{d}{dx}\refphi_r(X)\frac{d}{dx}\refphi_s(X)
\frac{h}{2} \dX,
\end{equation*}
$$
$$ \refphi_0(X)=\half(1-X),\quad\refphi_1(X)=\half(1+X)$$
$$ \frac{d\refphi_0}{dX} = -\half,\quad \frac{d\refphi_1}{dX} = \half $$
Z reguły łańcuchowej:
$$ \frac{d\refphi_r}{dx} = \frac{d\refphi_r}{dX}\frac{dX}{dx}
= \frac{2}{h}\frac{d\refphi_r}{dX}$$
$$
\begin{equation*}
A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} \basphi_i'(x)\basphi_j'(x) \dx
= \int_{-1}^1 \frac{2}{h}\frac{d\refphi_r}{dX}\frac{2}{h}\frac{d\refphi_s}{dX}
\frac{h}{2} \dX = \tilde A_{r,s}^{(e)}
\end{equation*}
$$
$$
\begin{equation*}
b_{i-1}^{(e)} = \int_{\Omega^{(e)}} 2\basphi_i(x) \dx =
\int_{-1}^12\refphi_r(X)\frac{h}{2} \dX = \tilde b_{r}^{(e)},
\quad i=q(e,r),\ r=0,1
\end{equation*}
$$
Obliczenia należy wykonać dla indeksów \( r,s \) przyjmujących wszystkie możliwe kombinacje wartości \( 0,1 \), obliczając za każdym razem jeden z elementów lokalnej macierzy oraz wektora:
$$
\begin{equation*}
\tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{rr}
1 & -1\\
-1 & 1
\end{array}\right),\quad
\tilde b^{(e)} = h\left(\begin{array}{c}
1\\
1
\end{array}\right)
\end{equation*}
$$
Przykład:
$$ \tilde A^{(e)}_{0,1} =
\int_{-1}^1 \frac{2}{h}\frac{d\refphi_0}{dX}\frac{2}{h}\frac{d\refphi_1}{dX}
\frac{h}{2} \dX
= \frac{2}{h}(-\half)\frac{2}{h}\half\frac{h}{2} \int_{-1}^1\dX
= -\frac{1}{h}
$$
Dla \( e=0 \) oraz \( e=N_e \):
$$
\tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{r}
1
\end{array}\right),\quad
\tilde b^{(e)} = h\left(\begin{array}{c}
1
\end{array}\right)
$$
Jeden stopień swobody ("węzeł") w skrajnych elementach (\( r=0 \) odpowiada jednemu stopniowi swobody)
4 elementy typu P1:
vertices = [0, 0.5, 1, 1.5, 2]
cells = [[0, 1], [1, 2], [2, 3], [3, 4]]
dof_map = [[0], [0, 1], [1, 2], [2]] # only 1 dof in elm 0, 3
Kod Pythona wykonujący assembling macierzy globalnej:
# Ae[e][r,s]: element matrix, be[e][r]: element vector
# A[i,j]: coefficient matrix, b[i]: right-hand side
for e in range(len(Ae)):
for r in range(Ae[e].shape[0]):
for s in range(Ae[e].shape[1]):
A[dof_map[e,r],dof_map[e,s]] += Ae[e][i,j]
b[dof_map[e,r]] += be[e][i,j]
W wyniku powstaje ten sam URL co w przypadku obliczeń we współrzędnych globalnych.
Niech
Ogólna postać funkcji \( B \):
$$
\begin{equation*}
B(x) = \sum_{j\in\Ifb} U_j\basphi_j(x)
\end{equation*}
$$
Niech dany będzie warunek Dirichleta \( u(\xno{k})=U_k \), \( k\in\Ifb \):
$$
u(\xno{k}) = \sum_{j\in\Ifb} U_j\underbrace{\basphi_j(x)}_{\neq 0
\hbox{ tylko dla }j=k} +
\sum_{j\in\If} c_j\underbrace{\basphi_{\nu(j)}(\xno{k})}_{=0,\ k\not\in\If}
= U_k $$
$$ -u''=2, \quad u(0)=C,\ u(L)=D $$
$$ \int_0^L u'v'\dx = \int_0^L2v\dx\quad\forall v\in V$$
$$ (u',v') = (2,v)\quad\forall v\in V$$
$$
\begin{equation*}
B(x) = \sum_{j\in\Ifb} U_j\basphi_j(x)
\end{equation*}
$$
W tym wypadku \( \Ifb = \{0,N_n-1\} \), \( U_0=C \), \( U_{N_n-1}=D \); \( \baspsi_i \) są funkcjami związanymi z węzłami wewnątrz obszaru \( \basphi_i \):
$$ \baspsi_i = \basphi_{\nu(i)}, \quad \nu(i)=i+1,\quad i\in\If =
\{0,\ldots,N=N_n-3\} $$
$$
\begin{align*}
u(x) &= \underbrace{C\cdot\basphi_0 + D\basphi_{N_n-1}}_{B(x)}
+ \sum_{j\in\If} c_j\basphi_{j+1}\\
&= C\cdot\basphi_0 + D\basphi_{N_n-1} + c_0\basphi_1 + c_1\basphi_2 +\cdots
+ c_N\basphi_{N_n-2}
\end{align*}
$$
Wstawmy \( u = B + \sum_j c_j\baspsi_j \) do sformułowania wariacyjnego:
$$ (u',v') = (2,v)\quad\Rightarrow\quad (\sum_jc_j\baspsi_j',\baspsi_i')
= (2-B',\baspsi_i)\quad \forall v\in V$$
$$
\begin{align*}
A_{i-1,j-1} &= \int_0^L \basphi_i'(x)\basphi_j'(x) \dx\\
b_{i-1} &= \int_0^L (f(x)\basphi_i(x) -
B'(x)\basphi_i'(x))\dx,\quad B'(x)=C\basphi_{0}'(x) + D\basphi_{N_n-1}'(x)
\end{align*}
$$
for \( i,j = 1,\ldots,N+1=N_n-1 \).
Przyczynki pochodzące od warunków brzegowych - od całki \( -\int B'\basphi_i'\dx \): \( C/h \) należy dodać do \( b_0 \), \( D/h \) należy dodać do \( b_N \).
\footnotesize
Dla pierwszego elementu:
$$
\begin{align*}
\tilde b_0^{(1)} = \int_{-1}^1 \left(f\refphi_1 -
C\frac{2}{h}\frac{d\refphi_0}{dX}\frac{2}{h}\frac{d\refphi_1}{dX}\right)
\frac{h}{2} \dX = \\
\frac{h}{2} 2\int_{-1}^1 \refphi_1 \dX
- C\frac{2}{h}(-\frac{1}{2})\frac{2}{h}\frac{1}{2}\frac{h}{2}\cdot 2
= h + C\frac{1}{h}\tp
\end{align*}
$$
Dla ostatniego elementu:
$$
\begin{align*}
\tilde b_0^{N_e} = \int_{-1}^1 \left(f\refphi_0 -
D\frac{2}{h}\frac{d\refphi_1}{dX}\frac{2}{h}\frac{d\refphi_0}{dX}\right)
\frac{h}{2} \dX = \\
\frac{h}{2} 2\int_{-1}^1 \refphi_0 \dX
- D\frac{2}{h}\frac{1}{2}\frac{2}{h}(-\frac{1}{2})\frac{h}{2}\cdot 2
= h + D\frac{1}{h}\tp
\end{align*}
$$
Metoda 2: zawsze wybieraj \( \baspsi_i = \basphi_i \) dla wszystkich \( i\in\If \):
$$
\begin{equation*}
u(x) = \sum_{j\in\If}c_j\basphi_j(x),\quad \If=\{0,\ldots,N=N_n-1\}
\end{equation*}
$$
Wszystkie Wartości \( u \) (również te na brzegu Dirichleta) można traktować jako niewiadome (wymagające obliczenia).
$$ -u''=2,\quad u(0)=0,\ u(L)=D$$
Macierz w postaci identycznej jakby warunku Dirichleta nie było:
$$
\begin{equation*}
\frac{1}{h}\left(
\begin{array}{ccccccccc}
1 & -1 & 0
&\cdots &
\cdots & \cdots & \cdots &
\cdots & 0 \\
-1 & 2 & -1 & \ddots & & & & & \vdots \\
0 & -1 & 2 & -1 &
\ddots & & & & \vdots \\
\vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\
\vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\
\vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\
\vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\
\vdots & & & & &\ddots & \ddots &\ddots & -1 \\
0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 1
\end{array}
\right)
\left(
\begin{array}{c}
c_0 \\
\vdots\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
\vdots\\
c_{N}
\end{array}
\right)
=
\left(
\begin{array}{c}
h \\
2h\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
2h\\
h
\end{array}
\right)
\end{equation*}
$$
$$
\begin{equation*}
\frac{1}{h}\left(
\begin{array}{ccccccccc}
h & 0 & 0
&\cdots &
\cdots & \cdots & \cdots &
\cdots & 0 \\
-1 & 2 & -1 & \ddots & & & & & \vdots \\
0 & -1 & 2 & -1 &
\ddots & & & & \vdots \\
\vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\
\vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\
\vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\
\vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\
\vdots & & & & &\ddots & \ddots &\ddots & -1 \\
0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h
\end{array}
\right)
\left(
\begin{array}{c}
c_0 \\
\vdots\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
\vdots\\
c_{N}
\end{array}
\right)
=
\left(
\begin{array}{c}
0 \\
2h\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
2h\\
D
\end{array}
\right)
\end{equation*}
$$
Na elemencie 0 znane jest \( u \) w węźle lokalnym o indeksie \( r=0 \)
->
Wymieńmy pierwsze równanie dla tego elementu na \( \tilde c_0 = 0 \):
$$
\begin{equation*}
\tilde A^{(0)} =
A = \frac{1}{h}\left(\begin{array}{rr}
h & 0\\
-1 & 1
\end{array}\right),\quad
\tilde b^{(0)} = \left(\begin{array}{c}
0\\
h
\end{array}\right)
\end{equation*}
$$
Na elemencie \( N_e \) znane jest \( u \) w węźle lokalnym o indeksie \( r=1 \)
->
Wymieńmy ostatnie równanie dla tego elementu na \( \tilde c_1 = D \):
$$
\begin{equation*}
\tilde A^{(N_e)} =
A = \frac{1}{h}\left(\begin{array}{rr}
1 & -1\\
0 & h
\end{array}\right),\quad
\tilde b^{(N_e)} = \left(\begin{array}{c}
h\\
D
\end{array}\right)
\end{equation*}
$$
Algorytm:
$$
\begin{equation*}
\frac{1}{h}\left(
\begin{array}{ccccccccc}
h & 0 & 0
&\cdots &
\cdots & \cdots & \cdots &
\cdots & 0 \\
0 & 2 & -1 & \ddots & & & & & \vdots \\
0 & -1 & 2 & -1 &
\ddots & & & & \vdots \\
\vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\
\vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\
\vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\
\vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\
\vdots & & & & &\ddots & \ddots &\ddots & 0 \\
0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h
\end{array}
\right)
\left(
\begin{array}{c}
c_0 \\
\vdots\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
\vdots\\
c_{N}
\end{array}
\right)
=
\left(
\begin{array}{c}
0 \\
2h\\
\vdots\\
\vdots \\
\vdots \\
\vdots \\
\vdots \\
2h +\frac{D}{h}\\
D
\end{array}
\right)
\end{equation*}
$$
Modyfikacja z zachowaniem symetrii zastosowana dla macierzy lokalnej \( \tilde A^{(N_e)} \):
$$
\begin{equation*}
\tilde A^{(N_e)} =
A = \frac{1}{h}\left(\begin{array}{rr}
1 & 0\\
0 & h
\end{array}\right),\quad
\tilde b^{(N_e)} = \left(\begin{array}{c}
h + D/h\\
D
\end{array}\right)
\end{equation*}
$$
Jak uwzględnić \( u'(0)=C \) w \( \FEM \)?
$$ -u''=f,\quad u'(0)=C,\ u(L)=D$$
Metoda Galerkina:
$$
\begin{equation*}
\int_0^L(u''(x)+f(x))\baspsi_i(x) dx = 0,\quad i\in\If
\end{equation*}
$$
Całkowanie \( u''\baspsi_i \) przez części:
$$
\begin{equation*}
\int_0^Lu'(x)\baspsi_i'(x) \dx -(u'(L)\baspsi_i(L) - u'(0)\baspsi_i(0)) -
\int_0^L f(x)\baspsi_i(x) \dx =0
\end{equation*}
$$
$$
\begin{equation*}
\int_0^Lu'(x)\basphi_i'(x) dx =
\int_0^L f(x)\basphi_i(x) dx - C\basphi_i(0),\quad i\in\If
\end{equation*}
$$
$$
\begin{equation*}
\sum_{j=0}^{N}\left(
\int_0^L \basphi_i'\basphi_j' dx \right)c_j =
\int_0^L\left(f\basphi_i -D\basphi_N'\basphi_i\right) dx
- C\basphi_i(0)
\end{equation*}
$$
for \( i=0,\ldots,N=N_n-2 \).
Można zatem pominąć człon \( u'(L)\basphi_i(L) \)!
Składowe \( u'\basphi_i \) w węzłach \( \xno{i} \), w których wymuszamy warunek Dirichleta, mogą zostać pominięte.
\footnotesize \vspace{-0.3cm}
$$
\begin{equation*}
u(x) = \sum_{j=0}^{N=N_n-1} c_j\basphi_j(x)
\end{equation*}
$$
$$
\begin{equation*}
\sum_{j=0}^{N=N_n-1}\left(
\int_0^L \basphi_i'(x)\basphi_j'(x) dx \right)c_j =
\int_0^L f(x)\basphi_i(x) dx - C\basphi_i(0)
\end{equation*}
$$
\normalsize
Trzeba obliczyć wszystkie elementy dla \( i=0,\ldots,N=N_n-1 \), a następnie zmodyfikować ostatnie równanie do \( c_N=D \).
Dodatkowy człon \( C\basphi_0(0) \) ma wpływ jedynie na wektor prawej strony pierwszego elementu (\( \basphi_0=0 \) na wszystkich pozostałych elementach).
$$
\begin{equation*}
\tilde A^{(0)} =
A = \frac{1}{h}\left(\begin{array}{rr}
1 & 1\\
-1 & 1
\end{array}\right),\quad
\tilde b^{(0)} = \left(\begin{array}{c}
h - C\\
h
\end{array}\right)
\end{equation*}
$$
Równanie różniczkowe definiuje całki sformułowania wariacyjnego.
Funkcje, które użytkownik musi podać na wejście programu:
integrand_lhs(phi, r, s, x)
boundary_lhs(phi, r, s, x)
integrand_rhs(phi, r, x)
boundary_rhs(phi, r, x)
(Ponadto potrzebna jest również informacja dotycząca siatki
zagadnienia, opisana przy pomocy struktur:
vertices, cells, oraz dof_map.
<Declare global matrix, global rhs: A, b>
# Loop over all cells
for e in range(len(cells)):
# Compute element matrix and vector
n = len(dof_map[e]) # no of dofs in this element
h = vertices[cells[e][1]] - vertices[cells[e][0]]
<Declare element matrix, element vector: A_e, b_e>
# Integrate over the reference cell
points, weights = <numerical integration rule>
for X, w in zip(points, weights):
phi = <basis functions + derivatives at X>
detJ = h/2
x = <affine mapping from X>
for r in range(n):
for s in range(n):
A_e[r,s] += integrand_lhs(phi, r, s, x)*detJ*w
b_e[r] += integrand_rhs(phi, r, x)*detJ*w
# Add boundary terms
for r in range(n):
for s in range(n):
A_e[r,s] += boundary_lhs(phi, r, s, x)*detJ*w
b_e[r] += boundary_rhs(phi, r, x)*detJ*w
for e in range(len(cells)):
...
# Incorporate essential boundary conditions
for r in range(n):
global_dof = dof_map[e][r]
if global_dof in essbc_dofs:
# dof r is subject to an essential condition
value = essbc_docs[global_dof]
# Symmetric modification
b_e -= value*A_e[:,r]
A_e[r,:] = 0
A_e[:,r] = 0
A_e[r,r] = 1
b_e[r] = value
# Assemble
for r in range(n):
for s in range(n):
A[dof_map[e][r], dof_map[e][r]] += A_e[r,s]
b[dof_map[e][r] += b_e[r]
<solve linear system>
$$
\begin{equation*}
-\int_{\Omega} \nabla\cdot (\dfc(\x)\nabla u) v\dx =
\int_{\Omega} \dfc(\x)\nabla u\cdot\nabla v \dx -
\int_{\partial\Omega} a\frac{\partial u}{\partial n} v \ds
\end{equation*}
$$
$$
\begin{align*}
\v\cdot\nabla u + \beta u &= \nabla\cdot\left( \dfc\nabla u\right) + f,
\quad & \x\in\Omega\\
u &= u_0,\quad &\x\in\partial\Omega_D\\
-\dfc\frac{\partial u}{\partial n} &= g,\quad &\x\in\partial\Omega_N
\end{align*}
$$
Metoda 1 (funkcja brzegowa i \( \baspsi_i=0 \) na \( \partial\Omega_D \)) gwarantuje spełnienie warunku \( u=u_0 \):
$$ u(\x) = B(\x) + \sum_{j\in\If} c_j\baspsi_j(\x),\quad B(\x)=u_0(\x) $$
Metoda Galerkina: pomnóż przez \( v\in V \) i scałkuj nad \( \Omega \),
$$
\int_{\Omega} (\v\cdot\nabla u + \beta u)v\dx =
\int_{\Omega} \nabla\cdot\left( \dfc\nabla u\right)v\dx + \int_{\Omega}fv \dx
$$
Całkowanie przez części pierwszej całki po prawej stronie zgodnie ze wzorem:
$$
\int_{\Omega} \nabla\cdot\left( \dfc\nabla u\right) v \dx =
-\int_{\Omega} \dfc\nabla u\cdot\nabla v\dx
+ \int_{\partial\Omega} \dfc\frac{\partial u}{\partial n} v\ds,
$$
daje ostatecznie
$$
\int_{\Omega} (\v\cdot\nabla u + \beta u)v\dx =
-\int_{\Omega} \dfc\nabla u\cdot\nabla v\dx
+ \int_{\partial\Omega} \dfc\frac{\partial u}{\partial n} v\ds
+ \int_{\Omega} fv \dx
$$
Uwaga: \( v\neq 0 \) jedynie na \( \partial\Omega_N \) (ponieważ \( v=0 \) na \( \partial\Omega_D \)):
$$ \int_{\partial\Omega} \dfc\frac{\partial u}{\partial n} v\ds
= \int_{\partial\Omega_N} \underbrace{\dfc\frac{\partial u}{\partial n}}_{-g} v\ds
= -\int_{\partial\Omega_N} gv\ds
$$
Ostateczna postać sformułowania wariacyjnego:
$$
\int_{\Omega} (\v\cdot\nabla u + \beta u)v\dx =
-\int_{\Omega} \dfc\nabla u\cdot\nabla v \dx
- \int_{\partial\Omega_N} g v\ds
+ \int_{\Omega} fv \dx
$$
Równoważnie w notacji abstrakcyjnej:
$$
(\v\cdot\nabla u, v) + (\beta u,v) =
- (\dfc\nabla u,\nabla v) - (g,v)_{N} + (f,v)
$$
\( (g,v)_{N} \): całka krzywoliniowa (2D) lub powierzchniowa (3D) nad \( \partial\Omega_N \).
$$
A_{i,j} = (\v\cdot\nabla \baspsi_j, \baspsi_i) +
(\beta \baspsi_j ,\baspsi_i) + (\dfc\nabla \baspsi_j,\nabla \baspsi_i)
$$
$$
b_i = (g,\baspsi_i)_{N} + (f,\baspsi_i) -
(\v\cdot\nabla u_0, \baspsi_i) + (\beta u_0 ,\baspsi_i) +
(\dfc\nabla u_0,\nabla \baspsi_i)
$$
Chcemy obliczyć całkę nad elementem całkując we współrzędnych unormowanych.
$$
\begin{equation*}
\int_{{\Omega}^{(e)}} \dfc(\x)\nabla\basphi_i\cdot\nabla\basphi_j\dx
\end{equation*}
$$
Odwzorowanie geometrii ze współrzędnych lokalnych (unormowanych) do współrzędnych globalnych:
$$ \x(\X) $$
gdzie Jakobian przekształcenia \( J \),
$$ J_{i,j}=\frac{\partial x_j}{\partial X_i} $$
Znajdujemy:
$$
\begin{align*}
\nabla_{\X}\refphi_r &= J\cdot\nabla_{\x}\basphi_i\\
\nabla_{\x}\basphi_i &= \nabla_{\x}\refphi_r(\X)
= J^{-1}\cdot\nabla_{\X}\refphi_r(\X)
\end{align*}
$$
Transformacja całki ze współrzędnych globalnych do lokalnych:
$$
\begin{equation*}
\int_{\Omega^{(e)}} \dfc(\x)\nabla_{\x}\basphi_i\cdot\nabla_{\x}\basphi_j\dx =
\int_{\tilde\Omega^r} \dfc(\x(\X))(J^{-1}\cdot\nabla_{\X}\refphi_r)\cdot
(J^{-1}\cdot\nabla\refphi_s)\det J\dX
\end{equation*}
$$
Całkowanie numeryczne na unormowanym elemencie (trójkącie lub czworościanie)
$$ \int_{\tilde\Omega^r} g\dX = \sum_{j=0}^{n-1} w_j g(\bar\X_j)$$
Moduł numint.py zawiera różne kwadratury całkujące:
>>> import numint
>>> x, w = numint.quadrature_for_triangles(num_points=3)
>>> x
[(0.16666666666666666, 0.16666666666666666),
(0.66666666666666666, 0.16666666666666666),
(0.16666666666666666, 0.66666666666666666)]
>>> w
[0.16666666666666666, 0.16666666666666666, 0.16666666666666666]
/