Wstęp

Ta prezentacja

• PDF
• HTML (jasny)
• reveal (jasny)
• reveal (ciemny)
• deck.js (jasny)

Kody Pythona

Hans Petter Langtangen (1962-2016)

• Strona domowa
• Github
• DocOnce i książki HPL
• Książka o FEM (EN)

Metoda elementów skończonych - wstęp

Finite Element Method, FEM, MES

Delfin

Rozwiązywanie \( \PDE \) przy pomocy FEM

Etapy nauki FEM

Aproksymacja w przestrzeniach wektorowych

Jak znaleźć wektor pewnej przestrzeni, który aproksymuje wektor przestrzeni o większym wymiarze?

Aproksymacja jako kombinacja liniowa założonych funkcji bazowych

Ogólna idea poszukiwania elementu (wektora/funkcji) \( u(x) \) pewnej przestrzeni przybliżającego zadany element \( f(x) \):

 
$$ u(x) = \sum_{i=0}^N c_I \baspsi_i(x) $$

gdzie

• \( \baspsi_i(x) \) założone funkcje
• \( c_i \), \( i=0,\ldots,N \) nieznane współczynniki (do wyznaczenia)

Trzy główne sposoby wyznaczania niewiadomych współczynników

• metoda najmniejszych kwadratów (ang. Least Squares Method LSM)
• metoda Galerkina
• metoda kolokacji

Aproksymacja wektorów: przykład – aproksymacja na płaszczyznie

Zadanie:

Znajdź przybliżenie wektora \( \f = (3,5) \) wzdłuż zadanego kierunku.

analogie: element - punkt - wektor - funkcja

Wektory 2D - rzut ortogonalny

Aproksymacja wektorów: przestrzenie wektorowe – terminologia

 
$$ V = \mbox{span}\,\{ \psib_0\} $$

• \( \psib_0 \) wektor bazowy przestrzeni \( V \)
• Znajdź \( \u = c_0\psib_0\in V \)
• Jak wyznaczyć \( c_0 \) tak, aby \( \u \) "najlepiej" przybliżało \( \f \)?
• Wizualnie rozwiązanie narzuca się samo
• Jak sformułować to ogólnie, w języku matematyki?

Niech

• \( \e = \f - \u \) to błąd
• \( (\u,\v) \) – iloczyn skalarny wektorów (-> sens geometryczny il.s. -> przykład?)
• \( ||\e|| = \sqrt{(\e, \e)} \) – norma błędu (jaka? (normy \( p \)-te))

Metoda najmniejszych kwadratów - idea

• Idea: znaleźć \( c_0 \) takie, aby \( ||\e|| \) był minimalizowany (jak najmniejszy/najkrótszy)
• Dla wygody (matematycznej): minimalizujemy \( E=||\e||^2 \)

 
$$ \begin{equation*} \frac{\partial E}{\partial c_0} = 0 \end{equation*} $$

Metoda najmniejszych kwadratów - obliczenia

 
$$ \begin{align*} E(c_0) &= (\e,\e) = (\f - \u, \f - \u) = (\f - c_0\psib_0, \f - c_0\psib_0)\\ &= (\f,\f) - 2c_0(\f,\psib_0) + c_0^2(\psib_0,\psib_0) \end{align*} $$

 
$$ \begin{equation} \frac{\partial E}{\partial c_0} = -2(\f,\psib_0) + 2c_0 (\psib_0,\psib_0) = 0 \tag{1} \end{equation} $$

 
$$ c_0 = \frac{(\f,\psib_0)}{(\psib_0,\psib_0)} = \frac{3a + 5b}{a^2 + b^2} $$

Spostrzeżenie (na później): warunek znikania pochodnej (1) można równoważnie zapisać jako:

 
$$ (\e, \psib_0) = 0 $$

 
$$ (\e, \psib_0) = (\f - \u, \psib_0) = (\f, \psib_0) - (\u, \psib_0) .... \quad\quad (*-2 ...) $$

 
$$ [\f-\u] \cdot \left[ \begin{array}{c} \\ \psib_0 \\ \\ \end{array} \right] = [\quad\quad \f \quad\quad] \cdot \left[ \begin{array}{c} \\ \psib_0 \\ \\ \end{array} \right] - [\quad\quad \u \quad\quad] \cdot \left[ \begin{array}{c} \\ \psib_0 \\ \\ \end{array} \right] $$

Metoda rzutu ortogonalnego (m. Galerkina)

Aproksymacja wektora przestrzeni dowolniewymiarowej

Dla danego wektora \( \f \), znajdź przybliżenie \( \u\in V \):

 
$$ \begin{equation*} V = \hbox{span}\,\{\psib_0,\ldots,\psib_N\} \end{equation*} $$

(\( \hbox{span} \) czyt. przestrzeń rozpięta na wektorach...)

Mając dany zbiór wektorów niezależnych liniowo \( \psib_0,\ldots,\psib_N \), dowolny wektor \( \u\in V \) można zapisać jako:

 
$$ \u = \sum_{j=0}^Nc_j\psib_j$$

Metoda najmniejszych kwadratów

Idea: znaleźć takie \( c_0,\ldots,c_N \), aby \( E= ||\e||^2 \) był minimalizowany, \( \e=\f-\u \).

 
$$ \begin{align*} E(c_0,\ldots,c_N) &= (\e,\e) = (\f -\sum_jc_j\psib_j,\f -\sum_jc_j\psib_j) \nonumber\\ &= (\f,\f) - 2\sum_{j=0}^Nc_j(\f,\psib_j) + \sum_{p=0}^N\sum_{q=0}^N c_pc_q(\psib_p,\psib_q) \end{align*} $$

 
$$ \begin{equation*} \frac{\partial E}{\partial c_i} = 0,\quad i=0,\ldots,N \end{equation*} $$

Po odrobinie obliczeń otrzymuje się układ równań liniowych:

 
$$ \begin{align} \sum_{j=0}^N A_{i,j}c_j &= b_i,\quad i=0,\ldots,N \tag{2}\\ A_{i,j} &= (\psib_i,\psib_j) \tag{3}\\ b_i &= (\psib_i, \f) \tag{4} \end{align} $$

Metoda Galerkina

Można pokazać, że poszukiwanie minimalnego \( ||\e|| \) jest równoważne poszukiwaniu \( \e \) ortogonalnego do wszystkich \( \v\in V \):

 
$$ (\e,\v)=0,\quad \forall\v\in V $$

co jest równoważne temu aby \( \e \) był ortogonalny do każdego wektora bazowego:

 
$$ (\e,\psib_i)=0,\quad i=0,\ldots,N $$

Warunek ortogonalności – podstawa metody Galerkina. Generuje ten sam układ równań liniowych co \( \LSM \).

Aproksymacja funkcji w przestrzeni funkcyjnej

Niech \( V \) będzie przestrzenią funkcyjną rozpiętą na zbiorze funkcji bazowych \( \baspsi_0,\ldots,\baspsi_N \),

 
$$ \begin{equation*} V = \hbox{span}\,\{\baspsi_0,\ldots,\baspsi_N\} \end{equation*} $$

Dowolną funkcję tej przestrzeni \( u\in V \) można przedstawić jako kombinację liniową funkcji bazowych:

 
$$ u = \sum_{j\in\If} c_j\baspsi_j,\quad\If = \{0,1,\ldots,N\} $$

Uogólnienie \( \LSM \) na przestrzenie funkcyjne

Tak jak dla przestrzeni wektorowych, minimalizujemy normę błędu \( E \), ze względu na współczynniki \( c_j \), \( j\in\If \):

 
$$ E = (e,e) = (f-u,f-u) = \left(f(x)-\sum_{j\in\If} c_j\baspsi_j(x), f(x)-\sum_{j\in\If} c_j\baspsi_j(x)\right) $$

 
$$ \frac{\partial E}{\partial c_i} = 0,\quad i=\in\If $$

Czym jest iloczyn skalarny jeśli \( \baspsi_i \) jest funkcją?

 
$$(f,g) = \int_\Omega f(x)g(x)\, dx$$

(iloczyn skalarny funkcji cg jako uogólnienie il. skalarnego funkcji dyskretnych \( (\u, \v) = \sum_j u_jv_j \))

Szczegóły \( \LSM \)

 
$$ \begin{align*} E(c_0,\ldots,c_N) &= (e,e) = (f-u,f-u) \\ &= (f,f) -2\sum_{j\in\If} c_j(f,\baspsi_i) + \sum_{p\in\If}\sum_{q\in\If} c_pc_q(\baspsi_p,\baspsi_q) \end{align*} $$

 
$$ \frac{\partial E}{\partial c_i} = 0,\quad i=\in\If $$

Obliczenia identyczne jak dla przypadku wektorowego -> w rezultacie otrzymujemy układ równań liniowych

 
$$ \sum_{j\in\If}^N A_{i,j}c_j = b_i,\ i\in\If,\quad A_{i,j} = (\baspsi_i,\baspsi_j),\ b_i = (f,\baspsi_i) $$

Metoda Galerkina

Jak poprzednio minimalizacja \( (e,e) \) jest równoważna

 
$$ (e,\baspsi_i)=0,\quad i\in\If \tag{5} $$

co z kolei jest równoważne

 
$$ (e,v)=0,\quad\forall v\in V \tag{6} $$

Równoważność wzorów jak dla przestrzeni wektorowych.

Równoważność wzorów jak dla wyprowadzenia przy pomocy \( \LSM \).

Przykład: aproksymacja paraboli funkcją liniową

 
$$ \begin{equation*} V = \hbox{span}\,\{1, x\} \end{equation*} $$

czyli \( \baspsi_0(x)=1 \), \( \baspsi_1(x)=x \) oraz \( N=1 \). Szukane

 
$$ \begin{equation*} u=c_0\baspsi_0(x) + c_1\baspsi_1(x) = c_0 + c_1x \end{equation*} $$

Przykład: aproksymacja paraboli funkcją liniową

 
$$ \begin{align*} A_{0,0} &= (\baspsi_0,\baspsi_0) = \int_1^21\cdot 1\, dx = 1\\ A_{0,1} &= (\baspsi_0,\baspsi_1) = \int_1^2 1\cdot x\, dx = 3/2\\ A_{1,0} &= A_{0,1} = 3/2\\ A_{1,1} &= (\baspsi_1,\baspsi_1) = \int_1^2 x\cdot x\,dx = 7/3\\ b_1 &= (f,\baspsi_0) = \int_1^2 (10(x-1)^2 - 1)\cdot 1 \, dx = 7/3\\ b_2 &= (f,\baspsi_1) = \int_1^2 (10(x-1)^2 - 1)\cdot x\, dx = 13/3 \end{align*} $$

Rozwiązanie układu równań 2x2:

 
$$ c_0 = -38/3,\quad c_1 = 10,\quad u(x) = 10x - \frac{38}{3} $$

Przykład: aproksymacja paraboli funkcją liniową

 
$$ \begin{align*} \left[ \begin{array}{cc} 1 & \frac{3}{2} \\ \frac{3}{2} & \frac{7}{3} \\ \end{array} \right] \left[ \begin{array}{c} c_0 \\ c_1 \\ \end{array} \right] = \left[ \begin{array}{c} \frac{7}{3} \\ \frac{13}{3} \\ \end{array} \right] \quad \rightarrow \quad \left[ \begin{array}{c} c_0 \\ c_1 \\ \end{array} \right] = \left[ \begin{array}{c} -38/3 \\ 10 \\ \end{array} \right] \end{align*} $$

 
$$ u(x) = 10x - 12 \frac{2}{3} $$

Symboliczna realizacja algorytmu \( \LSM \)

Problem: napisać program/funkcję, który przeprowadzi obliczenia (obliczenie całek i rozwiązanie układu równań liniowych) i zwróci rozwiązanie postaci n \( u(x)=\sum_jc_j\baspsi_j(x) \).

Niech

• \( f(x) \) będzie dane przez funkcję sympy oznaczoną symbolem f (funkcję zmiennej (symbolu) x)
• psi będzie listą funkcji \( \sequencei{\baspsi} \),
• Omega będzie dwuelementową krotką/listą zawierającą początek i koniec przedziału \( \Omega \)

\( \LSM \) symbolicznie: podejście nr 1

import sympy as sym

def least_squares(f, psi, Omega):
    N = len(psi) - 1
    A = sym.zeros((N+1, N+1))
    b = sym.zeros((N+1, 1))
    x = sym.Symbol('x')
    for i in range(N+1):
        for j in range(i, N+1):
            A[i,j] = sym.integrate(psi[i]*psi[j],
                                  (x, Omega[0], Omega[1]))
            A[j,i] = A[i,j]
        b[i,0] = sym.integrate(psi[i]*f, (x, Omega[0], Omega[1]))
    c = A.LUsolve(b)
    u = 0
    for i in range(len(psi)):
        u += c[i,0]*psi[i]
    return u, c

Spostrzeżenie: macierz układu jest symetryczna, dzięki czemu można zoptymalizować proces wyznaczania jej elementów

• Może się zdażyć, że obliczanie całki się nie powiedzie (skomplikowana funkcja f) , sym.integrate zwróci wtedy obiekt typu sym.Integral. Ulepszenie kodu: sprawdzenie czy takie zdarzenie wystąpiło i ew. obliczenia numeryczne.
• Ulepszenie 2: Dodatkowa flaga przy pomocy, której użytkownik może zdecydować bezpośrednio jaki rodzaj całkowania (symboliczny czy numeryczny) ma zostać wykorzystany.

\( \LSM \) symbolicznie: podejście nr 2

def least_squares(f, psi, Omega, symbolic=True):
    ...
    for i in range(N+1):
        for j in range(i, N+1):
            integrand = psi[i]*psi[j]
            if symbolic:
                I = sym.integrate(integrand, (x, Omega[0], Omega[1]))
            if not symbolic or isinstance(I, sym.Integral):
                # Could not integrate symbolically,
                # fall back on numerical integration
                integrand = sym.lambdify([x], integrand)
                I = sym.mpmath.quad(integrand, [Omega[0], Omega[1]])
            A[i,j] = A[j,i] = I

        integrand = psi[i]*f
        if symbolic:
            I = sym.integrate(integrand, (x, Omega[0], Omega[1]))
        if not symbolic or isinstance(I, sym.Integral):
            integrand = sym.lambdify([x], integrand)
            I = sym.mpmath.quad(integrand, [Omega[0], Omega[1]])
        b[i,0] = I
    ...

Prezentacja rozwiązania

Graficzne porównanie \( f \) i \( u \):

def comparison_plot(f, u, Omega, filename='tmp.pdf'):
    x = sym.Symbol('x')
    # Turn f and u to ordinary Python functions
    f = sym.lambdify([x], f, modules="numpy")
    u = sym.lambdify([x], u, modules="numpy")
    resolution = 401  # no of points in plot
    xcoor  = linspace(Omega[0], Omega[1], resolution)
    exact  = f(xcoor)
    approx = u(xcoor)
    plot(xcoor, approx)
    hold('on')
    plot(xcoor, exact)
    legend(['approximation', 'exact'])
    savefig(filename)

Zastosowanie kodu

>>> from approx1D import *
>>> x = sym.Symbol('x')
>>> f = 10*(x-1)**2-1
>>> u, c = least_squares(f=f, psi=[1, x], Omega=[1, 2])
>>> comparison_plot(f, u, Omega=[1, 2])

Przypadek aproksymacji funkcji \( f\in V \)

• Rozszerzmy zbiór funkcji bazowych przestrzeni \( V \) o funkcję \( \baspsi_2=x^2 \), wciąż poszukując przybliżenia dla funkcji \( f=10(x-1)^2-1 \) w przestrzeni \( V \)
• -> przybliżenie paraboli pewną parabolą \ldots
• Rozwiązanie odwzoruje \( f \) ściśle!

>>> from approx1D import *
>>> x = sym.Symbol('x')
>>> f = 10*(x-1)**2-1
>>> u, c = least_squares(f=f, psi=[1, x, x**2], Omega=[1, 2])
>>> print u
10*x**2 - 20*x + 9
>>> print sym.expand(f)
10*x**2 - 20*x + 9

Rozwiązanie przybliżone \( \equiv \) rozwiązanie dokładne, jeśli \( f \in V \)!

Uogólnienie: Przypadek aproksymacji funkcji \( f\in V \)

• A co jeśli baza to \( \psi_i(x)=x^i \) dla \( i=0,\ldots,N=40 \)?
• Wynik funkcji least_squares: dla \( i>2 \), \( c_i=0 \)

Dlaczego dla \( f\in V \) aproksymacja jest bezbłędna? Dowód:

Jeśli \( f\in V \), wtedy \( f=\sum_{j\in\If}d_j\baspsi_j \), dla pewnego \( \sequencei{d} \). Wtedy

 
$$ \begin{equation*} b_i = (f,\baspsi_i) = \sum_{j\in\If}d_j(\baspsi_j, \baspsi_i) = \sum_{j\in\If} d_jA_{i,j} \end{equation*} $$

a URL \( \sum_j A_{i,j}c_j = b_i \), \( i\in\If \), przedstawia sie:

 
$$ \begin{equation*} \sum_{j\in\If}c_jA_{i,j} = \sum_{j\in\If}d_jA_{i,j},\quad i\in\If \end{equation*} $$

co oznacza, że \( c_i=d_i \) dla \( i\in\If \), czyli \( u \) jest tożsame z \( f \).

Skończona precyzja obliczeń numerycznych

Poprzednie wnioski -> teoria i obliczenia symboliczne \ldots

Co w przypadku obliczeń numerycznych? -> (rozwiązanie URL macierzami liczb zmiennoprzecinkowych)

\( f \) to wciąż funkcja kwadratowa przybliżana przez

 
$$ u(x) = c_0 + c_1x + c_2x^2 + c_3x^3 +\cdots + c_Nx^N $$

Oczekiwane: \( c_2=c_3=\cdots=c_N=0 \), skoro \( f\in V \) oznacza \( u=f \).

A naprawdę?

Skończona prezycja obliczeń numerycznych – wyniki

• Kolumna 2: matrix oraz lu_solve z biblioteki sympy.mpmath.fp
• Kolumna 3: numpy 4B liczby zmiennoprzecinkowe
• Kolumna 4: numpy 8B liczby zmiennoprzecinkowe

Złe uwarunkowanie URL - ''liniowa zależność'' w bazie

• Znaczne błędy zaokrągleń rozwiązania numerycznego (!)
• Jednocześnie ''na oko'' (graficzne) rozwiązanie wygląda w porządku (!)

Źródło kłopotów: funkcje \( x^i \) dla bardzo dużych \( i \) stają się praktycznie liniowo zależne

Złe uwarunkowanie URL: wnioski

• Prawie liniowa zależność funkcji bazowych skutkuje bliskoosobliwymi macierzami
• macierz prawie osobliwa \( \equiv \) macierz źle uwarunkowana -> problemy w trakcie m.elim. Gaussa
• Baza wielomianów \( 1, x, x^2, x^3, x^4, \ldots \) to ''nienajszczęśliwszy'' wybór
• Istnieją lepsze bazy (nawet wielomianowe), ale im bardziej ortogonalne te bazy są, tym lepiej (\( (\baspsi_i,\baspsi_j)\approx 0 \))

Aproksymacja szeregami Fouriera; problem and code

Aproksymacja funkcji \( f \) szeregiem Fouriera

 
$$ u(x) = \sum_i a_i\sin i\pi x = \sum_{j=0}^Nc_j\sin((j+1)\pi x) $$

to tylko ''zmiana bazy'':

 
$$ \begin{equation*} V = \hbox{span}\,\{ \sin \pi x, \sin 2\pi x,\ldots,\sin (N+1)\pi x\} \end{equation*} $$

Obliczenia z wykorzystaniem funkcji least_squares:

N = 3
from sympy import sin, pi
psi = [sin(pi*(i+1)*x) for i in range(N+1)]
f = 10*(x-1)**2 - 1
Omega = [0, 1]
u, c = least_squares(f, psi, Omega)
comparison_plot(f, u, Omega)

Aproksymacja szeregami Fouriera; wykres

L: \( N=3 \), P: \( N=11 \):

Aproksymacja szeregami Fouriera; ulepszenie

• Znaczna poprawa aproksymacja dla \( N=11 \) wyrazów, pomimo niepożądanych rozbieżności w \( x=0 \) i \( x=1 \)
• Możliwe rozwiązanie: dodać składową, która pozwoli na odwzorowanie właściwych wartości na brzegu

 
$$ u(x) = {\color{red}f(0)(1-x) + xf(1)} + \sum_{j\in\If} c_j\baspsi_j(x) $$

Dodatkowy wyraz nie tylko zapewnia \( u(0)=f(0) \) oraz \( u(1)=f(1) \), ale także zaskakująco dobrze poprawia jakość aproksymacji!

Aproksymacja szeregami Fouriera; wyniki

\( N=3 \) vs \( N=11 \):

Bazy funkcji ortogonalnych

Zalety wyboru funkcji sinus jako funkcji bazowych:

• funkcje bazowe są parami ortogonalne: \( (\baspsi_i,\baspsi_j)=0 \) (jedynie \( (\baspsi_i,\baspsi_j) \neq 0 \))) dzięki czemu
• macierz \( A_{i,j} \) jest diagonalna, dzięki czemu
• nie ma potrzeby rozwiązywać URL! Rozwiązanie sprowadza się do obliczenia: \( c_i = 2\int_0^1 f(x)\sin ((i+1)\pi x) dx \)
• wynik: rozwinięcie funkcji \( f \) w szereg Fouriera

W ogólnym przypadku, dla baz ortogonalnych, \( A_{i,j} \) jest macierzą diagonalną, a nieznane współczynniki \( c_i \) można łatwo obliczyć:

 
$$ c_i = \frac{b_i}{A_{i,i}} = \frac{(f,\baspsi_i)}{(\baspsi_i,\baspsi_i)} $$

Implementacja \( \LSM \) dla ortogonalnych funkcji bazowych

def least_squares_orth(f, psi, Omega):
    N = len(psi) - 1
    A = [0]*(N+1)
    b = [0]*(N+1)
    x = sym.Symbol('x')
    for i in range(N+1):
        A[i] = sym.integrate(psi[i]**2, (x, Omega[0], Omega[1]))
        b[i] = sym.integrate(psi[i]*f,  (x, Omega[0], Omega[1]))
    c = [b[i]/A[i] for i in range(len(b))]
    u = 0
    for i in range(len(psi)):
        u += c[i]*psi[i]
    return u, c

Implementacja \( \LSM \) dla ortogonalnych funkcji bazowych: całkowanie symboliczne i numeryczne

• Uwzględnienie parametru sterującego wyborem rodzaju całkowania (argument symbolic).
• W przypadku gdy całkowanie symboliczne zawiedzie (sym.integrate zwraca sym.Integral), obliczenia wykonywane numerycznie (w przypadku funkcji sinus nie powinno być problemów z symbolicznym obliczeniem \( \int_\Omega\basphi_i^2dx \))

def least_squares_orth(f, psi, Omega, symbolic=True):
    ...
    for i in range(N+1):
        # Diagonal matrix term
        A[i] = sym.integrate(psi[i]**2, (x, Omega[0], Omega[1]))

        # Right-hand side term
        integrand = psi[i]*f
        if symbolic:
            I = sym.integrate(integrand,  (x, Omega[0], Omega[1]))
        if not symbolic or isinstance(I, sym.Integral):
            print 'numerical integration of', integrand
            integrand = sym.lambdify([x], integrand)
            I = sym.mpmath.quad(integrand, [Omega[0], Omega[1]])
        b[i] = I
    ...

Metoda kolokacji (interpolacji); idea i teoria

Inny sposób znalezienia przybliżenia \( f(x) \) przez \( u(x)=\sum_jc_j\baspsi_j \):

• Wymuszamy \( u(\xno{i}) = f(\xno{i}) \) w pewnych wybranych punktach \( \sequencei{x} \) (punktach kolokacji)
• \( u \) interpoluje \( f \)
• metoda znana jako metoda kolokacji (interpolacji)

 
$$ u(\xno{i}) = \sum_{j\in\If} c_j \baspsi_j(\xno{i}) = f(\xno{i}) \quad i\in\If,N $$

Współczynniki wygenerowanego układu równań to po prostu wartości funkcji, nie ma potrzeby całkowania!

 
$$ \begin{align} \sum_{j\in\If} A_{i,j}c_j &= b_i,\quad i\in\If \tag{7}\\ A_{i,j} &= \baspsi_j(\xno{i}) \tag{8}\\ b_i &= f(\xno{i}) \tag{9} \end{align} $$

W ogólnym przypadku macierz wynikowa niesymetryczna: \( \baspsi_j(\xno{i})\neq \baspsi_i(\xno{j}) \)

Metoda kolokacji – implementacja

Zmienna points przechowuje punkty kolokacji

def interpolation(f, psi, points):
    N = len(psi) - 1
    A = sym.zeros((N+1, N+1))
    b = sym.zeros((N+1, 1))
    x = sym.Symbol('x')
    # Turn psi and f into Python functions
    psi = [sym.lambdify([x], psi[i]) for i in range(N+1)]
    f = sym.lambdify([x], f)
    for i in range(N+1):
        for j in range(N+1):
            A[i,j] = psi[j](points[i])
        b[i,0] = f(points[i])
    c = A.LUsolve(b)
    u = 0
    for i in range(len(psi)):
        u += c[i,0]*psi[i](x)
    return u

Metoda kolokacji: przybliżenie paraboli funkcją liniową

• Problem: jak wybrać \( \xno{i} \)?
• Wynik zależy od położenia punktów kolokacji!

\( (4/3,5/3) \) vs \( (1,2) \):

Regresja

• Idea: metoda kolokacji dla \( m \gg N+1 \) punktów
• Problem: Więcej równań niż niewiadomych
• Znana np. ze statystyki regresja

Regresja – nadokreślone URL

 
$$ u(\xno{i}) = \sum_{j\in\If} c_j \baspsi_j(\xno{i}) = f(\xno{i}), \quad i=0,1,\ldots,m $$

 
$$ \sum_{j\in\If} A_{i,j}c_j = b_i,\quad i=0,1,\ldots,m $$

 
$$ A_{i,j} = \baspsi_j(\xno{i}),\quad b_i = f(\xno{i})$$

Rozwiązywanie nadokreślonych URL przy pomocy \( \LSM \)

• Jak rozwiązać \( Ac=b \) jeśli jest więcej równań niż niewiadomych?
• Idea: Poszukiwanie rozwiązania minimalizującego \( r=b-Ac \)
• Rezultat: układ równań normalnych \( A^TAc=A^Tb \)
• Zapiszmy układ w postaci \( Bc=d \)
• \( B = A^TA \) już kwadratowe: układ równań o rozmiarach \( (N+1)\times(N+1) \)

 
$$ \begin{align*} B_{i,j} &= \sum_k A^T{i,k}A_{k,j} = \sum_k A{k,i}A_{k,j} =\sum_{k=0}^m\baspsi_i(\xno{k}) \baspsi_j(\xno{k}) \\ d_i &=\sum_k A^T_{i,k}b_k = \sum_k A_{k,i}b_k =\sum_{k=0}^m \baspsi_i(\xno{k})f(\xno{k}) \end{align*} $$

Implementacja

def regression(f, psi, points):
    N = len(psi) - 1
    m = len(points)
    # Use numpy arrays and numerical computing
    B = np.zeros((N+1, N+1))
    d = np.zeros(N+1)
    # Wrap psi and f in Python functions rather than expressions
    # so that we can evaluate psi at points[i]
    x = sym.Symbol('x')
    psi_sym = psi  # save symbolic expression
    psi = [sym.lambdify([x], psi[i]) for i in range(N+1)]
    f = sym.lambdify([x], f)
    for i in range(N+1):
        for j in range(N+1):
            B[i,j] = 0
            for k in range(m+1):
                B[i,j] += psi[i](points[k])*psi[j](points[k])
        d[i] = 0
        for k in range(m+1):
            d[i] += psi[i](points[k])*f(points[k])
    c = np.linalg.solve(B, d)
    u = sum(c[i]*psi_sym[i] for i in range(N+1))
    return u, c

Przykład zastosowania – kod

• Zadanie: Dokonać aproksymacji funkcji \( f(x)=10(x-1)^2-1 \) na przedziale \( \Omega=[1,2] \) przy pomocy funkcji liniowej.

import sympy as sym
x = sym.Symbol('x')
f = 10*(x-1)**2 - 1
psi = [1, x]
Omega = [1, 2]
m_values = [2-1, 8-1, 64-1]
# Create m+3 points and use the inner m+1 points
for m in m_values:
    points = np.linspace(Omega[0], Omega[1], m+3)[1:-1]
    u, c = regression(f, psi, points)
    comparison_plot(f, u, Omega, points=points,
        points_legend='%d interpolation points' % (m+1))

Przykład zastosowania – wyniki

 
$$ \begin{align*} u(x) &= 10x - 13.2,\quad 2\hbox{ punkty}\\ u(x) &= 10x - 12.7,\quad 8\hbox{ punktów}\\ u(x) &= 10x - 12.7,\quad 64\hbox{ punkty} \end{align*} $$

Wielomiany Lagrange'a

Motywacja::

• Metoda kolokacji pozwala uniknąć całkowania
• Dla macierzy diagonalnej \( A_{i,j} = \baspsi_j(\xno{i}) \) rozwiązanie URL jest banalnie proste

Własność wielomianów Lagrange'a \( \baspsi_j \):

 
$$ \baspsi_i(\xno{j}) =\delta_{ij},\quad \delta_{ij} = \left\lbrace\begin{array}{ll} 1, & i=j\\ 0, & i\neq j \end{array}\right. $$

Zatem, \( c_i = f(x_i) \) and

 
$$ u(x) = \sum_{j\in\If} f(\xno{i})\baspsi_i(x) $$

• Wielomiany Lagrange'a w połączeniu z metodą kolokacji są niezwykle wygodne
• Często stosowane w \( \FEM \)

Wielomiany Lagrange'a – wzór i implementacja

 
$$ \baspsi_i(x) = \prod_{j=0,j\neq i}^N \frac{x-\xno{j}}{\xno{i}-\xno{j}} = \frac{x-x_0}{\xno{i}-x_0}\cdots\frac{x-\xno{i-1}}{\xno{i}-\xno{i-1}}\frac{x-\xno{i+1}}{\xno{i}-\xno{i+1}} \cdots\frac{x-x_N}{\xno{i}-x_N} $$

def Lagrange_polynomial(x, i, points):
    p = 1
    for k in range(len(points)):
        if k != i:
            p *= (x - points[k])/(points[i] - points[k])
    return p

Wielomiany Lagrange'a – zachęcający przykład zastosowania

Wielomiany Lagrange'a – mniej zachęcający przykład zastosowania

Wielomiany Lagrange'a – efekt Runge'go

12 punktów, dwa wielomiany stopnia 11 (Uwaga!: \( \psi_2(x_2) \neq 0 \) i \( \psi_7(x_7) \neq 0 \)

Problem: oscylacje w okolicach krańców przedziałów dla większych \( N \).

Wielomiany Lagrange'a: jak zapobiec oscylacjom?

Odpowiedni dobór węzłów interpolacji – węzły Czebyszewa:

 
$$ \xno{i} = \half (a+b) + \half(b-a)\cos\left( \frac{2i+1}{2(N+1)}\pi\right),\quad i=0\ldots,N $$

na przedziale \( [a,b] \).

Wielomiany Lagrange'a + węzły Czebyszewa

Wielomiany Lagrange'a + węzły Czebyszewa

12 punktów, dwa wielomiany stopnia 11.

Uwaga!: Tym razem węzły są inaczej rozmieszczone!

Mniej oscylacyjny charakter wielomianów w porównaniu do węzłów równomiernie rozmieszczonych.

Funkcje bazowe elementów skończonych

Funkcje bazowe o nośniku nieograniczonym: \( \baspsi_i(x) \neq 0 \) prawie w całym przedziale określoności

Nośnik funkcji: domknięcie zbioru argumentów funkcji, dla których ma ona wartość różną od zera (takie iksy dla których \( f(x) \neq 0 \))

Funkcje bazowe o nośniku ograniczonym – \( \FEM \)

• Nośnik zwarty (Local support): domknięcie zbioru tych \( x \)-ów, dla których \( \baspsi_i(x) \neq 0 \)
• Typowe dla 1D - funkcje trójkątne (hat-shaped)
• \( u(x) \) zbudowana przy pomocy takich funkcji \( \baspsi_i \) będzie funkcją przedziałami liniową
• Niech symbol \( \basphi_i \) oznacza odtąd tego typu funkcję trójkątną (przyjmijmy również \( \baspsi_i=\basphi_i \))

Kombinacja liniowa funkcji trójkątnych

Elementy i węzły

Podzielmy \( \Omega \) na \( N_e \) rozdzielnych podobszarów podobszarów – elementów:

 
$$ \Omega = \Omega^{(0)}\cup \cdots \cup \Omega^{(N_e)} $$

Na każdym elemencie wprowadzamy \( N_n \) wezłów (punktów): \( \xno{0},\ldots,\xno{N_n-1} \)

• \( \basphi_i(x) \) – \( i \)-ta funkcja bazowa
• \( \basphi_i=1 \) w węźle \( i \) i \( \basphi_i=0 \) w pozostałych węzłach
• \( \basphi_i \) to wielomian Lagrange'a na każdym elemencie
• Dla węzłów granicznych, leżących w punktach łączących dwa elementy funkcja \( \basphi_i \) jest zbudowane z wielomianów Lagrange'a na obu elementach

Przykład: obszar podzielony na elementy dwuwęzłowe (elementy typu P1)

Struktura nodes – współrzędne węzłów.

Struktura elements – numery (globalne) węzłów tworzących odpowiedni element.

nodes = [0, 1.2, 2.4, 3.6, 4.8, 5]
elements = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5]]

