Ecuación de ondas

Resumen: En este ejemplo se considera un problema relativo a la ecuación de ondas.

Problema considerado

$\Omega = B\big((0,0);R\big) \setminus B\big( (0,0); r \big) \subset \mathbb{R}^2$ $T>0$ $u:\Omega\times(0,T) \to \mathbb{R}$ tal que

\begin{matrix} (1) & {\begin{cases} \partial_{t t} u + a \partial_{t} u - c^{2} Δ u + m u = f & en Ω \times (0, T), \\ u = 0 & sobre \partial Ω \times (0, T), \\ u (x, y, 0) = u_{0} (x, y) & en Ω, \\ \partial_{t} u (x, y, 0) = u_{1} (x, y) & en Ω, \end{cases} \end{matrix}

$a,c\in\mathbb{R}$ $a>0$ $m: \Omega \to \mathbb{R}$ $f:\Omega\times(0,T) \to \mathbb{R}$ $u_0, u_1: \Omega\to \mathbb{R}$ son funciones dadas.

Para calcular una solución aproximada, procedemos como sigue:

$t$ $t$ $\Omega$ .
Cada problema estacionario se discretiza mediante el método de los elementos finitos sobre una triangulación dada, como en los ejemplos de problemas elípticos considerados con anterioridad.

Discretización en tiempo

$[0,T]$ $M+1$ $0=t_0 <\dots < t_M=T$ $[0,T]$ $M$ $\Delta t$ :

\begin{matrix} Δ t = \frac{T}{M}, t_{n} = n Δ t, n = 0, 1, \dots, M . \end{matrix}

$n = 1,\dots, M$ $\eqref{eq.wave2D.1}$ $t = t_n$ :

\begin{matrix} (2) & {\begin{cases} \partial_{t t} u (x, y, t_{n}) + a \partial_{t} u (x, y, t_{n}) - c^{2} Δ u (x, y, t_{n}) + m (x, y) u (x, y, t_{n}) = f (x, y, t_{n}) & en Ω, \\ u (x, y, t_{n}) = 0 sobre \partial Ω . \end{cases} \end{matrix}

Denotamos

\begin{matrix} u^{n} (x, y) : = u (x, y, t_{n}), f^{n} (x, y) : = f (x, y, t_{n}), n = 0, 1, \dots, M . \end{matrix}

$n = 0$ , se tiene

\begin{matrix} u^{0} (x, y) = u (x, y, t_{0}) = u (x, y, 0) = u_{0} (x, y) . \end{matrix}

$n=1$ , se tiene, utilizando un desarrollo de Taylor

\begin{matrix} u^{1} (x, y) = u (x, y, Δ t) \approx u (x, y, 0) + \partial_{t} u (x, y, 0) Δ t = u_{0} (x, y) + Δ t u_{1} (x, y) \end{matrix}

$n\geq 2$ $\eqref{eq.wave2D.2}$ $u$ $t$ por aproximaciones en diferencias:

\begin{matrix} \partial_{t t} u (x, y, t_{n}) \approx \frac{u (x, y, t_{n}) - 2 u (x, y, t_{n - 1}) + u (x, y, t_{n - 2})}{(Δ t)^{2}} \end{matrix}

\begin{matrix} \partial_{t} u (x, y, t_{n}) \approx \frac{u (x, y, t_{n}) - u (x, y, t_{n - 1})}{Δ t} \end{matrix}

se tiene la siguiente ecuación (aproximada)

\begin{matrix} \frac{u^{n} - 2 u^{n - 1} + u^{n - 2}}{(Δ t)^{2}} + a \frac{u^{n} - u^{n - 1}}{Δ t} - c^{2} Δ u^{n} + m u^{n} = f^{n} en Ω, \end{matrix}

O, equivalentemente,

\begin{matrix} \frac{1}{(Δ t)^{2}} u^{n} + \frac{a}{Δ t} u^{n} - c^{2} Δ u^{n} + m u^{n} = f^{n} - \frac{1}{(Δ t)^{2}} u^{n - 2} + (\frac{2}{(Δ t)^{2}} + \frac{a}{Δ t}) u^{n - 1} en Ω . \end{matrix}

$n= 2,\dots,M$ $u^n$ $u^{n-1}$ $u^{n-2}$ han sido calculadas en las etapas anteriores). Por tanto, el esquema discreto en tiempo es como sigue:

$u^0 = u_0$ .
$u^1 = u_0 + \Delta t u_1$ .
$n=2,\dots, M$ $u^n:\Omega \to \mathbb{R}$ solución de
$\begin{matrix} (3) & {\begin{cases} \frac{1}{(Δ t)^{2}} u^{n} + \frac{a}{Δ t} u^{n} - c^{2} Δ u^{n} + m u^{n} \\ = f^{n} - \frac{1}{(Δ t)^{2}} u^{n - 2} + (\frac{2}{(Δ t)^{2}} + \frac{a}{Δ t}) u^{n - 1} en Ω, \\ u^{n} = 0 sobre \partial Ω . \end{cases} \end{matrix}$

$n$ $\eqref{eq.wave2D.3}$ es un problema que no depende del tiempo, que podemos aproximar usando el método de los elementos finitos razonando como en los ejemplos anteriores.

Formulación débil o variacional

$\eqref{eq.wave2D.3}$ $H^1_0(\Omega)$ $v\in H_0^1(\Omega)$ $\Omega$ $v$ $\partial\Omega$ , obtenemos:

\begin{matrix} (PV) & {\begin{cases} Hallar u^{n} \in H_{0}^{1} (Ω) tal que \\ c^{2} \int_{Ω} \nabla u^{n} \cdot \nabla v d x d y + \int_{Ω} (\frac{1}{(Δ t)^{2}} + \frac{a}{Δ t} + m) u^{n} v d x d y \\ = \int_{Ω} (f^{n} - \frac{1}{(Δ t)^{2}} u^{n - 2} + (\frac{2}{(Δ t)^{2}} + \frac{a}{Δ t}) u^{n - 1}) v d x d y \\ \forall v \in H_{0}^{1} (Ω) . \end{cases} \end{matrix}

$\eqref{eq.wave2D.PV}$ $n=2,\dots,M$ , comparten la misma forma bilineal

\begin{matrix} a (u, v) = c^{2} \int_{Ω} \nabla u^{n} \cdot \nabla v d x d y + \int_{Ω} (\frac{1}{(Δ t)^{2}} + \frac{a}{Δ t} + m) u^{n} v d x d y \end{matrix}

y que solo difieren en la forma lineal. En la resolución numérica, esto se traduce en el hecho de que todos los sistemas lineales a resolver tienen la misma matriz, que, por lo tanto, solo habrá que construir y factorizar la primera vez.

Código


x
//
//  Problem Data
//
real RR = 5., rr = 1.;                // Outer and inner radius
real T = 15.;                         // Final time
real a = 0.1, c = 0.5;                // edp coefficients

func KGm = 0.1*(1.+x*x)*(1.+y*y);     // Klein-Gordon function
func force = 0.;                      // External force
//func u0 = x*x - y*y;                  // Initial u
func u0 = x*x + y*y - 25.;            // Initial u
func u1 = 0.;                         // Initial u_t

//
//  Data for mesh and time discretization
//
int nxo  = 200;                                  // Discret. of the outer border
int nxi  = 50;                                   // Discret. of the inner border
int nt   = 100;                                  // Time steps
int LO = 7;                                      // Label of outer border
int LI = 12;                                     // Label of inner border

//
//  Coefficients of stationary problems
//
real c2   = c*c;
real dt   = T/nt;
real dt1  = 1./dt;
real dt2  = dt1*dt1;
real adt0 = a*dt1;
real adt1 = dt2 + adt0;
real adt2 = 2*dt2 + adt0;

//
//  Mesh
//
border CO(t=0, 2*pi){x=RR*cos(t); y=RR*sin(t); label = LO;}    // outer
border CI(t=0, 2*pi){x=rr*cos(t); y=rr*sin(t); label = LI;}    // inner
plot( CO(nxo) + CI(-nxi), dim = 2, wait = true );
mesh Th = buildmesh( CO(nxo) + CI(-nxi) );
plot(Th, cmm = " Initial mesh ", wait = true);

//
// FE space 
//
fespace Vh(Th, P1);
Vh u, v;
Vh uold2 = u0;                   // uold2 --> u^{i-2}
Vh uold1 = u0 + dt*u1;           // uold1 --> u^{i-1}
Vh alpha = adt1 + KGm;           // u coefficient = adt1+m(x)
Vh f = force;
Vh F;                            // RHS; to be defined at each iteration

//
//  Linear stationary problems
//
problem StWave(u, v, init=1) = 
        int2d(Th)(alpha*u*v + c2*(dx(u)*dx(v)+dy(u)*dy(v)))
      - int2d(Th)(F*v)
    + on(LI, u = 0.) + on(LO, u = 0.);

//
//  Time iterates
//
plot(uold2, dim = 3, nbiso = 40, cmm = " Time t = "+ 0., 
            value = true, fill = true, wait = true);
plot(uold1, dim = 3, nbiso = 40, cmm = " Time t = "+ dt, 
        value = true, fill = true, prev = true, wait=true);
real t;
for (int k = 2; k <= nt; k ++)
{
  t = k*dt;
  F = f + adt2 * uold1 - dt2 * uold2;     // current RHS
  StWave;
  plot(u, dim = 3, nbiso = 40, cmm = " Time t = "+ t, 
          value = true, fill = true, prev = true, wait=true);
    uold2 = uold1;
    uold1 = u;
}

Ejercicios

$\eqref{eq.wave2D.1}$ , con las siguientes variantes:

Almacenando todas las iteraciones en un "array" y dibujando la solución al finalizar bucle.
Almacenando todas las iteraciones en un fichero. Escribir otro programa que lea los datos almacenados y dibuje con ellos la solución.
$f(x,y,t) = \sin(t)$ .

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla