Elementos Finitos en \(d=1\)

Métodos Numéricos Avanzados. 2022-2023

Author

David Gómez-Castro

Formulación débil

Integración por partes

Recordamos que \[ \frac{\partial}{\partial x} \left( v \frac{\partial u}{\partial x} \right) = \frac{\partial v}{\partial x} \frac{\partial u}{\partial x} + v \frac{\partial^2 u}{\partial x^2} \]

Integrando \[ v(1) \frac{\partial u}{\partial x}(t,1) - v(-1) \frac{\partial u}{\partial x}(t,-1) = \int_{-1}^1 \frac{\partial v}{\partial x} (x) \frac{\partial u}{\partial x}(t,x) \diff x + \int_{-1}^1 v(x) \frac{\partial^2 u}{\partial x^2} (t,x) \diff x \]

Un momento de teoría

Si \(v\) es una función Lipschitz, entonces es derivable en casi todo punto. Y denotamos \(\frac{\diff v}{\diff x}\).

Esta derivada (a veces llamada “débil”) sigue teniendo la propiedad de integración por partes: si \(u\) es de clase \(C^2\) entonces \[ v(1) \frac{\partial u}{\partial x}(t,1) - v(-1) \frac{\partial u}{\partial x}(t,-1) = \int_{-1}^1 \frac{\partial v}{\partial x} (x) \frac{\partial u}{\partial x}(t,x) \diff x + \int_{-1}^1 v(x) \frac{\partial^2 u}{\partial x^2} (t,x) \diff x \]

El espacio de funciones Lipschitz es un ejemplo de los llamados espacios de Sobolev. Denotaremos \[ W^{1,\infty}_0 (-1,1) = \{ v \in C([-1,1]) \text{ Lipschitz tal que } v(-1) = v(1) = 0 \}. \]

Ejemplo. Ecuación de calor con condiciones Dirichlet

Volvemos sobre la ecuación de calor \[ \tag{C} \label{eq:calor} \frac{\partial u}{\partial t} = \Delta u + f(t,x) \]

Tomemos \(\eqref{eq:calor}\) en \((-1,1)\) con condiciones de Dirichlet homogéneas \(u(t,-1) = u(t,1) = 0\). La solución exacta de la ecuación satisface que para toda \(v\) Lipscthitz se tiene

\[ \int_{-1}^1 \frac{\partial u}{\partial t} v = \int_{-1}^1 (\Delta u) v + \int_{-1}^1 f v \]

Si \(v(\pm 1) = 0\) podemos integrar por partes llegamos a la llamada llamada formulación débil \[ \int_{-1}^1 \frac{\partial u}{\partial t} v = - \int_{-1}^1 \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v , \qquad \forall v \in W_0^{1,\infty} (\Omega). \]

Ejemplo. Ecuación de calor con condiciones Neumann

Tomemos \(\eqref{eq:calor}\) en \((-1,1)\) con condiciones de Neumann homogéneas \(\frac{\partial u}{\partial x}(t,-1) = \frac{\partial u}{\partial x}(t,1) = 0\). La solución exacta de la ecuación satisface que para toda \(v\) Lipscthitz se tiene

\[ \int_{-1}^1 \frac{\partial u}{\partial t} v = \int_{-1}^1 (\Delta u) v + \int_{-1}^1 f v \]

Podemos integrar por partes incluso si \(v(\pm 1) \ne 0\) a la llamada llamada formulación débil \[ \int_{-1}^1 \frac{\partial u}{\partial t} v = - \int_{-1}^1 \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v , \qquad \forall v \in W^{1,\infty} (\Omega). \]

Método de Galerkin

Un problema en dimensión finita

Elegimos funciones de ``base’’ \[ \varphi^h_1, \cdots, \varphi^h_N \]

y queremos buscar una aproximación \(u(t,x) \sim u^{h} (t,x)\) donde \[ \tag{U} \label{eq:uh} u^{h} (t,x)= \sum_{j=1}^N u_j^h (t) \varphi^h_j (x). \]

Estas funciones \(\varphi^h_i\) generan un sub-espacio vectorial en el espacio de funciones, que llamamos \(V^h\). Se busca \[ u^h (t) \in V^h = \mathrm{span} \{ \varphi^h_1, \cdots, \varphi^h_N \}. \]

Sin pérdida de generalidad consideramos que son linealmente independiente, eliminando en las que sean necesario para construir un base, y cambiándolo los índices.

Ejemplo. Ecuación de calor con condiciones Dirichlet

Volvemos sobre \(\eqref{eq:calor}\) en \(t > 0\) y \(x \in (-1,1)\) con condiciones de Dirichlet homogéneas \(u(t,-1) = u(t,1) = 0\).

La solución exacta es la única curva \(t \mapsto u(t) \in C^2([-1,1] ) \cap C_0([-1,1])\) tal que \[ \int_{-1}^1 \frac{\partial u}{\partial t} v = - \int_{-1}^1 \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v , \qquad \forall v \in W_0^{1,\infty}([-1,1]) \]

Como aproximación numérica buscamos un espacio \(V^h \subset W_0^{1,\infty} (-1,1)\) y la única curva \(t \mapsto u^h(t) \in V^h\) tal que \[ \int_{-1}^1 \frac{\partial u^h}{\partial t} v = - \int_{-1}^1 \frac{\partial u^h}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v \diff x , \qquad \forall v \in V^h. \]

Como \(V^h\) es un espacio de dimensión finita, esto es equivalente a que

\[ \int_{-1}^1 \frac{\partial u^h}{\partial t} \varphi_i^h = - \int_{-1}^1 \frac{\partial u^h}{\partial x} \frac{\partial \varphi_i^h}{\partial x} + \int_{-1}^1 f \varphi_i^h \diff x , \qquad \forall i = 1, \cdots, N. \]

Volviendo sobre \(\eqref{eq:uh}\) tenemos \[ \sum_{j=1}^N \frac{\diff u_j^h}{\diff t} (t) \int_{-1}^1 \varphi_j^h \varphi_i^h = - \sum_{j=1}^N u_j^h(t) \int_{-1}^1 \frac{\partial \varphi_j^h}{\partial x} \frac{\partial \varphi_i^h}{\partial x} + \int_{-1}^1 f(t,x) \varphi_i^h (x) \diff x. \]

Llegamos a un sistema de la forma \[ M^h \frac{\diff }{\diff t} \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = - A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} + b^h (t) \]

donde \[ M^h_{ij} = \int_{-1}^1 \varphi_j^h \varphi_i^h , \qquad A_{ij}^h = \int_{-1}^1 \frac{\diff \varphi_j^h}{\diff x} \frac{\diff \varphi_i^h}{\diff x}, \qquad b_i^h (t) = \int_{-1}^1 f(t,x) \varphi_i^h (x) \diff x . \]

Dado que los elementos de la base son linealmente independiente, la matriz \(M^h\) es no singular. Por tanto, podemos escribir \[ \frac{\diff }{\diff t} \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = - (M^h)^{-1} A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} + (M^h)^{-1} b^h(t) \] Este problema evolutivo tiene siempre solución.

Concluímos así la semi-discretización espacial. Vamos a estudiar las buenas propiedades que tiene este método, empezando por el problema estacionario asociado.

Ejemplo: ecuación de Laplace con condiciones Dirichlet

Consiste en encontrar soluciones estacionarias (e.d. independientes del tiempo) de la ecuación de calor, es decir buscamos un único vector \(u^h \in \mathbb R^N\) tal que \[ A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = b^h \]

Marco abstracto

Sea \(A: H \to H\) un operador lineal, donde \(H\) tiene un producto escalar \((\cdot, \cdot)_H\).

Para estudiar el problema evolutivo

\[ \frac{\diff}{\diff t} u = - A u \]

Siempre podemos escoger una base y escribir \(\eqref{eq:uh}\).

Repitiendo el razonamiento anterior, donde \(Au = - \Delta u - f\), tenemos \[ M^h \frac{\diff }{\diff t} \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = - A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} \]

La matriz \(M^h\) se suele llamar matriz de masa, y \(A^h\) se llama matriz de rigidez. Esto tiene que ver con los problemas de cuarto orden que aparecen en elasticidad.

donde \[ M^h_{ij} = (\varphi_j^h, \varphi_i^h)_H , \qquad A_{ij}^h = - (A \varphi_j, \varphi_i)_H. \]

Observación Se puede extender esta formulación a un problema no-lineal.

Método de elementos finitos

El método de elementos finitos consiste en:

  • Dividir el dominio en poliedros (es decir, crear una malla)

  • Tomar como funciones de base funciones que son polinomios de un grado \(g\) dado en cada elemento

Para esta sección recomendamos la bibliografía: Ramos (2012), Süli (2020)

En dimensión 1

En dimensión 1, esto significa tomar puntos \[ a = x_0^h < x_1^h < \dots < x_{M+1}^h = b. \]

Sobre esta malla se toma el espacio de las funciones de polinomios de grado \(g\) a trozos.

Si introducimos condiciones de tipo Dirichlet, entonces añadimos esta condiciones al espacio: \[ V^h = \{ v^h \in C_0 ([-1,1]) : v^h|_{[x_i, x_{i+1}]} \text{ es polinomio de grado } g \} \] Este espacio se denota, a veces, como \(V_g^h\).

La base de funciones de funciones lineales a trozos

Es muy habitual tomar \(g = 1\), es decir funciones continuas y lineales a trozos.

Como las funciones son lineales, para determinar una de ellas nos basta con determinar los valores \(v^h(x_i)\).

Elegimos la base de funciones tales que \(\varphi_i(x_j) = \delta_{ij}\).

En tal caso, se cogen las funciones de base lineales a trozos \[ \varphi_i (x) = \begin{dcases} \left( \frac{x-x_{i-1}}{x_i - x_{i-1}} \right)_+ & x \le x_i , \\ \left( \frac{x_{i+1}-x}{x_{i+1} - x_i} \right)_+ & x > x_i , \\ \end{dcases}, \]

Matriz de rigidez para el Laplaciano con condiciones Dirichlet

Como las condiciones son Dirichlet podemos tomar el espacio \[ V^h = \mathrm{span} \{ \varphi_1, \cdots, \varphi_M \}. \]

Si la malla es uniforme \(x_{i+1} - x_i = \Delta x\) y tomamos condiciones Dichlet \[ \int_{-1}^1 \frac{\partial \varphi_i^h}{\partial x} \frac{\partial \varphi_j^h}{\partial x} = \frac 1 {\Delta x} \begin{cases} 2 & i = j \\ {-1} & |i-j| = 1, \\ 0 & |i-j| > 1. \end{cases} \]

Acabamos con la famosa matriz para el problema de Dichilet \[ A^h = \frac{1}{\Delta x} \begin{pmatrix} 2 & - 1 & \\ -1 & 2 & -1 \\ && \ddots \\ && -1 & 2 \end{pmatrix} \]

No calcularemos analíticamente la matriz \(M^h\) ahora.

Observación. En el método de diferencias finitas la matriz tiene factor \((\Delta x)^2\). Pero el vector \[ b_i^h(t) = \int_{-1}^1 f(t,x) \varphi_i^h(x) \diff x \sim \Delta x \]

Matriz de rigidez para el Laplaciano con condiciones Neumann

Como las condiciones son Neumann podemos tomar el espacio \[ V^h = \mathrm{span} \{ \varphi_0, \varphi_1, \cdots, \varphi_M, \varphi_{M+1} \}. \]

Para los puntos interiores los la integrales son iguales. Para los puntos de borde \[ \int_{-1}^1 \frac{\partial \varphi_0}{\partial x}\frac{\partial \varphi_0}{\partial x} = \frac 1 {\Delta x} \qquad \int_{-1}^1 \frac{\partial \varphi_0}{\partial x}\frac{\partial \varphi_1}{\partial x} = \frac {-1} {\Delta x} \]

Se tiene la matriz \[ A^h = \frac{1}{\Delta x} \begin{pmatrix} 1 & - 1 & \\ -1 & 2 & -1 \\ && \ddots \\ && -1 & 1 \end{pmatrix} \]

Cálculo simbólico

El paquete SymPy.jl

Utilizaremos una librería de Python sympy. Debemos tener instalado python en el sistema

En MacOS / Linux ejecutar en terminal

$ pip3 install sympy

Y luego en la interfaz en Julia

using Pkg; Pkg.add("SymPy")

En Windows si Python está instalado mediante Anaconda, el cálculo simbólico explicado abajo sólo funcionará en libreta de Jupyter. En una libreta de Jupyter escribir

using Pkg; Pkg.add("SymPy")

Función lineal a trozos para puntos equiespaciados

Son un desplazamiento de

using SymPy, Plots, LaTeXStrings

x  = symbols("x", real=true)

Φ = sympy.Piecewise( 
        ( 1 + x , And(Ge(x,-1),Le(x, 0))) , 
        ( 1 - x, And(Le(x,1),Gt(x, 0))), 
        (0, True))

display(Φ)

plot(Φ,xlim=(-2.0, 2.0), label=L"\Phi", thickness_scaling=1.5, linewidth=2)
⎧x + 1  for x ≥ -1 ∧ x ≤ 0
⎪                         
⎨1 - x  for x ≤ 1 ∧ x > 0 
⎪                         
⎩  0        otherwise

Problema con condiciones Dirichlet

Funciones de base lineales a trozos con condiciones de Dirichlet

Tomamos la malla \[ x_0 = -1 < x_1 = -1+\Delta x < \cdots < x_{M+1} = 1 \] Las condiciones de Dirichlet homogéneas nos dicen que los valores \(u^h(t,\pm 1)=0\), luego las incógnitas son \(u_1^h, \cdots, u_N^h\).

Construímos las funciones de base lineales a trozos en los puntos que son incógnitas

Δx = 2.0^(-2)
xx = (-1.0 + Δx):Δx:(1.0-Δx)
NN = length(xx);

φ = subs.(Φ , x ,  x/Δx .- xx/Δx )