Przykład: dwie funkcje bazowe na siatce

Przykład: elementy niejednorodne o trzech węzłach (elementy typu P2)

nodes = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
elements = [[0, 1, 2], [2, 3, 4], [4, 5, 6], [6, 7, 8]]

Przykład: funkcje bazowe na siatce (elementy typu P2)

Przykład: elementy typu P3 (o czterech węzłach interpolacji)

d = 3  # d+1 nodes per element
num_elements = 4
num_nodes = num_elements*d + 1
nodes = [i*0.5 for i in range(num_nodes)]
elements = [[i*d+j for j in range(d+1)] for i in range(num_elements)]

Przykład: funkcje bazowe na siatce (elementy typu P3)

Indeksacja nieregularna

nodes = [1.5, 5.5, 4.2, 0.3, 2.2, 3.1]
elements = [[2, 1], [4, 5], [0, 4], [3, 0], [5, 2]]

Współczynniki \( c_i \) – interpretacja

Ważna własność: \( c_i \) to wartość funkcji \( u \) w węźle \( i \), \( \xno{i} \):

 
$$ u(\xno{i}) = \sum_{j\in\If} c_j\basphi_j(\xno{i}) = c_i\basphi_i(\xno{i}) = c_i $$

Powód: \( \basphi_j(\xno{i}) =0 \) jeśli \( i\neq j \) i \( \basphi_i(\xno{i}) =1 \)

Własności funkcji bazowych

• \( \basphi_i(x) \neq 0 \) jedynie na tych elementach, które zawierają węzeł

  1. globalnym indeksie \( i \)

• \( \basphi_i(x)\basphi_j(x) \neq 0 \) wtedy i tylko wtedy gdy węzly \( i \) oraz \( j \) leżą na tym samym elemencie

Ponieważ \( A_{i,j}=\int\basphi_i\basphi_j\dx \), większość współczynników macierzy będzie równa zero -> macierze rzadkie

Konstrukcja kwadratowych \( \basphi_i \) (elementy typu P2)

1. każdy węzeł elementu ma przypisany wielomian Lagrange'a
2. Wielomian o wartości 1 na brzegu elementu należy ''połączyć'' z wielomianem z sąsiedniego elementu, który ma wartość 1 w tym samym punkcie

Liniowe \( \basphi_i \) (elementy typu P1)

 
$$ \basphi_i(x) = \left\lbrace\begin{array}{ll} 0, & x < \xno{i-1}\\ (x - \xno{i-1})/h & \xno{i-1} \leq x < \xno{i}\\ 1 - (x - x_{i})/h, & \xno{i} \leq x < \xno{i+1}\\ 0, & x\geq \xno{i+1} \end{array} \right. $$

Sześcienne \( \basphi_i \) (elementy typu P3)

Generowanie URL

Przykład 1: Obliczenie wartości (niediagonalnego) elementu macierzy

Uproszczenie: elementy jednakowej długości.

\( A_{2,3}=\int_\Omega\basphi_2\basphi_3 dx \): \( \basphi_2\basphi_3\neq 0 \) jedynie na elemencie 2. Dla tego elementu:

 
$$ \basphi_3(x) = (x-x_2)/h,\quad \basphi_2(x) = 1- (x-x_2)/h$$

 
$$ A_{2,3} = \int_\Omega \basphi_2\basphi_{3}\dx = \int_{\xno{2}}^{\xno{3}} \left(1 - \frac{x - \xno{2}}{h}\right) \frac{x - x_{2}}{h} \dx = \frac{h}{6} $$

Przykład 2: Obliczenie wartości (diagonalnego) elementu macierzy

 
$$ A_{2,2} = \int_{\xno{1}}^{\xno{2}} \left(\frac{x - \xno{1}}{h}\right)^2\dx + \int_{\xno{2}}^{\xno{3}} \left(1 - \frac{x - \xno{2}}{h}\right)^2\dx = \frac{2h}{3} $$

Ogólna postać wzoru na wartość elementu \( A_{ij} \) - rysunek

 
$$ A_{i,i-1} = \int_\Omega \basphi_i\basphi_{i-1}\dx = \hbox{?}$$

Ogólna postać wzoru na wartość elementu \( A_{ij} \) - obliczenia

 
$$ \begin{align*} A_{i,i-1} &= \int_\Omega \basphi_i\basphi_{i-1}\dx\\ &= \underbrace{\int_{\xno{i-2}}^{\xno{i-1}} \basphi_i\basphi_{i-1}\dx}_{\basphi_i=0} + \int_{\xno{i-1}}^{\xno{i}} \basphi_i\basphi_{i-1}\dx + \underbrace{\int_{\xno{i}}^{\xno{i+1}} \basphi_i\basphi_{i-1}\dx}_{\basphi_{i-1}=0}\\ &= \int_{\xno{i-1}}^{\xno{i}} \underbrace{\left(\frac{x - x_{i}}{h}\right)}_{\basphi_i(x)} \underbrace{\left(1 - \frac{x - \xno{i-1}}{h}\right)}_{\basphi_{i-1}(x)} \dx = \frac{h}{6} \end{align*} $$

• \( A_{i,i+1}=A_{i,i-1} \) ze względu na symetrię
• \( A_{i,i}=2h/3 \) (obliczenia jak w przypadku \( A_{2,2} \)), z wyjątkiem:
• \( A_{0,0}=A_{N,N}=h/3 \) (całka tylko na jednym elemencie)

Obliczenia dla prawej strony równania

 
$$ b_i = \int_\Omega\basphi_i(x)f(x)\dx = \int_{\xno{i-1}}^{\xno{i}} \frac{x - \xno{i-1}}{h} f(x)\dx + \int_{x_{i}}^{\xno{i+1}} \left(1 - \frac{x - x_{i}}{h}\right) f(x) \dx $$

Do dalszych obliczeń potrzebna konkretna postać \( f(x) \) ...

Przykład: rozwiązanie dla obszaru dwu-elementowego – URL i rozwiązanie

• \( f(x)=x(1-x) \) na \( \Omega=[0,1] \)
• Dwa elementy o jednakowej długości: \( [0,0.5] \) oraz \( [0.5,1] \)

 
$$ \begin{equation*} A = \frac{h}{6}\left(\begin{array}{ccc} 2 & 1 & 0\\ 1 & 4 & 1\\ 0 & 1 & 2 \end{array}\right),\quad b = \frac{h^2}{12}\left(\begin{array}{c} 2 - 3h\\ 12 - 14h\\ 10 -17h \end{array}\right) \end{equation*} $$

 
$$ \begin{equation*} c_0 = \frac{h^2}{6},\quad c_1 = h - \frac{5}{6}h^2,\quad c_2 = 2h - \frac{23}{6}h^2 \end{equation*} $$

Przykład: rozwiązanie dla obszaru dwu-elementowego – rozwiązanie-rysunek

 
$$ \begin{equation*} u(x)=c_0\basphi_0(x) + c_1\basphi_1(x) + c_2\basphi_2(x)\end{equation*} $$

Przykład: Rozwiązanie dla obszaru 4-elementowego – rozwiązanie-rysunek

Przykład: elementy typu P2

Generowanie macierzy globalnej - logika obliczeń

(ang. assemble - gromadzić, składać, zbierać)

assembling, assemblacja?

Całkowanie z perspektywy elementu

 
$$ A_{i,j} = \int_\Omega\basphi_i\basphi_jdx = \sum_{e} \int_{\Omega^{(e)}} \basphi_i\basphi_jdx,\quad A^{(e)}_{i,j}=\int_{\Omega^{(e)}} \basphi_i\basphi_jdx $$

Ważne spostrzeżenia:

• \( A^{(e)}_{i,j}\neq 0 \) wtedy i tylko wtedy gdy węzeły \( i \) oraz \( j \) leżą na tym samym elemencie \( e \) (w przeciwnym wypadku nośniki funkcji to zbiory rozłączne)
• Wszystkie niezerowe współczynniki danego elementu \( A^{(e)}_{i,j} \) tworzą lokalną macierz dla danego elementu (element matrix)
• ''Wkład'' w macierz lokalną elementu mają wyłącznie funkcje bazowe związane z węzłami leżącymi na tym elemencie
• Wygodne rozwiązanie: wprowadzenie indeksacji lokalnej węzłów leżących na danym elemencie: \( 0,1,\ldots,d \)

Macierz elementów: indeksacja lokalna/indeksacja globalna

 
$$ \tilde A^{(e)} = \{ \tilde A^{(e)}_{r,s}\},\quad \tilde A^{(e)}_{r,s} = \int_{\Omega^{(e)}}\basphi_{q(e,r)}\basphi_{q(e,s)}dx, \quad r,s\in\Ifd=\{0,\ldots,d\} $$

Przykład: assembling macierzy dla kolejno ponumerowanych elementów P1

TODO

Przykład: assembling macierzy dla kolejno ponumerowanych elementów P3

TODO

Przykład: assembling macierzy dla nieregularnej siatki elementów P1

TODO

Assembling prawej strony układu

 
$$ b_i = \int_\Omega f(x)\basphi_i(x)dx = \sum_{e} \int_{\Omega^{(e)}} f(x)\basphi_i(x)dx,\quad b^{(e)}_{i}=\int_{\Omega^{(e)}} f(x)\basphi_i(x)dx $$

Ważne spostrzeżenia:

• \( b_i^{(e)}\neq 0 \) wtedy i tylko wtedy gdy węzeł globalny \( i \) leży na danym elemencie \( e \) (w przeciwnym przypadku \( \basphi_i=0 \))
• \( d+1 \) niezerowych wartości \( b_i^{(e)} \) może być zgromadzonych w lokalnym wektorze elementu e. \( \tilde b_r^{(e)}=\{ \tilde b_r^{(e)}\} \), \( r\in\Ifd \)

Assembling:

 
$$ b_{q(e,r)} := b_{q(e,r)} + \tilde b^{(e)}_{r},\quad r\in\Ifd $$

Transformacja współrzędnych globalnych do współrzędnych unormowanych

Normalizacja współrzędnych położenia:

Zamiast całkować w granicach \( [x_L, x_R] \)

 
$$ \begin{equation*} \tilde A^{(e)}_{r,s} = \int_{\Omega^{(e)}}\basphi_{q(e,r)}(x)\basphi_{q(e,s)}(x)dx = \int_{x_L}^{x_R}\basphi_{q(e,r)}(x)\basphi_{q(e,s)}(x)dx \end{equation*} $$

można transformować przedział \( [x_L, x_R] \) na przedział unormowany \( [-1,1] \) o współrzędnej lokalnej \( X \)

Transformacja liniowa \( X\in [-1,1] \) w \( x\in [x_L,x_R] \)

(Transformacja afiniczna)

 
$$ x = \half (x_L + x_R) + \half (x_R - x_L)X $$

inaczej

 
$$ x = x_m + {\half}hX, \qquad x_m=(x_L+x_R)/2,\quad h=x_R-x_L $$

Transformacja odwrotna:

 
$$ X = \frac{2 x + (x_L + x_R)}{(x_R - x_L)} $$

Transformacja całki

Zmiana granic całkowania -> całkowanie na przedziale unormowanym: podstawienie \( x(X) \) w miejsce \( x \).

Lokalne funkcje bazowe we współrzędnych unormowanych:

 
$$ \refphi_r(X) = \basphi_{q(e,r)}(x(X)) $$

 
$$ \begin{array}{c} x = \half (x_L + x_R) + \half (x_R - x_L)X \\ \downarrow \\ dx = \half (x_R - x_L) dX \end{array} $$

 
$$ \tilde A^{(e)}_{r,s} = \int_{\Omega^{(e)}}\basphi_{q(e,r)}(x)\basphi_{q(e,s)}(x)dx = \int\limits_{-1}^1 \refphi_r(X)\refphi_s(X)\underbrace{\frac{dx}{dX}}_{\det J= h/2}dX $$

 
$$ \tilde A^{(e)}_{r,s} = \int\limits_{-1}^1 \refphi_r(X)\refphi_s(X)\det J\,dX $$

 
$$ \tilde b^{(e)}_{r} = \int_{\Omega^{(e)}}f(x)\basphi_{q(e,r)}(x)dx = \int\limits_{-1}^1 f(x(X))\refphi_r(X)\det J\,dX $$

Zalety całkowania na przedziale unormowanym

• Całkowanie zawsze w tych samych granicach całkowania \( [-1,1] \)
• Potrzebne wzory tylko dla \( \refphi_r(X) \) na jednym elemencie (brak funkcji definiowanych na przedziałach (piecewise polynomial))
• Funkcja \( \refphi_r(X) \) jest taka sama dla wszystkich elementów niezależnie od ich położenia i rozmiarów (długości). Długość odcinka jest uwzględniona poprzez jakobian \( \det J \)

Funkcje bazowe P1 na elemencie unormowanym

 
$$ \begin{align} \refphi_0(X) &= \half (1 - X) \tag{10}\\ \refphi_1(X) &= \half (1 + X) \tag{11} \end{align} $$

(proste funkcje wielomianowe zamiast definicji funkcji na podprzedziałach

Funkcje bazowe P2 na elemencie unormowanym

 
$$ \begin{align} \refphi_0(X) &= \half (X-1)X \tag{12}\\ \refphi_1(X) &= 1 - X^2 \tag{13}\\ \refphi_2(X) &= \half (X+1)X \tag{14} \end{align} $$

Łatwość wygenerowania elementów dowolnego rzędu... Jak?

Sposoby znalezienia wzorów na funkcje bazowe

1. Transformacja globalnych funkcji bazowych \( \basphi_i(x) \) na element unormowany ze współrzędną \( X \)
2. Obliczenie \( \refphi_r(X) \)

• dla zadanego stopnia \( d \) szukamy wielomianów opartych o węzły wewnątrz przedziału \( [-1,1] \) o własności

• \( \refphi_r(X)=1 \) w węźle \( r \)
• \( \refphi_r(X)=0 \) we wszystkich pozostałych \( d \) węzłach

3. Wykorzystanie wzoru interpolacyjnego Lagrange'a

Całkowanie po elemencie unormowanym - lokalna macierz

Założenie: elementy typu P1, oraz funkcja \( f(x)=x(1-x) \).

 
$$ \begin{align} \tilde A^{(e)}_{0,0} &= \int_{-1}^1 \refphi_0(X)\refphi_0(X)\frac{h}{2} dX\nonumber\\ &=\int_{-1}^1 \half(1-X)\half(1-X) \frac{h}{2} dX = \frac{h}{8}\int_{-1}^1 (1-X)^2 dX = \frac{h}{3} \tag{15}\\ \tilde A^{(e)}_{1,0} &= \int_{-1}^1 \refphi_1(X)\refphi_0(X)\frac{h}{2} dX\nonumber\\ &=\int_{-1}^1 \half(1+X)\half(1-X) \frac{h}{2} dX = \frac{h}{8}\int_{-1}^1 (1-X^2) dX = \frac{h}{6} \tag{16}\\ \tilde A^{(e)}_{0,1} &= \tilde A^{(e)}_{1,0} \tag{17}\\ \tilde A^{(e)}_{1,1} &= \int_{-1}^1 \refphi_1(X)\refphi_1(X)\frac{h}{2} dX\nonumber\\ &=\int_{-1}^1 \half(1+X)\half(1+X) \frac{h}{2} dX = \frac{h}{8}\int_{-1}^1 (1+X)^2 dX = \frac{h}{3} \tag{18} \end{align} $$

Całkowanie po elemencie unormowanym - wektor prawej strony

 
$$ \begin{align} \tilde b^{(e)}_{0} &= \int_{-1}^1 f(x(X))\refphi_0(X)\frac{h}{2} dX\nonumber\\ &= \int_{-1}^1 (x_m + \half hX)(1-(x_m + \half hX)) \half(1-X)\frac{h}{2} dX \nonumber\\ &= - \frac{1}{24} h^{3} + \frac{1}{6} h^{2} x_{m} - \frac{1}{12} h^{2} - \half h x_{m}^{2} + \half h x_{m} \tag{19}\\ \tilde b^{(e)}_{1} &= \int_{-1}^1 f(x(X))\refphi_1(X)\frac{h}{2} dX\nonumber\\ &= \int_{-1}^1 (x_m + \half hX)(1-(x_m + \half hX)) \half(1+X)\frac{h}{2} dX \nonumber\\ &= - \frac{1}{24} h^{3} - \frac{1}{6} h^{2} x_{m} + \frac{1}{12} h^{2} - \half h x_{m}^{2} + \half h x_{m} \tag{20} \end{align} $$

\( x_m \): środek elementu

Obliczenia symboliczne zamiast żmudnego liczenia na kartce...

>>> import sympy as sym
>>> x, x_m, h, X = sym.symbols('x x_m h X')
>>> sym.integrate(h/8*(1-X)**2, (X, -1, 1))
h/3
>>> sym.integrate(h/8*(1+X)*(1-X), (X, -1, 1))
h/6
>>> x = x_m + h/2*X
>>> b_0 = sym.integrate(h/4*x*(1-x)*(1-X), (X, -1, 1))
>>> print b_0
-h**3/24 + h**2*x_m/6 - h**2/12 - h*x_m**2/2 + h*x_m/2

Implementacja

• Funkcje przedstawione na kolejnych slajdach znajdują sie w module fe_approx1D.py
• Przedstawione funkcje działają w trybie symbolicznym, jak i numerycznym
• Kod zawiera wszystkie kroku obliczeń elementami skończonymi.

Generowanie funkcji bazowych na przedziale unormowanymi

Niech \( \refphi_r(X) \) będzie wielomianem Lagrange'a stopnia d:

import sympy as sym
import numpy as np

def phi_r(r, X, d):
    if isinstance(X, sym.Symbol):
        h = sym.Rational(1, d)  # node spacing
        nodes = [2*i*h - 1 for i in range(d+1)]
    else:
        # assume X is numeric: use floats for nodes
        nodes = np.linspace(-1, 1, d+1)
    return Lagrange_polynomial(X, r, nodes)

def Lagrange_polynomial(x, i, points):
    p = 1
    for k in range(len(points)):
        if k != i:
            p *= (x - points[k])/(points[i] - points[k])
    return p

def basis(d=1):
    """Return the complete basis."""
    X = sym.Symbol('X')
    phi = [phi_r(r, X, d) for r in range(d+1)]
    return phi

Obliczanie współczynników macierzy

def element_matrix(phi, Omega_e, symbolic=True):
    n = len(phi)
    A_e = sym.zeros((n, n))
    X = sym.Symbol('X')
    if symbolic:
        h = sym.Symbol('h')
    else:
        h = Omega_e[1] - Omega_e[0]
    detJ = h/2  # dx/dX
    for r in range(n):
        for s in range(r, n):
            A_e[r,s] = sym.integrate(phi[r]*phi[s]*detJ, (X, -1, 1))
            A_e[s,r] = A_e[r,s]
    return A_e

Przykład: Obliczenia macierzy współczynników: symbolicznie vs numerycznie

>>> from fe_approx1D import *
>>> phi = basis(d=1)
>>> phi
[1/2 - X/2, 1/2 + X/2]
>>> element_matrix(phi, Omega_e=[0.1, 0.2], symbolic=True)
[h/3, h/6]
[h/6, h/3]
>>> element_matrix(phi, Omega_e=[0.1, 0.2], symbolic=False)
[0.0333333333333333, 0.0166666666666667]
[0.0166666666666667, 0.0333333333333333]

Obliczenia współczynników wektora prawej strony

def element_vector(f, phi, Omega_e, symbolic=True):
    n = len(phi)
    b_e = sym.zeros((n, 1))
    # Make f a function of X
    X = sym.Symbol('X')
    if symbolic:
        h = sym.Symbol('h')
    else:
        h = Omega_e[1] - Omega_e[0]
    x = (Omega_e[0] + Omega_e[1])/2 + h/2*X  # mapping
    f = f.subs('x', x)  # substitute mapping formula for x
    detJ = h/2  # dx/dX
    for r in range(n):
        b_e[r] = sym.integrate(f*phi[r]*detJ, (X, -1, 1))
    return b_e

Zwróć uwagę na f.subs('x', x) -> podstawienie \( x(X) \) za x w formule na f (od teraz f jest funkcją \( f(X) \))

Powrót do całkowania numerycznego w razie niepowodzenia całkowania symbolicznego \( \int f\refphi_r dx \)

• Macierz lewej strony: tylko wielomiany -> sympy zawsze da radę
• Wektor prawej strony: całkowanie \( \int f\refphi \dx \) może się nie powieść (sympy zwróci obiekt typu Integral zamiast liczby)

def element_vector(f, phi, Omega_e, symbolic=True):
        ...
        I = sym.integrate(f*phi[r]*detJ, (X, -1, 1))  # try...
        if isinstance(I, sym.Integral):
            h = Omega_e[1] - Omega_e[0]  # Ensure h is numerical
            detJ = h/2
            integrand = sym.lambdify([X], f*phi[r]*detJ)
            I = sym.mpmath.quad(integrand, [-1, 1])
        b_e[r] = I
        ...

Assembling URL i rozwiązanie

def assemble(nodes, elements, phi, f, symbolic=True):
    N_n, N_e = len(nodes), len(elements)
    zeros = sym.zeros if symbolic else np.zeros
    A = zeros((N_n, N_n))
    b = zeros((N_n, 1))
    for e in range(N_e):
        Omega_e = [nodes[elements[e][0]], nodes[elements[e][-1]]]

        A_e = element_matrix(phi, Omega_e, symbolic)
        b_e = element_vector(f, phi, Omega_e, symbolic)

        for r in range(len(elements[e])):
            for s in range(len(elements[e])):
                A[elements[e][r],elements[e][s]] += A_e[r,s]
            b[elements[e][r]] += b_e[r]
    return A, b

Rozwiązanie URL

if symbolic:
    c = A.LUsolve(b)           # sympy arrays, symbolic Gaussian elim.
else:
    c = np.linalg.solve(A, b)  # numpy arrays, numerical solve

Uwaga: obliczanie współczynników macierzy A, b oraz rozwiązanie URL A.LUsolve(b) może być baaardzo czasochłonne$\ldots$

Przykład: generowanie macierzy symbolicznie

>>> h, x = sym.symbols('h x')
>>> nodes = [0, h, 2*h]
>>> elements = [[0, 1], [1, 2]]
>>> phi = basis(d=1)
>>> f = x*(1-x)
>>> A, b = assemble(nodes, elements, phi, f, symbolic=True)
>>> A
[h/3,   h/6,   0]
[h/6, 2*h/3, h/6]
[  0,   h/6, h/3]
>>> b
[     h**2/6 - h**3/12]
[      h**2 - 7*h**3/6]
[5*h**2/6 - 17*h**3/12]
>>> c = A.LUsolve(b)
>>> c
[                           h**2/6]
[12*(7*h**2/12 - 35*h**3/72)/(7*h)]
[  7*(4*h**2/7 - 23*h**3/21)/(2*h)]

Przykład: generowanie macierzy numerycznie

>>> nodes = [0, 0.5, 1]
>>> elements = [[0, 1], [1, 2]]
>>> phi = basis(d=1)
>>> x = sym.Symbol('x')
>>> f = x*(1-x)
>>> A, b = assemble(nodes, elements, phi, f, symbolic=False)
>>> A
[ 0.166666666666667, 0.0833333333333333,                  0]
[0.0833333333333333,  0.333333333333333, 0.0833333333333333]
[                 0, 0.0833333333333333,  0.166666666666667]
>>> b
[          0.03125]
[0.104166666666667]
[          0.03125]
>>> c = A.LUsolve(b)
>>> c
[0.0416666666666666]
[ 0.291666666666667]
[0.0416666666666666]

Struktura macierzy współczynników

>>> d=1; N_e=8; Omega=[0,1]  # 8 linear elements on [0,1]
>>> phi = basis(d)
>>> f = x*(1-x)
>>> nodes, elements = mesh_symbolic(N_e, d, Omega)
>>> A, b = assemble(nodes, elements, phi, f, symbolic=True)
>>> A
[h/3,   h/6,     0,     0,     0,     0,     0,     0,   0]
[h/6, 2*h/3,   h/6,     0,     0,     0,     0,     0,   0]
[  0,   h/6, 2*h/3,   h/6,     0,     0,     0,     0,   0]
[  0,     0,   h/6, 2*h/3,   h/6,     0,     0,     0,   0]
[  0,     0,     0,   h/6, 2*h/3,   h/6,     0,     0,   0]
[  0,     0,     0,     0,   h/6, 2*h/3,   h/6,     0,   0]
[  0,     0,     0,     0,     0,   h/6, 2*h/3,   h/6,   0]
[  0,     0,     0,     0,     0,     0,   h/6, 2*h/3, h/6]
[  0,     0,     0,     0,     0,     0,     0,   h/6, h/3]

Uwaga (zadanie domowe): Wykonaj obliczenia na kartce papieru w celu potwierdzenia wartości poszczególnych elementów powyższej macierzy (pomocne w zrozumieniu materiału).

Wynik w przypadku ogólnym (\( N \) jednakowych elementów)

• Macierz rzadka -> większość współczynników to zera
• Przykład dla elementów typu P1, siatka regularna

 
$$ A = \frac{h}{6} \left( \begin{array}{cccccccccc} 2 & 1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ 1 & 4 & 1 & \ddots & & & & & \vdots \\ 0 & 1 & 4 & 1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & 1 & 4 & 1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & 1 & 4 & 1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 1 & 2 \end{array} \right) $$

Macierz rzadka dla elementów typu P2 (siatka regularna)

 
$$ A = \frac{h}{30} \left( \begin{array}{ccccccccc} 4 & 2 & - 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 2 & 16 & 2 & 0 & 0 & 0 & 0 & 0 & 0\\- 1 & 2 & 8 & 2 & - 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 16 & 2 & 0 & 0 & 0 & 0\\ 0 & 0 & - 1 & 2 & 8 & 2 & - 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 2 & 16 & 2 & 0 & 0\\ 0 & 0 & 0 & 0 & - 1 & 2 & 8 & 2 & - 1 \\0 & 0 & 0 & 0 & 0 & 0 & 2 & 16 & 2\\0 & 0 & 0 & 0 & 0 & 0 & - 1 & 2 & 4 \end{array} \right) $$

Macierz rzadka dla siatek regularnych/indeksowanych losowo dla elementów P1

• Po lewej: węzły i elementy ideksowane od lewej do prawej
• Po prawej: węzły i elementy indeksowane ''losowo''

Macierz rzadka dla siatek regularnych/indeksowanych losowo dla elementów P3

• Po lewej: węzły i elementy ideksowane od lewej do prawej
• Po prawej: węzły i elementy indeksowane ''losowo''

Macierze rzadkie – podsumowanie

Postać specyficznych macierzy \( A_{i,j} \):

• Elementy P1: 3 niezerowe elementy w wierszu
• Elementy P2: 5 niezerowe elementy w wierszu
• Elementy P3: 7 niezerowe elementy w wierszu

Wskazówki:

• Należy używać specjalne techniki przechowywania takich macierzy w pamięci i specjalnych solverów dla macierzy rzadkich
• W Pythonie: pakiet scipy.sparse

Przykład: przybliżenie funkcji \( f\sim x^9 \) elementami różnego typu; kod

Zadanie: Porównać rozwiązanie zadania przybliżenia funkcji \( f(x) \) przy pomocy siatki \( N_e \) elementów skończonych o funkcjach bazowych rzędu \( d \).

import sympy as sym
from fe_approx1D import approximate
x = sym.Symbol('x')

approximate(f=x*(1-x)**8, symbolic=False, d=1, N_e=4)
approximate(f=x*(1-x)**8, symbolic=False, d=2, N_e=2)
approximate(f=x*(1-x)**8, symbolic=False, d=1, N_e=8)
approximate(f=x*(1-x)**8, symbolic=False, d=2, N_e=4)

Przykład: przybliżenie funkcji \( f\sim x^9 \) elementami różnego typu; rysunki

Ograniczenia zaprezentowanego podejścia elementów skończonych

Najczęstsza interpretacja:

• Węzły: punkty potrzebne do zdefiniowania \( \basphi_i \) i obliczania

wartości \( u \) (wymagane do geometrii i aproksymacji funkcji)

• Elementy: podobszary (zawierające kilka węzłów)

Problem:

• węzły na brzegu potrzebne przy warunkach brzegowych, a nie zawsze tak musi być dla szczególnych rodzajów interpolacji (np. elementy stałe)
• Trzeba wymyśleć coś lepszego...

Uogólnienie koncepcji elementu skończonego (komórki, wierzchołki, węzły, stopnie swobody)

• Rozdzielenie aproksymacji geometrii (obszaru) od aproksymacji funkcji ''nad'' obszarem
• Nowe pojęcia: komórka (ang. cell) – podobszar, element, kawałek obszaru
• Komórka zbudowana jest z wierzchołków (ang. vertices <- vertex) – (krańców przedziału w 1D)
• Węzły (ang. nodes) p- punkty, w których należy wyznaczyć wartość poszukiwanej funkcji (nie muszą pokrywać się z wierzchołkami!, ale mogą\ldots)
• Stopnie swobody (ang. degrees of freedom)

  – wielkości reprezentowane przez \( c_j \) (niewiadome w URL) -> najczęściej: wartości funkcji w węźle \( \sum_{j\in\If} c_j \basphi_j(\xno{i}) = c_i \)

wierzchołki -> komórki -> interpolacja geometrii

węzły, stopnie swobody -> interpolacja funkcji

Pojęcie elementu skończonego

1. komórka odniesienia z unormowanym, lokalnym układem współrzędnych
2. zbiór funkcji bazowych \( \refphi_r \) dla komórki
3. zbiór stopni swobody (t.j. wartości funkcji), jednoznacznie wyznaczający funkcje bazowe, dobrane tak aby \( \refphi_r=1 \) dla \( r \)-tego stopnia swobody oraz \( \refphi_r=0 \) dla wszystkich pozostałych stopni swobody
4. odwzorowanie (ang. mapping) pomiędzy lokalną a globalną indeksacją (transformacja numeracji) stopni swobody (odwzorowanie dof – dof map)
5. odwzorowanie komórki unormowanej na komórkę rzeczywistego obszaru (w 1D: \( [-1,1]\ \Rightarrow\ [x_L,x_R] \))

Struktury danych: vertices, cells, dof_map

• Współrzędne wierzchołków komórek: vertices (równoważne strukturze nodes dla elementów P1)
• Wierzchołki dla elementów (komórek): cells[e][r] numer globalny dla wierchołka r elementu e (równoważne strukturze elements dla elementów typu P1)
• dof_map[e,r] odwzorowanie lokalnego indeksu stopnia swobody r elementu e na number globalny (równoważne strukturze elements dla elementów typu Pd)

W trakcie assemblingu należy skorzystać ze struktury dof_map:

A[dof_map[e][r], dof_map[e][s]] += A_e[r,s]
b[dof_map[e][r]] += b_e[r]

vertices = [0, 0.4, 1]
cells = [[0, 1], [1, 2]]
dof_map = [[0, 1, 2], [2, 3, 4]]

Przykład: elementy P0

Przykład: Ta sama siatka, ale \( u \) to funkcja stała na każdej komórce (przedziałami stała) -> elementy typu P0.

Te same struktury vertices i cells, ale dodatkowo

dof_map = [[0], [1]]

Można traktować te elementy jak elementy z interpolacją opartą na węźle znajdującym się pośrodku elementu.

Szkielet programu

# Use modified fe_approx1D module
from fe_approx1D_numint import *

x = sym.Symbol('x')
f = x*(1 - x)

N_e = 10
# Create mesh with P3 (cubic) elements
vertices, cells, dof_map = mesh_uniform(N_e, d=3, Omega=[0,1])

# Create basis functions on the mesh
phi = [basis(len(dof_map[e])-1) for e in range(N_e)]

# Create linear system and solve it
A, b = assemble(vertices, cells, dof_map, phi, f)
c = np.linalg.solve(A, b)

# Make very fine mesh and sample u(x) on this mesh for plotting
x_u, u = u_glob(c, vertices, cells, dof_map,
                resolution_per_element=51)
plot(x_u, u)

Przybliżenie paraboli elementami P0

Funkcja approximate ''opakowuje'' polecenia z poprzedniego slajdu:

from fe_approx1D_numint import *
x=sym.Symbol("x")
for N_e in 4, 8:
    approximate(x*(1-x), d=0, N_e=N_e, Omega=[0,1])

Obliczanie błędów aproksymacji; uwagi ogólne

Błąd jako funkcja:

 
$$ e(x) = f(x) - u(x) $$

Błąd – dyskretna wartość -> normy:

 
$$ L^2 \hbox{ error: }\quad ||e||_{L^2} = \left(\int_{\Omega} e^2 dx\right)^{1/2}$$

Szacowanie całki:

• dokładne, analityczne (symboliczne) - nieuniwersalne -> kwadratury
• odpowiednio dokładne spróbkowanie \( u(x) \) w wielu punktach każdego elementu (np. poprzez wywołanie u_glob, które zwróci x i u), a następnie
• scałkowanie metodą trapezów
• Uwaga! Ważne! Całka powinna być policzona dokładnie ''po elementach'' (zmienność funkcji \( f \))

Obliczanie błędów aproksymacji; szczegóły

# Given c, compute x and u values on a very fine mesh
x, u = u_glob(c, vertices, cells, dof_map,
              resolution_per_element=101)
# Compute the error on the very fine mesh
e = f(x) - u
e2 = e**2
# Vectorized Trapezoidal rule
E = np.sqrt(0.5*np.sum((e2[:-1] + e2[1:])*(x[1:] - x[:-1]))

Zależność błędu od \( h \) i \( d \)

Teoria i eksperymenty pokazują, że aplikacja \( \LSM \) czy metody Galerkina dla elementów skończonych typu Pd o tej samej długości \( h \) daje błąd:

 
$$ ||e||_{L^2} = C | f^{d+1} | h^{d+1} $$

gdzie \( C \) zależy od \( d \) i \( \Omega = [0, L] \) ale nie zależy od \( h \), oraz

 
$$ |f^{d+1}|^{2} = \int_0^L \left( \frac{d^{d+1}f}{d x^{d+1}} \right)^2 dx $$

Kubiczne wielomiany Hermite'a - definicja

• Czy da się skonstruować \( \basphi_i(x) \) z ciągłą pochodną? Tak!

Niech dana będzie unormowana komórka \( [-1,1] \) z dwoma węzłami \( X=-1 \) i \( X=1 \). Stopnie swobody:

• 0: wartość funkcji w \( X=-1 \)
• 1: wartość pierwszej pochodnej w \( X=-1 \)
• 2: wartość funkcji w \( X=1 \)
• 3: wartość pierwszej pochodnej w \( X=1 \)

Kubiczne wielomiany Hermite'a - wyprowadzenie

4 ogranicznia na \( \refphi_r \) (1 dla stopnia swobody \( r \), 0 dla pozostałych):

• \( \refphi_0(\Xno{0}) = 1 \), \( \refphi_0(\Xno{1}) = 0 \), \( \refphi_0'(\Xno{0}) = 0 \), \( \refphi_0'(\Xno{1}) = 0 \)
• \( \refphi_1'(\Xno{0}) = 1 \), \( \refphi_1'(\Xno{1}) = 0 \), \( \refphi_1(\Xno{0}) = 0 \), \( \refphi_1(\Xno{1}) = 0 \)
• \( \refphi_2(\Xno{1}) = 1 \), \( \refphi_2(\Xno{0}) = 0 \), \( \refphi_2'(\Xno{0}) = 0 \), \( \refphi_2'(\Xno{1}) = 0 \)
• \( \refphi_3'(\Xno{1}) = 1 \), \( \refphi_3'(\Xno{0}) = 0 \), \( \refphi_3(\Xno{0}) = 0 \), \( \refphi_3(\Xno{1}) = 0 \)

Cztery układy równań liniowych z 4 niewiadomymi - współczynnikami wielomianów 3 stopnia.

Kubiczne wielomiany Hermite'a - wynik

 
$$ \begin{align} \refphi_0(X) &= 1 - \frac{3}{4}(X+1)^2 + \frac{1}{4}(X+1)^3 \tag{21}\\ \refphi_1(X) &= -(X+1)(1 - \half(X+1))^2 \tag{22}\\ \refphi_2(X) &= \frac{3}{4}(X+1)^2 - \half(X+1)^3 \tag{23}\\ \refphi_3(X) &= -\half(X+1)(\half(X+1)^2 - (X+1)) \tag{24}\\ \tag{25} \end{align} $$

Kubiczne wielomiany Hermite'a - sprawdzenie

# definition of the interval ends
x = np.array([-1, 1]) 

C = []  # list of polynomials stored as coefficients 
B = []  # list of basis functions
dB = [] # list of the derivatives of basis functions

for k in np.arange(0,4):
    A = np.array(  [[   x[0]**3,   x[0]**2, x[0], 1 ],
                    [ 3*x[0]**2, 2*x[0],       1, 0 ],
                    [   x[1]**3,   x[1]**2, x[1], 1 ],
                    [ 3*x[1]**2, 2*x[1],       1, 0 ]])

    b = np.zeros( (4,1) ); b[k] = 1 
    
    c = np.linalg.solve(A, b); C.append( c )

    B.append( lambda x:     C[k][0,0] * x**3 +   C[k][1,0] * x**2 + C[k][2,0] * x + C[k][3,0] )
    dB.append( lambda x: 3* C[k][0,0] * x**2 + 2*C[k][1,0] * x + C[k][2,0] )

Kubiczne wielomiany Hermite'a - sprawdzenie

# Check numerically that resulting cubic polynomial
# fulfills imposed requirements
A = [1,       1,       2,      1]      # basis function coefficients 
#    U(x[0])  dU(x[0]) U(x[1]) dU(x[1])

xx = np.arange(-1,1, 0.001)  
U = np.zeros(xx.shape)

for k in np.arange(0,4):
    U = U + A[k] * B[k](xx)

# numerical approximation of the derivatives at the ends of the interval
dl = (U[1]-U[0])/(xx[1]-xx[0])
dr = (U[-1]-U[-2])/(xx[-1]-xx[-2])

numericalApproximationOfA = [ U[0], dl, U[-1], dr]

print(numericalApproximationOfA)

Wynik działania skryptu:

[1.0, 0.9992502500000269, 1.999000749750002, 0.9977517500000526]

Kubiczne wielomiany Hermite'a - wyniki

Całkowanie numeryczne

• \( \int_\Omega f\basphi_idx \) - konieczność całkowania numerycznego
• Współczynniki macierzy lewej strony - zwykle również numerycznie (bo wygodnie)

Ogólna postać kwadratury

 
$$ \int_{-1}^{1} g(X)dX \approx \sum_{j=0}^M w_jg(\bar X_j), $$

gdzie

• \( \bar X_j \) to węzły kwadratury
• \( w_j \) – wagi kwadratury

Różne metody -> różny wybór węzłów i wag

Wzór prostokątów

(ang. midpoint rule) – metoda punktu środkowego

Najprostsza metoda

 
$$ \int_{-1}^{1} g(X)dX \approx 2g(0),\quad \bar X_0=0,\ w_0=2, $$

Dokładna dla funkcji podcałkowych będących wielomianami 1. stopnia

Metody Newtona-Cotesa

• Idea: węzły kwadratury równomiernie rozmieszczone na \( [-1,1] \)
• Węzły kwadratury często pokrywają się węzłami siatki

Wzór trapezów:

 
$$ \int_{-1}^{1} g(X)dX \approx g(-1) + g(1),\quad \bar X_0=-1,\ \bar X_1=1,\ w_0=w_1=1, $$

Wzór Simpsona (parabol):

 
$$ \int_{-1}^{1} g(X)dX \approx \frac{1}{3}\left(g(-1) + 4g(0) + g(1)\right), $$

gdzie

 
$$ \bar X_0=-1,\ \bar X_1=0,\ \bar X_2=1,\ w_0=w_2=\frac{1}{3},\ w_1=\frac{4}{3} $$

Metoda Gaussa-Legendre'a

• optymalne położenie węzłów kwadratury -> wyższa dokładność
• Kwadratury Gaussa-Legendre'a -> dobranie położenia węzłów oraz wag tak, aby całkować z jak najlepszą dokładnością

 
$$ \begin{align} M=1&:\quad \bar X_0=-\frac{1}{\sqrt{3}},\ \bar X_1=\frac{1}{\sqrt{3}},\ w_0=w_1=1 \tag{26}\\ M=2&:\quad \bar X_0=-\sqrt{\frac{3}{{5}}},\ \bar X_0=0,\ \bar X_2= \sqrt{\frac{3}{{5}}},\ w_0=w_2=\frac{5}{9},\ w_1=\frac{8}{9} \tag{27} \end{align} $$

• \( M=1 \): dokładna dla wielomianów 3. stopnia
• \( M=2 \): dokładna dla wielomianów 5. stopnia
• W ogólności, \( M \)-punktowy wzór Gaussa-Legendre'a jest dokładny dla wielomianów stopnia \( 2M+1 \).

Plik `numint.py zawiera zbiór węzłów i wag dla metody Gaussa-Legendre'a.

Aproksymacja funkcji w 2D

Krótkie omówienie zagadnienia w 2D

Iloczyn skalarny w 2D:

 
$$ (f,g) = \int_\Omega f(x,y)g(x,y) dx dy $$

Zastosowanie \( \LSM \) lub metody Galerkina da URL:

 
$$ \begin{align*} \sum_{j\in\If} A_{i,j}c_j &= b_i,\quad i\in\If\\ A_{i,j} &= (\baspsi_i,\baspsi_j)\\ b_i &= (f,\baspsi_i) \end{align*} $$

Problem: Jak skonstruować dwuwymiarowe funkcje bazowe \( \baspsi_i(x,y) \)?

Funkcje bazowe 2D jako iloczyn tensorowy funkcji 1D

Korzystając z funkcji bazowych 1D zmiennej \( x \) oraz funkcji bazowych 1D zmiennej \( y \):

 
$$ \begin{align} V_x &= \mbox{span}\{ \hat\baspsi_0(x),\ldots,\hat\baspsi_{N_x}(x)\} \tag{28}\\ V_y &= \mbox{span}\{ \hat\baspsi_0(y),\ldots,\hat\baspsi_{N_y}(y)\} \tag{29} \end{align} $$

Przestrzeń wektorowa 2D może być zdefinowana jako iloczyn tensorowy \( V = V_x\otimes V_y \) z funkcjami bazowymi:

 
$$ \baspsi_{p,q}(x,y) = \hat\baspsi_p(x)\hat\baspsi_q(y) \quad p\in\Ix,q\in\Iy\tp $$

Iloczyn tensorowy

Niech dane będą dwa wektory \( a=(a_0,\ldots,a_M) \) i \( b=(b_0,\ldots,b_N) \). Ich zewnętrznym iloczynem tensorowym ( iloczynem diadycznym jeśli \( N=M \)) jest \( p=a\otimes b \) zdefiniowane jako:

 
$$ p_{i,j}=a_ib_j,\quad i=0,\ldots,M,\ j=0,\ldots,N\tp$$

Uwaga: \( p \) to macierz/tablica dwuwymiarowa

Przykład: baza 2D jako iloczyn tensorowy przestrzeni 1D:

 
$$ \baspsi_{p,q}(x,y) = \hat\baspsi_p(x)\hat\baspsi_q(y), \quad p\in\Ix,q\in\Iy$$

 
iloczyn tensorowy macierzy

(dowolnych wymiarów) ->

(dowolnych wymiarów) ->

(tego samego wymiaru)

Równoważność notacji z dwoma lub jednym indeksem

Baza przestrzeni 2D wymaga dwóch indeksów (i podwójnego sumowania) :

 
$$ u = \sum_{p\in\Ix}\sum_{q\in\Iy} c_{p,q}\baspsi_{p,q}(x,y) $$

Lub tylko jednego indeksu

 
$$ u = \sum_{j\in\If} c_j\baspsi_j(x,y)$$

jeśli posiadamy odwzorowanie \( (p,q)\rightarrow i \):

 
$$ \baspsi_i(x,y) = \hat\baspsi_p(x)\hat\baspsi_q(y), \quad i=p (N_y+1) + q\hbox{ or } i=q (N_x+1) + p $$

Przykładowa baza przestrzeni 2D; wzory

Dla dwuelementowej bazy 1D

 
$$ \{ 1, x \} $$

iloczyn tensorowy (wszystkie kombinacje) generuje bazę przestrzeni 2D:

 
$$ \baspsi_{0,0}=1,\quad \baspsi_{1,0}=x, \quad \baspsi_{0,1}=y, \quad \baspsi_{1,1}=xy $$

W notacji jednoindeksowej:

 
$$ \baspsi_0=1,\quad \baspsi_1=x, \quad \baspsi_2=y,\quad\baspsi_3 =xy $$

Przykładowa baza przestrzeni 2D; zastosowanie

 
$$ \baspsi_0=1,\quad \baspsi_1=x, \quad \baspsi_2=y,\quad\baspsi_3 =xy $$

 
$$ \int_0^{L_y} \int_0^{L_x} \left\{ \left[ \begin{array}{cccc} 1\cdot 1 & 1 \cdot x & 1 \cdot y & 1 \cdot xy \\ x\cdot 1 & x \cdot x & x \cdot y & x \cdot xy \\ y\cdot 1 & y \cdot x & y \cdot y & y \cdot xy \\ xy\cdot 1 & xy \cdot x & xy \cdot y & xy \cdot xy \end{array} \right] \left[ \begin{array}{c} c_0 \\ c_1 \\ c_2 \\ c_3 \end{array} \right] = \left[ \begin{array}{c} 1\cdot f(x,y) \\ x\cdot f(x,y) \\ y\cdot f(x,y) \\ xy\cdot f(x,y) \end{array} \right] \right\} dxdy $$

Funkcja aproksymowana (kwadratowa) \( f(x,y) = (1+x^2)(1+2y^2) \) (po lewej), funkcja aproksymująca (biliniowa) \( u \) (po prawej) (\( x^2y^2 \) vs \( xy \)):

Implementacja; principal changes to the 1D code

Zmiany w kodzie w stosunku do wersji 1D (approx1D.py):

• Omega = [[0, L_x], [0, L_y]]
• Całkowanie symboliczne w 2D
• Generowanie funkcji bazowych 2D (jako iloczynów tensorowych)

Implementacja 2D: całkowanie

import sympy as sym

integrand = psi[i]*psi[j]
I = sym.integrate(integrand,
                 (x, Omega[0][0], Omega[0][1]),
                 (y, Omega[1][0], Omega[1][1]))

# Fall back on numerical integration if symbolic integration
# was unsuccessful
if isinstance(I, sym.Integral):
    integrand = sym.lambdify([x,y], integrand)
    I = sym.mpmath.quad(integrand,
                       [Omega[0][0], Omega[0][1]],
                       [Omega[1][0], Omega[1][1]])

Implementacja 2D: funkcje bazowe

Iloczyn tensorowy bazy potęgowej \( x^i \) (bazy Taylora):

def taylor(x, y, Nx, Ny):
    return [x**i*y**j for i in range(Nx+1) for j in range(Ny+1)]

Iloczyn tensorowy bazy sinusoidalnej \( \sin((i+1)\pi x) \):

def sines(x, y, Nx, Ny):
    return [sym.sin(sym.pi*(i+1)*x)*sym.sin(sym.pi*(j+1)*y)
            for i in range(Nx+1) for j in range(Ny+1)]

Cały kod w approx2D.py.

Implementacja 2D: zastosowanie

>>> from approx2D import *
>>> f = (1+x**2)*(1+2*y**2)
>>> psi = taylor(x, y, 1, 1)
>>> Omega = [[0, 2], [0, 2]]
>>> u, c = least_squares(f, psi, Omega)
>>> print u
8*x*y - 2*x/3 + 4*y/3 - 1/9
>>> print sym.expand(f)
2*x**2*y**2 + x**2 + 2*y**2 + 1

Implementacja 2D: przykład zastosowanie bazy umożliwiającej konstrukcję rozwiązania dokładnego

Dodajemy funkcje bazowe o wyższych potęgach tak, aby \( f\in V \). Spodziewany wynik: \( u=f \)

>>> psi = taylor(x, y, 2, 2)
>>> u, c = least_squares(f, psi, Omega)
>>> print u
2*x**2*y**2 + x**2 + 2*y**2 + 1
>>> print u-f
0

Uogólnienie do zagadnień 3D

Kluczowa idea:

 
$$ V = V_x\otimes V_y\otimes V_z$$

W szczególności dla 3D:

 
$$ \begin{align*} \baspsi_{p,q,r}(x,y,z) &= \hat\baspsi_p(x)\hat\baspsi_q(y)\hat\baspsi_r(z)\\ u(x,y,z) &= \sum_{p\in\Ix}\sum_{q\in\Iy}\sum_{r\in\Iz} c_{p,q,r} \baspsi_{p,q,r}(x,y,z) \end{align*} $$

Elementy skończone w 2D i 3D

Zalety \( \FEM \) w zastosowaniach 2D i 3D:

• łatwość aproksymowania skomplikowanych geometrii
• łatwość generowania wielomianów (funkcji bazowych) wyższych rzędów w celu zwiększenia dokładności aproksymacji funkcji

\( \FEM \) w 1D: głównie dla celów dydaktycznych, debugowania

Przykłady komórek 2D i 3D

2D:

• trójkąty (triangles)
• czworokąty (quadrilaterals)

3D:

• czworościany (tetrahedra)
• sześciościany (hexahedra)

Obszar prostokątny (2D) zbudowany z elementów typu P1

Nieregularny obszar 2D zbudowany z elementów typu P1

Obszar prostokątny (2D) zbudowany z elementów typu Q1

Aproksymacja funkcji 2D na siatce elementów trójkątnych

Element trójkątny typu P1: aproksymacja \( u \) na każdym elemencie (komórce) funkcją liniową \( ax + by + c \)

Własności elementów 2D typu P1

• Komórki = trójkąty
• Wierzchołki = wierzchołki komórek
• węzły = wierzchołki trójkąta
• Stopnie swobody = wartości funkcji w węzłach
• \( \refphi_r(X,Y) \) jest funkcją liniową na komórce unormowanej
• \( \basphi_i(x,y) \) jest odzorowaniem \( \refphi_r(X,Y) \) na komórce rzeczywistej

Odwzorowanie liniowe elementu unormowanego na komórkę trójkątną

\( \basphi_i \): funkcja-piramida

• \( \basphi_i(x,y) \) – zmienność liniowa na poszczególnych komórkach
• \( \basphi_i=1 \) w wierzchołku (węźle) \( i \), 0 w pozostałych wierzchołkach (węzłach)

Elementy macierzy i wektora prawej strony

• Jak w 1D, wkład pojedynczej komórki do macierzy globalnego URL ogranicza się do kilku wartości w macierzy i wektorze wyrazów wolnych
• \( \basphi_i\basphi_j\neq 0 \) wtedy i tylko wtedy gdy \( i \) oraz \( j \) są stopniami swobody (wierzchołkami/węzłami) na tym samym elemencie
• Lokalna macierz trójkątengo elementu P1 to macierz o rozmiarach \( 3\times 3 \)

Funkcje bazowe na unormowanym elemencie trójkątnym

 
$$ \begin{align} \refphi_0(X,Y) &= 1 - X - Y \tag{30}\\ \refphi_1(X,Y) &= X \tag{31}\\ \refphi_2(X,Y) &= Y \tag{32} \end{align} $$

Funkcje bazowe \( \refphi_r \) wyższych stopni opierają się na większej liczbie węzłów (stopni swobody)

Elementy trójkątne typu P1, P2, P3, P4, P5, P6 przestrzeni 2D

Elementy P1 przestrzeni 1D, 2D i 3D

Elementy P2 przestrzeni 1D, 2D i 3D

• Interval, triangle, tetrahedron: sympleksy (ang. simplex -> *simplices*/*simplexes*)
• ściana (ang. face) – bok komórki (ścianka/krawedź/punkt)
• W czworościanie również krawędzie (edges)

Odwzorowanie afiniczne komórki unormowanej – wzór

Transformacja (Odwzorowanie) komórki we współrzędnych unormowanych

do komórki we współrzędnych globalnych:

 
$$ \begin{equation} \x = \sum_{r} \refphi_r^{(1)}(\X)\xdno{q(e,r)} \tag{33} \end{equation} $$

gdzie

• \( r \) przebiega przez wszystkie wierzchołki komórki
• \( \xdno{i} \) to globalne współrzędne \( (x,y) \) wierzchołka \( i \)
• \( \refphi_r^{(1)} \) to funkcja bazowa typu P1

Odwzorowanie zachowuje liniowość ścian i krawędzi.

• TODO (Przykład rachunkowy)

Odwzorowanie afiniczne komórki unormowanej

Komórki izoparametryczne

Idea: Wykorzystanie funkcji bazowych elementu (nie tylko funkcji typu P1 ale i wyższych rzędów) do odwzorowania geometrii:

 
$$ \x = \sum_{r} \refphi_r(\X)\xdno{q(e,r)} $$

Zaleta: pozwala generować elementy o geomtrii nielinowej

Obliczanie całek

Wymagana transformacja całek z \( \Omega^{(e)} \) (komórka we współrządnych globalnych) w \( \tilde\Omega^r \) (komórka unormowana/odniesienia):

 
$$ \begin{align} \int_{\Omega^{(e)}}\basphi_i (\x) \basphi_j (\x) \dx &= \int_{\tilde\Omega^r} \refphi_i (\X) \refphi_j (\X) \det J\, \dX \tag{34}\\ \int_{\Omega^{(e)}}\basphi_i (\x) f(\x) \dx &= \int_{\tilde\Omega^r} \refphi_i (\X) f(\x(\X)) \det J\, \dX \tag{35} \end{align} $$

gdzie \( \dx = dx dy \) lub \( \dx = dxdydz \) oraz \( \det J \) to wyznacznik jakobianu odwzorowania \( \x(\X) \).

 
$$ J = \left[\begin{array}{cc} \frac{\partial x}{\partial X} & \frac{\partial x}{\partial Y}\\ \frac{\partial y}{\partial X} & \frac{\partial y}{\partial Y} \end{array}\right], \quad \det J = \frac{\partial x}{\partial X}\frac{\partial y}{\partial Y} - \frac{\partial x}{\partial Y}\frac{\partial y}{\partial X} $$

Odwzorowanie afiniczne (33): \( \det J=2\Delta \), \( \Delta = \hbox{powierzchnia komórki/elementu} \)

Uwaga dot. uogólnienia FEM z 1D do 2D/3D