Elementos finitos

Queremos resolver una ecuación del tipo: $$ -u'' + Vu = f $$ o, escribiendo la variable $x$ $$ -u''(x) + V(x)u(x) = f(x) $$ en un intervalo $[a,b]$, con condiciones de frontera $u(a)=u(b)=0$.

Fijamos una serie de nodos $x_0=a,x_1,\dots,x_{n}=b$. Buscamos una solución aproximada en el espacio $V$ de funciones lineales en cada intervalo $[x_i,x_{i+1}]$, que además valen 0 en los extremos: $u(x_0)=u(x_n)=0$.

En este espacio, la ecuación anterior es cierta si y sólo si el producto de $-u''(x) + V(x)u(x) - f(x)$ por cualquier función $v\in V$ tiene integral 0: $$ \int_a^b \left(-u''(x) + V(x)u(x) - f(x) \right) v(x) = 0\;\text{ para cualquier }v\in V $$

Esa integral se puede simplificar un poco integrando por partes (usando que $v(x_0)=v(x_n)=0$). $$ \int_a^b u'(x)v'(x) + V(x)u(x)v(x) - f(x) v(x) = 0\;\text{ para cualquier }v\in V $$ y vemos que esta integral es exactamente la derivada del funcional $E$ en el "punto" $u$ y en la dirección $v$, donde $E$ es un funcional de energía definido en el espacio $V$ (a cada función de $V$ le asigna un número real): $$ E(u)=\int_a^b (u'(x))^2 + V(x)u(x)^2 - f(x) u(x) = 0. $$

Cualquier solución de la ecuación diferencial es un punto crítico del funcional de energía, y viceversa.

$$ -u'' + Vu = f\; \text{ en sentido "débil", en }V\Leftrightarrow \nabla_u E(v)=0 \text{ para cualquier }v\in V\Leftrightarrow \nabla_u E=0 $$

Cualquier función de $V$ se puede escribir en la forma $$ u(x) = \sum_{i=1}^n c_i\phi_i(x) $$ Donde $\phi_i$ son funciones que forman una base de $V$. La función $\phi_i$ es la única función de $V$ que vale $1$ en $x_i$, y $0$ en todos los demás nodos

finite_element_d1.png finite_element_d1_bis.png

In [2]:
function f=FE(x1,x2,x3)
% Return a piecewise linear "basis" function

% piecewise inline function
% https://stackoverflow.com/questions/796072/how-can-i-create-a-piecewise-inline-function-in-matlab
%Heaviside function
    H = @(x) (x >= 0);
%H(x-x1).*H(x2-x) vale 1 si x esta en el intervalo [x1, x2], y 0 en otro caso...
    f = @(x)(H(x-x1).*H(x2-x).*(x-x1)./(x2-x1) + H(x-x2).*H(x3-x).*(x3-x)./(x3-x2));
end
In [4]:
phi = FE(-0.25, 0, 0.25);
x_eval = linspace(-1,1,200);
y_eval = phi(x_eval);
plot(x_eval, y_eval);

Si desarrollamos el funcional de energía en esta base tenemos: $$ \begin{array}{rcll} E(u)&=&\int_a^b &(u'(x))^2 + V(x)u(x)^2 - f(x) u(x) = 0\\ E\left(\sum_{i=1}^n c_i\phi_i(x)\right)&=&\Large\int_a^b& \left(\sum_{i=1}^n c_i\phi_i'(x)\sum_{j=1}^n c_j\phi_j'(x)\right) \\ &&&+ V(x)\left(\sum_{i=1}^n c_i\phi_i(x)\sum_{j=1}^n c_j\phi_j(x)\right)\\ &&&- f(x)\left(\sum_{i=1}^n c_i\phi_i(x)\right) = 0 \end{array} $$ $$ \begin{array}{rcl} E\left(\sum_{i=1}^n c_i\phi_i(x)\right)&=& \left(\sum_{i=1}^n \sum_{j=1}^n c_jc_i\int_a^b\phi_i'(x)\phi_j'(x)\right) \\ &&+\left(\sum_{i=1}^n \sum_{j=1}^n c_jc_i\int_a^b V(x)\phi_i(x)\phi_j(x)\right)\\ &&- \sum_{i=1}^n c_i\int_a^b f(x)\phi_i(x) = 0 \end{array} $$

Definimos:

  • Una matriz $A_1$ cuya entrada $i,j$ es $\int_a^b\phi_i'(x)\phi_j'(x)\:dx$
  • Una matriz $A_2$ cuya entrada $i,j$ es $\int_a^b V(x)\phi_i(x)\phi_j(x)\:dx$
  • La matriz $A$ es $A=A_1+A_2$.
  • Un vector $\mathbf{b}$ cuya entrada $i$ es $\int_a^b f(x)\phi_i(x)\:dx$
  • Un vector $\mathbf{c}$ cuya entrada $i$ es la incógnita $c_i$. y el funcional $E$ toma la forma: $$ \begin{array}{rcll} E(u)=E\left(\sum_{i=1}^n c_i\phi_i(x)\right)=E(\mathbf{c})&=&\mathbf{c}^T\cdot A\cdot \mathbf{c} - \mathbf{c}^T\cdot\mathbf{b}\\ \end{array} $$

Es importante observar que las integrales que definen $A$ se anulan si j no es i, i-1 ó i+1, porque $\phi_{i}(x)$ se anula si $x$ no está en $[x_{i-1}, x_{i+1}]$, y $\phi_{j}(x)$ se anula si $x$ no está en $[x_{j-1}, x_{j+1}]$, luego el producto $\phi_{i}(x)\phi_{j}(x)$ se anula si $x$ no está en la intersección de los dos intervalos $[x_{i-1}, x_{i+1}]$ y $[x_{j-1}, x_{j+1}]$.

Es decir, la matriz $A$ es tridiagonal.

El funcional $E$ es cuadrático y definido positivo, y por tanto tiene un mínimo absoluto en su único punto crítico, que se obtiene derivando: $$ \nabla_{\mathbf{c}} E = A\cdot \mathbf{c} - \mathbf{b} = \mathbf{0} $$ y finalmente encontramos una solución de nuestro problema variacional resolviendo un sistema de ecuaciones lineales.

In [38]:
l = 120;
Q = 100;
E = 3e7;
S = 1e3;
I = 625;

%Esta linea define x0, xf (números) y x_span (una lista con esos dos mismos números);
x0 = 0;
xf = l;
x_span = [0, l];
alpha = 0;
beta = 0;

% -y''(x) + V(x) y(x) = f (x);
V = @(x)(S/(E*I));
f = @(x)((Q/(2*E*I)).*x.*(x-l));

% Probamos distintos valores;
n = 20;

% Mallado;
h = (xf - x0)/n;
xs = linspace(x0,xf,n+1);
% puntos interiores;
xs_int = xs(2:n);

% Matriz banda A;
A = zeros(n-1,n-1);
% Vector b;
b = zeros(n-1,1);
for i = 1:n-1;
    xi = xs_int(i);
    xi_menos = xi - h;
    xi_mas   = xi + h;
    phi_i = FE(xi_menos, xi, xi_mas);
    %El ultimo argumento al llamar a quad es la tolerancia
    %Puede parecer que es una tolerancia innecesariamente baja
    % dada la escasa precision que estamos usando para los elementos
    % finitos, pero quad de matlab calcula las integrales bastante
    % mal si ponemos una tolerancia tan baja como 1e-7 :-O !
    %Seguramente se debe a que Simpson adaptativo puede
    % basarse en estimaciones del error muy malas cuando el intervalo
    % es grande, como vimos al estudiar la integracion adaptativa.
    %Es importante recordar que Simpson adaptativo nos proporciona
    % estimaciones del error, pero no cotas garantizadas.
    %En octave no es necesario exigir esa tolerancia para tener un
    % resultado decente, porque quad de octave usa otro metodo numerico.
    A(i,i) = 2/h + quad(@(x)(V(x).*phi_i(x).^2), xi_menos,xi_mas, 1e-10);
    b(i) = -quad(@(x)(f(x).*phi_i(x)), xi_menos,xi_mas, 1e-10);
    if i>1;
        phi_i_menos = FE(xi-2*h, xi-h, xi);
        A(i,i-1) = -1/h + quad(@(x)(V(x).*phi_i(x).*phi_i_menos(x)), xi_menos,xi, 1e-10);
    end
    if i<n-1;
        phi_i_mas = FE(xi, xi+h, xi+2*h);
        A(i,i+1) = -1/h + quad(@(x)(V(x).*phi_i(x).*phi_i_mas(x)), xi,xi_mas, 1e-10);
    end
end

% Solucion;
y_int = A\b;
ys_FEM = zeros(n+1,1);
ys_FEM(2:n) = y_int;

Recordemos que nuestra aproximación a la solución es $$ u(x)=\sum_{i=1}^n c_i\phi_i(x) $$ pero como $u$ es lineal a trozos, podemos hacer la gráfica de u usando el comando plot, que interpola linealmente. Los valores de $u$ en los nodos son exactamente: $$ u(x_j)=\sum_{i=1}^n c_i\phi_i(x_j) = c_j $$ porque $\phi_i(x_j)$ vale 0 cuando $i\neq j$ y 1 cuando $i=j$.

In [39]:
%Ponemos signo negativo al desplazamiento, ya que es más intuitivo 
%graficar el desplazamiento hacia abajo
hold on;
plot(xs,-ys_FEM);
xlabel('x')
ylabel('desplazamiento vertical')
title('Ecuacion de segundo orden lineal, con elementos finitos')
%Se ve a ojo que el máximo se alcanza en el punto medio
%También se puede deducir de la simetría del problema
%Pero este código permite encontrar el máximo de ys 
%y el punto donde se alcanza el máximo en general
[ymax, argmax] = max(ys_FEM);
xmax = xs(argmax);
xmax, ymax
plot([xmax], [-ymax], 'r*')
legend('desplazamiento', sprintf('desplazamiento maximo %.3f', ymax))
hold off;

xmax =  60
ymax =  0.014399

Hacer las integrales con quad no es eficiente, porque vamos a hacer los cálculos con un error mucho mayor.

Podemos sustituir esas integrales con una regla de Newton-Cotes (como la regla de Simpson) evaluando en los nodos $x_i$, pero los buenos resultados no son muy buenos. Por ejemplo, comprobamos que las aproximaciones por Simpson para $\int_a^b f(x)\phi_i(x)$ y para $\int_a^b f(x)\phi_i(x)^2$ son exactamente iguales, a pesar de que la segunda integral es claramente inferior: $$ \int_a^b f(x)\phi_i(x) = \int_{x_{i-1}}^{x^{i+1}} f(x)\phi_i(x) \approx (h/6)\left(f(x_{i-1})\phi(x_{i-1}) + 4 f(x_i)\phi(x_i) + f(x_{i+1})\phi(x_{i+1}) \right) =(4h/6)f(x_i) $$ $$ \int_a^b f(x)\phi_i(x)^2 = \int_{x_{i-1}}^{x^{i+1}} f(x)\phi_i(x)^2 \approx (h/6)\left(f(x_{i-1})\phi(x_{i-1})^2 + 4 f(x_i)\phi(x_i)^2 + f(x_{i+1})\phi(x_{i+1})^2 \right) =(4h/6)f(x_i) $$

In [41]:
xs = linspace(-1,1,300);
x0 = 0;
h = 0.25;
phi = FE(x0-h,x0,x0+h);
integrando1 = @(x)(V(x).*phi(x));
integrando2 = @(x)(V(x).*phi(x).^2);
ys1 = integrando1(xs);
ys2 = integrando2(xs);
hold on;
plot(xs, ys1)
plot(xs, ys2)
title('V\cdot\phi versus V\cdot\phi^2')
legend('V\cdot\phi', 'V\cdot\phi^2')
hold off;

Una idea mejor es sustituir $f$ por su aproximación en el espacio $V$ (funciones lineales en cada intervalo $[x_i,x_{i+1}]$): $$ f(x)\approx\sum_{i=1}^n f(x_i)\phi_i(x) $$ Las funciones $f(x)$ y $\sum_{i=1}^n f(x_i)\phi_i(x)$ toman los mismos valores en todos los nodos $x_i$.

Hacemos las integrales: $$ \int_a^b f(x)\phi_j(x) \approx \int_a^b \sum_{i=1}^n f(x_i)\phi_i(x) \phi_j(x) = \sum_{i=1}^n f(x_i)\int_a^b \phi_i(x) \phi_j(x) $$ Las integrales $\int_a^b \phi_i(x) \phi_j(x)$ son fáciles de hacer a mano (y recordamos que sólo tres de ellas son distintas de cero). Resumimos los resultados:

  • $\int_a^b f(x)\phi_j(x)\approx (h/6)\left(f(x_{j-1}) + 4 f(x_j) + f(x_{j+1})\right)$
  • $\int_a^b \phi_j'(x)^2\approx (2/h)$
  • $\int_a^b \phi_j'(x)\phi_{j+1}'(x)\approx (-1/h)$
  • $\int_a^b V(x)\phi_j(x)^2\approx (h/12)\left(V(x_{j-1}) + 6 V(x_j) + V(x_{j+1})\right)$
  • $\int_a^b V(x)\phi_j(x)\phi_{j+1}(x)\approx (h/12)\left(V(x_j) + V(x_{j+1})\right)$

Ejercicio

Reescribe el método de los elementos finitos sustituyendo todas las integrales con quad por las aproximaciones anteriores.

In [ ]:
In [ ]:
In [ ]:
In [ ]: