θ (generic function with 1 method)
Métodos Numéricos Avanzados. 2023-2024
Hasta ahora, hemos trabajado con simplexes (segmentos, triángulos, …) y con funciones lineales en cada simplex.
La malla podría no de simplexes. Por ejemplo, en dimensión 2, podría ser de cuadrados.
Pero si \(K\) es un cuadrado de vértices \(\mathbf x_i\), no podemos encontrar una función lineal tal que \(\varphi(\mathbf x_j) = \delta_{1j}\).
Pero podemos encontrar \(\varphi\) polinomios de orden superior.
Vamos a tomar un elemento “normalizado”, y a utilizarlo para construir las funciones de base. Digamos que la malla \(\mathcal T\) se compone de subconjunto \(K_i \in \mathbb T\), todos ellos transformaciones de uno de conjunto de referencia \(K_a\) (si hay varios \(K_a\) se habla de mallas mixtas).
La idea es construir un conjunto de polinomios \(P\) definidos \(K_a\), que se pueda “pegar” para construir una base \[ V_h = \{ v \in V : v|_{K_i} \circ \mathbf x_{K_i} \in P \} \] donde \(\mathbf x_{K_i} : K_a \to K_i\).
En [Ciarlet, 2002] se hace estudia la noción rigurosa de “elemento”, y los elementos válidos y no válidos.
En esta sección buscamos discretizar \[ V_h = \{ v \in H^1(\Omega) : v|_{K_i} \circ \mathbf x_{K_i} \in P \} \]
donde \(P\) son polinomios de grado \(p\).
Dado que en el interior de cada \(K_i\) tenemos polinomios, basta pegar de manera continua en las aristas.
Tomamos un elemento de referencia \(K_a = (-1,1)\). Vamos a construir polinomios de orden \(p\) que se anulen en ciertos puntos \[ X_0 = -1, \qquad X_1 = 1, \qquad X_{j} = -1 + \frac{2(j-1)}{p} \text{para} j = 2, \cdots, p \]
Así \[ \theta_i (\xi) = \frac{ \xi - X_1 }{X_i - X_1} \frac{ \xi - X_2 }{X_i - X_2} \cdots \frac{\xi - X_{i-1}} {X_i - X_{i-1}} \frac{\xi - X_{i+1}} {X_i - X_{i+1}} \cdots \frac{\xi - X_{p+1}} {X_i - X_{p+1}} \]
La idea de las familias jerárquicas es encontrar bases \(\mathcal B_p \subset \mathcal B_{p+1}\), donde \(\mathcal B_p\) es una base de \(\mathbb R_p [X]\).
Para buscar bases de orden superior, se toman los polinomios de Legendre:
\[ L_n (x) = \frac{(-1)^n}{2^n n!} \frac{d^n}{dx^n} \left[ (1-x)^n(1+x)^n \right] \]
Tienen la propiedad
\[ \begin{aligned} L_0(x) &= 1, L_1(x) = x, \\ L_k(x) &= \frac{2k-1}{k} x L_{k-1} (x) - \frac{k-1}{k} L_{k-2}(x), \forall k \ge 2. \end{aligned} \]
y los polinomios de Lobatto
\[ \begin{aligned} l_0(x) &= \frac{1-x}{2}, \\ l_1(x) &= \frac{x+1}{2}, \\ l_k(x) &= \sqrt{\frac{2k-1}{2}} \int_{-1}^x L_{k-1} (\xi) d \xi. \end{aligned} \]
Las funciones \(l_0, \cdots, l_p\) forman una base de \(\mathbb R_p[X]\).
Además \(l_k(-1) = l_k(1) = 0\) para \(k \ge 2\).
En este contexto \(\theta_k = l_k\) para \(k = 0, \cdots, p\).
Fijamos un intervalo \((a,b)\) y una malla \(K_i = (x_{i}, x_{i+1})\) con \(i=1, \cdots I\) donde es construir la aplicación afín \(\mathbf x_i : (-1,1) \to (x_{i}, x_{i+1})\).
Tomemos \(p\) fijo y una elección de \(\theta_i\) una de las elecciones anteriores (nodal o jerárquica). Observamos que \[ \begin{aligned} \theta_0 (-1) = 1 &, \qquad \theta_0(1) = 0, \\ \theta_1 (-1) = 0 &, \qquad \theta_1(1) = 1, \\ \theta_k(-1) = 0 &, \qquad \theta_k(1) = 0. \end{aligned} \]
Para los problemas elípticos que hemos visto, queremos que la base para Galerkin sea de funciones continuas tenemos que pegar con continuidad.
Esto nos permite construir funciones base de tipo vértice (con \(v_{i,0} (x_j) = \delta_{ij}\)) \[ v_{i,0}(x) = \begin{cases} (\theta_1 \circ \mathbf x_{K_i}^{-1}) (x) , & x \in K_i \\ (\theta_0 \circ \mathbf x_{K_{i+1}}^{-1}) (x) , & x \in K_{i+1} \\ 0, & x \notin K_{i} \cup K_{i+1}. \end{cases} \]
Ejemplo. Tomando como \(\theta_i\) las funciones nodales
mesh = [0.1, 0.2, 0.4]
p = 3
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
function v(x::Real)
if x >= mesh[1] && x < mesh[2]
ξ = -1 + 2*(x-mesh[1])/(mesh[2] - mesh[1])
return θ(1,ξ,X)
elseif x >= mesh[2] && x < mesh[3]
ξ = -1 + 2*(x-mesh[2])/(mesh[3] - mesh[2])
return θ(0,ξ,X)
else
return 0.0
end
end
plot(x -> v(x), label=L"v", xlim=(0,0.5), xlabel=L"x")
scatter!(mesh, zeros(size(mesh)), label="")
También podemos construir funciones base de tipo burbuja, para \(k \ge 2\)
\[ v_{i,k-1}(x) = \begin{cases} \theta_{k} \circ \mathbf x_{K_i}^{-1} & x \in K_i \\ 0 & x \notin K_i \end{cases} \]
Estas no existen para \(p = 1\).
Ejemplo. Tomando como \(\theta_i\) las funciones nodales
mesh = [0.1, 0.2, 0.4]
p = 3
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
function w1(x::Real)
if x >= mesh[1] && x < mesh[2]
ξ = -1 + 2*(x-mesh[1])/(mesh[2] - mesh[1])
return θ(2,ξ,X)
else
return 0.0
end
end
function w2(x::Real)
if x >= mesh[1] && x < mesh[2]
ξ = -1 + 2*(x-mesh[1])/(mesh[2] - mesh[1])
return θ(3,ξ,X)
else
return 0.0
end
end
plot(x -> w1(x), label=L"w_1", xlim=(0,0.5))
plot!(x -> w2(x), label=L"w_2", xlim=(0,0.5))
scatter!(mesh, zeros(size(mesh)), label="")
xlabel!(L"x")
Con estas elecciones tomamos como espacio de funciones base para la aproximación de Galerkin \[ V_h = \mathrm{span} \{ v_{i,k} : i = 1, \cdots, I \text{ y } k = 0, \cdots, p \}. \]
Este mismo planteamiento se puede hacer en dimensión superior, cuando \(K\) no es un triángulo/tetrahedro: [Šolin, Segeth & Doležel, 2004]
SymPy.jl
using SymPy, Plots, LaTeXStrings
Δx = 1.0
xx = -1.0:Δx:1.0
ξ = symbols("ξ", real=true)
x = symbols("x", real=true)
p = 3
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
θ(i,ξ,X) = prod([
(ξ - X[j])/(X[i+1] - X[j])
for j in setdiff(1:length(X),i+1)
] )
v(i,x) = sympy.Piecewise(
( θ(1,-1 + 2*(x-xx[i-1])/Δx, X), And(Ge(x,xx[i-1]),Le(x, xx[i]))) ,
( θ(0,-1 + 2*(x-xx[i])/Δx, X), And(Ge(x,xx[i]),Le(x, xx[i+1]))) ,
(0, True))
w1(i,x) = sympy.Piecewise(
( θ(2,-1 + 2*(x-xx[i-1])/Δx, X), And(Ge(x,xx[i-1]),Le(x, xx[i]))) ,
(0, True))
w2(i,x) = sympy.Piecewise(
( θ(3,-1 + 2*(x-xx[i-1])/Δx, X), And(Ge(x,xx[i-1]),Le(x, xx[i]))) ,
(0, True))
φ = [ [ v(i,x) for i in 2:length(xx)-1 ];
[ w1(i,x) for i in 2:length(xx) ] ;
[ w2(i,x) for i in 2:length(xx) ] ;
];
plot(legend_position=:outerright,xlim=(-1.01,1.01))
scatter!(xx,zeros(length(xx)),label="")
for i in 1:length(φ)
plot!(φ[i],label=latexstring("\\varphi_{$i}"))
end
plot!()
xlabel!(L"x")
SymPy.jl
Dφ = diff.(φ,x)
NN = length(φ)
A = zeros(NN,NN)
for i=1:NN, j=1:NN
A[i,j] = integrate(Dφ[j]*Dφ[i],(x,-1.0,1.0))
end
A
u_analytic = 1 - x^4
f_analytic = 12*x^2
b = [integrate(f_analytic*φ[i],(x,-1.0,1.0)) for i=1:NN]
uh_coord = A \ b;
uh(x) = sum([uh_coord[i]*φ[i](x) for i=1:length(φ)])
plot(x->uh(x), label="Sol numérica", thickness_scaling=1.5,linewidth=2,legend=:outerright,xlim=(-1.0001,1.0001))
plot!(u_analytic, label="Sol exacta",linewidth=2)
scatter!(xx,zeros(length(xx)),label="Malla")
xlabel!(L"x")
ylabel!(L"u")
Las construcciones anteriores son solo continuas (no \(C^1\)) en los puntos de la malla. Esto es suficiente para problemas de segundo orden como los que hemos visto: \(-\frac{\partial^2 u}{\partial x^2} = f\).
Es muy habitual, por ejemplo en elasticidad, encontrar problemas de cuarto ordern \(\frac{\partial^4 u}{\partial x^4} = f\). La formulación débil es, entonces (si \(I = (-1,1)\)) \[ \int_{-1}^1 \frac{\partial^2 u}{\partial x^2} \frac{\partial^2 \varphi}{\partial x^2} = \int_{-1}^1 f \varphi + \text{ términos en $-1$ y $1$ (condiciones de contorno)} \]
En este caso buscaremos discretizar \[ V_h = \{ v \in H^2(\Omega) : v|_{K_i} \circ \mathbf x_{K_i} \in P \} \]
Debemos “pegar” de manera \(C^1\) en las aristas de la malla.
Para “pegar” de forma \(C^1\) debemos tener cuidado con los valores de \(\varphi'\) en los puntos de la malla.
Si trabajamos en el intervalo de referencia \(K = [0,1]\) esto significa tener cuidado con \(\varphi'(0), \varphi'(1)\).
Podemos construir primero dos polinomios de base cuyas derivadas se anulen en los extremos ellas o sus derivadas \[ \theta_0(0) = 1, \qquad \theta_0(1) = 0, \qquad \theta_0'(0) = 0, \qquad \theta_0'(1) = 0 \]
\[ \theta_1(0) = 0, \qquad \theta_1(1) = 1, \qquad \theta_1'(0) = 0, \qquad \theta_1'(1) = 0 \]
\[ \theta_2(0) = 0, \qquad \theta_2(1) = 0, \qquad \theta_2'(0) =1, \qquad \theta_2'(1) = 0 \]
\[ \theta_3(0) = 0, \qquad \theta_3(1) = 0, \qquad \theta_3'(0) =0, \qquad \theta_3'(1) = 1 \]
Esto requiere, al menos, grado 3, y en tal caso
Ahora pegamos haciendo \(v_{i}(x_j) = \delta_{ij}, v'_i (x_j) = 0\) y \(w_i (x_j) = 0, w_i'(x_j) = \delta_{ij}\) \[ v_{i} (x) = \begin{cases} \theta_1 \left( \frac{x-x_{i-1}}{x_{i}-x_{i-1}} \right) , & x \in [x_{i-1},x_{i}) \\ \theta_0 (\frac{x-x_{i}}{x_{i+1}-x_{i}}) , & x \in x \in [x_{i},x_{i+1})\\ 0, & \text{otherwise}. \end{cases} \qquad w_{i} (x) = \begin{cases} (x_{i}-x_{i-1}) \theta_4 \left( \frac{x-x_{i-1}}{x_{i}-x_{i-1}} \right) , & x \in [x_{i-1},x_{i}) \\ (x_{i+1}-x_{i}) \theta_3 (\frac{x-x_{i}}{x_{i+1}-x_{i}}) , & x \in x \in [x_{i},x_{i+1})\\ 0, & \text{otherwise}. \end{cases} \]
mesh = [-Inf, 0.0, 0.1, 0.2, 0.4, 0.5, +Inf]
function v(i::Integer,x::Real)
if x >= mesh[i-1] && x < mesh[i]
ξ = (x-mesh[i-1])/(mesh[i] - mesh[i-1])
return θHermite[2](ξ)
elseif x >= mesh[i] && x < mesh[i+1]
ξ = (x-mesh[i])/(mesh[i+1] - mesh[i])
return θHermite[1](ξ)
else
return 0.0
end
end
scatter(mesh[2:end-1], zeros(size(mesh[2:end-1])),
label="",
leg=:outerright)
for i=2:length(mesh)-1
plot!(x -> v(i,x), label=latexstring("v_{$(i-1)}"))
end
xlabel!(L"x")
mesh = [-Inf, 0.0, 0.1, 0.2, 0.4, 0.5, +Inf]
function v(i::Integer,x::Real)
if x >= mesh[i-1] && x < mesh[i]
dx = mesh[i] - mesh[i-1]
ξ = (x-mesh[i-1])/dx
return dx*θHermite[4](ξ)
elseif x >= mesh[i] && x < mesh[i+1]
dx = mesh[i+1] - mesh[i]
ξ = (x-mesh[i])/dx
return dx*θHermite[3](ξ)
else
return 0.0
end
end
scatter(mesh[2:end-1], zeros(size(mesh[2:end-1])),
label="",
leg=:outerright)
for i=2:length(mesh)-1
plot!(x -> v(i,x), label=latexstring("w_{$(i-1)}"))
end
xlabel!(L"x")
Tomando \(K = [0,1] \times [0,1]\) podemos construir \[ \theta_{ij}(x,y) = \theta_i(x) \theta_j(y), \qquad i,j \in \{0,\cdots,3\} \]
También se pueden construir elementos de tipo \(C^1\) en el triángulo de referencia, e.g., triángulo de Argyris.
Para pegar estos elementos en una malla hay que tener cuidado con rescalar las orientación de los bordes (ver [Domínguez & Sayas, 2008], [Valdman, 2020]).
Más información [Ciarlet & Lions, 1990, Capítulo VII]
Métodos Numéricos Avanzados. David Gómez-Castro