display(φ[1]);
⎧4.0⋅x + 4.0   for 4.0⋅x + 3.0 ≥ -1 ∧ 4.0⋅x + 3.0 ≤ 0
⎪                                                    
⎨-4.0⋅x - 2.0  for 4.0⋅x + 3.0 ≤ 1 ∧ 4.0⋅x + 3.0 > 0 
⎪                                                    
⎩     0                      otherwise               
p = plot(xlim=(-1,1),thickness_scaling=1.5,legend=:outerright,
    xlabel=L"x",ylabel=L"u")
for i=1:length(xx)
    p=plot!(φ[i], label=latexstring("\\varphi_{$(i)}"),linewidth=2)
end
plot(p)

Matriz de masa con condiciones de Dirichlet

Con esta elección, el cálculo de las matrices \[ M_{ij} = \int_a^b \varphi_j \varphi_i , \qquad A_{ij} = \int_a^b \frac{\diff \varphi_j}{\diff x} \frac{\diff \varphi_i}{\diff x} \] puede hacerse de forma analítica.

Primero la matriz de masa

M = zeros(NN,NN)
for i=1:NN, j=maximum([i-1,1]):minimum([i+1,NN])
        M[j,i] = integrate(φ[i]*φ[j],(x,-1.0,1.0))
end
M
7×7 Matrix{Float64}:
 0.166667   0.0416667  0.0        0.0        0.0        0.0        0.0
 0.0416667  0.166667   0.0416667  0.0        0.0        0.0        0.0
 0.0        0.0416667  0.166667   0.0416667  0.0        0.0        0.0
 0.0        0.0        0.0416667  0.166667   0.0416667  0.0        0.0
 0.0        0.0        0.0        0.0416667  0.166667   0.0416667  0.0
 0.0        0.0        0.0        0.0        0.0416667  0.166667   0.0416667
 0.0        0.0        0.0        0.0        0.0        0.0416667  0.166667

Matriz de rigidez con condiciones de Dirichlet

Para construir la matriz de rigidez calculamos primero las derivadas simbólicas

Dφ = diff.(φ,x)
display(Dφ[1])
⎧4.0   for 4.0⋅x + 3.0 ≥ -1 ∧ 4.0⋅x + 3.0 ≤ 0
⎪                                            
⎨-4.0  for 4.0⋅x + 3.0 ≤ 1 ∧ 4.0⋅x + 3.0 > 0 
⎪                                            
⎩ 0                  otherwise

Calculamos la matriz de rigidez

A = zeros(NN,NN)
for i=1:NN, j=maximum([i-1,1]):minimum([i+1,NN])
    A[j,i] = integrate(Dφ[i]*Dφ[j],(x,-1.0,1.0))
end
A
7×7 Matrix{Float64}:
  8.0  -4.0   0.0   0.0   0.0   0.0   0.0
 -4.0   8.0  -4.0   0.0   0.0   0.0   0.0
  0.0  -4.0   8.0  -4.0   0.0   0.0   0.0
  0.0   0.0  -4.0   8.0  -4.0   0.0   0.0
  0.0   0.0   0.0  -4.0   8.0  -4.0   0.0
  0.0   0.0   0.0   0.0  -4.0   8.0  -4.0
  0.0   0.0   0.0   0.0   0.0  -4.0   8.0

Resolución de la ecuación de Laplace con condición Dirichlet

Resolvemos la ecuación de Laplace \[ \begin{cases} -\frac{\partial^2 u}{\partial x^2} = 12 x^2 & \text{en } (-1,1) \\ u(1) = u(-1) = 0 \end{cases} \] que tiene solución exacta \(u(x) = 1-x^4\).

u = 1 - x^4
f = 12*x^2

b = zeros(NN,1)
for i=1:NN
    b[i] = integrate(f*φ[i],(x,-1.0,1.0))
end

uh = A \ b;

plot([-1;xx;1],[0;uh;0], label="Sol numérica", thickness_scaling=1.5,linewidth=2,legend=:outerright)
plot!(u, label="Sol exacta",linewidth=2)

Problemas con condición Neumann

Funciones de base lineales a trozos en el caso Neumann

Con la malla \[ x_0 = -1 < x_1 = -1+\Delta x < \cdots < x_{M+1} = 1 \]

Ahora debemos tomar el espacio \[ V^h = \mathrm{span} \{ \varphi_0, \cdots, \varphi_{M+1} \} \]

Es decir

Code 
Δx = 2.0^(-2)
xx = -1.0:Δx:1.0

φ = subs.(Φ , x ,  x/Δx .- xx/Δx )

p = plot(xlim=(-1,1),thickness_scaling=1.5,legend=:outerright,
    xlabel=L"x",ylabel=L"u")
for i=1:length(xx)
    p=plot!(φ[i], label=latexstring("\\varphi_{$(i-1)}"),linewidth=2)
end
plot(p)