Stability and dispersion of the problem of seismic wave propagation in 2-D using the generalized finite difference method

==Resumen==

Este artículo muestra la resolución del problema de propagación de ondas sísmicas en 2-D, mediante la utilización de esquemas explícitos en diferencias finitas generalizadas (GFD), lo que permite la utilización de mallas regulares e irregulares.

Puesto que se utiliza un método explícito, es necesario obtener la condicióń de estabilidad, lo que se ha realizado mediante un análisis de von Neumann. También se obtienen las relaciones de dispersión en la estrella de las velocidades de fase para las ondas P y S, así como las de las velocidades de grupo.

Dada la importancia que en la aplicación del método tiene el control sobre la irregularidad de la malla, se han definido unos índices de irregularidad para la estrella (IIS) y la malla (IIC), analizándose su relación con la dispersión y el paso de tiempo utilizado en los cálculos.

==Abstract==

This paper shows the solution to the problem of seismic wave propagation in 2-D using generalized finite difference (GFD) explicit schemes. Regular and irregular meshes can be used with this method.

As we are using an explicit method, it is necessary to obtain the stability condition by using the von Neumann analysis. We also obtained the star dispersion formulas for the phase velocities for the P and S waves, as well as the ones for the group velocities.

As the control over the irregularity in the mesh is very important in the application of this method, we have defined an index of irregularity for the star (IIS) and another for the cloud (IIC), analyzing its relationship with the dispersion and time step used in the calculations.

==Palabras clave==

Métodos sin malla ; Diferencias finitas generalizadas ; Ondas ; Estabilidad ; Dispersión ; Mallas irregulares ; Método explícito

==Keywords==

Meshless methods ; Generalized finite difference ; Waves ; Stability ; Dispersion ; Irregular meshes ; Explicit method

==1. Introducción==

Son muchos los trabajos publicados sobre propagación de ondas sísmicas, incluso centrándonos en los relativos a la utilización del método de Diferencias Finitas (FD) y dentro de ellos los que tratan aspectos relativos al tipo de malla, la estabilidad y dispersión, que son los aspectos en que fundamentalmente se centra el presente trabajo. Así, durante más de diez a<math display="inline">\tilde{n}</math> os se utilizarón mallas convencionales en FD, pudiéndose destacar los trabajos de Altermann y Karal [[#bib0005|[1]]]  y Kelly et al. [[#bib0010|[2]]] .      

Los esquemas Staggered-grid (S-G) en FD fueron introducidos por Madariaga [[#bib0015|[3]]] , popularizándose su uso en problemas de propagación de ondas con los trabajos de de Virieux [[#bib0020|[4]]] , [[#bib0025|[5]]]  en los que se usa un esquema de segundo orden para la formulación en velocidad-tensión con ondas SH y P-SV respectivamente.      

Otros esquemas han sido presentados por Levander [[#bib0030|[6]]]  para una formulación en cuarto orden para ondas P-SV, Graves [[#bib0035|[7]]] , Graves et al. [[#bib0040|[8]]] , Ohnimato and Chovet [[#bib0045|[9]]]  que lo aplican en un esquema de segundo orden en la formulación desplazamiento-tensión. También pueden destacarse los trabajos de Kristek et al. [[#bib0050|[10]]] , [[#bib0055|[11]]] .      

A pesar de los trabajos citados de Virieux y Levander sobre dispersión y Crase et al. [[#bib0060|[12]]]  que obtienen condiciones de estabilidad para ondas P-SV para un orden arbitrario de aproximación, Moczo, Kristek y Halada [[#bib0065|[13]]]  se dan cuenta de que la dispersión de los esquemas S-G en 3-D no están suficientemente investigados y lo estudian para esquemas de cuarto orden en el espacio y segundo orden en el tiempo con validez para formulaciones en desplazamiento-tensión, velocidad-tensión y desplazamiento-velocidad-tensión. También Graves et al. [[#bib0040|[8]]]  han estudiado la estabilidad y Mozco et al. [[#bib0070|[14]]]  presentan trabajos recopilatorios de gran interés.      

Las mallas parcialmente Staggered(P-S-G) fueron utilizadas por primera vez por Andrews [[#bib0075|[15]]]  para medios anisótropos y Zhang [[#bib0080|[16]]]  las utiliza para la formulación velocidad-tensión en 2-D. Cruz-Atienza y Virieux [[#bib0085|[17]]]  las utilizan para la simulación de una rotura de falla en 2-D.      

El primero en utilizar una malla con espaciado variable fue Boore [[#bib0090|[18]]]  para un caso veritcal 1-D, Mikuano y Miyatarke [[#bib0095|[19]]]  usan un espaciado variable en un medio elástico en 2-D para analizar el proceso de rotura de falla. Pitarka [[#bib0100|[20]]]  presenta un esquema de cuarto orden para la formulación en velocidad-tensión sobre una malla rectangular con espaciado variable y Wang et al. [[#bib0105|[21]]]  combinan mallas con diferentes espaciados.      

Los métodos que en general se pueden denominar métodos sin malla, surgieron hace aproximadamente treinta a<math display="inline">\tilde{n}</math> os, aunque el interés y el esfuerzo investigador ha sido mínimo hasta hace poco tiempo. Así pues, la referencia más alejada en el tiempo y que puede condiderarse como punto de partida es el método de partículas (SPHM) que fue desarrollado por en 1977 por Lucy [[#bib0110|[22]]] . En 1982 y 1988 Monaghan ([[#bib0115|[23]]]  y [[#bib0120|[24]]] ) desarrolla una explicación más rigurosa del método.      

Una camino paralelo al indicado en el desarrollo de aproximaciones sin malla se basa en la idea de realizar una aproximación local mediante el método de mínimos cuadrados móviles de Lancaster y Salkauskas (1981) ([[#bib0125|[25]]] ). Nayroles, Touzot y Villon (1992) ([[#bib0130|[26]]] ) fueron los primeros en utilizar una aproximación de este tipo en un método de Galerkin, denominado método de los elementos difusos (DEM).      

Belytschko, Lu y Gu (1994) ([[#bib0135|[27]]] ) refinaron y desarrollaron una implementación alternativa del método de los elementos difusos, que mejoró la precisión al utilizar multiplicadores de Lagrange para imponer las condiciones de contorno esenciales y un orden de cuadratura mayor en la integración. Ellos denominaron al método Galerkin libre de elementos (EFG).      

Un salto importante en el desarrollo de estos métodos se ha producido a partir de los trabajos de Duarte y Oden (1995) ([[#bib0140|[28]]] ) y de Babuska y Melenk (1995) ([[#bib0145|[29]]] ) que realizan una generalización muy interesante de la aproximación por mínimos cuadrados móviles usando el concepto de partición de la unidad, que los primeros denominaron nubes hp (hp clouds) y de partición de la unidad (PUFEM) los segundos.      

Otro camino en la evolución de los métodos sin malla se ha constituido mediante el desarrollo del método de diferencias finitas generalizadas (GFDM) que surge como consecuencia de la evolución del método de diferencias finitas clásicas. Los trabajos de Jensen (1972) ([[#bib0150|[30]]] ) y Perrone y Kao (1975) ([[#bib0155|[31]]] ) pusieron las bases del método. Liszka y Orkisz (1980) ([[#bib0160|[32]]] ) y Orkisz (1998) ([[#bib0165|[33]]] ) han realizado importantes contribuciones a la mejora y desarrollo del método.      

El método de diferencias finitas generalizadas (GFDM) se basa en la utilización de mínimos cuadrados móviles, lo que permite aplicar esquemas en diferencias finitas en dominios irregulares. Al desarrollo de este método han contribuido los trabajos de los autores. Así, se han obtenido fórmulas explícitas en diferencias para el caso de estrellas irregulares, analizandose la influencia de los parámetros fundamentales que aperecen en la formulación [[#bib0170|[34]]] , se ha obtenido un método adaptativo para la mejora de la solución aproximada utilizando GFDM [[#bib0175|[35]]] , en [[#bib0180|[36]]]  se muestra la aplicación del método adaptativo en GFDM a la resolución de ecuaciones diferenciales en derivadas parciales de segundo orden y [[#bib0185|[37]]]  es un resumen global de GFDM que incluye las aportaciones realizadas por los autores y su aplicación a problemas elásticos.      

Los artículos [[#bib0190|[38]]]  y [[#bib0195|[39]]]  muestran la aplicación del método de diferencias finitas generalizadas a la resolución de ecuaciones en derivadas parciales dependientes del tiempo.      

Los métodos numéricos constituyen una herramienta irremplazable en la investigacíón de los movimientos sísmicos y entre ellos el actualmente dominante es el DF, es más Moczo et al. [[#bib0200|[40]]]  opinan que está por llegar el mejor momento de la aplicación del método Df en sismología. Tras estas consideraciones parece oportuno presentar las aplicación del GFDM a este tipo de problemas.      

Con el presente artículo se pretende establecer una base sólida para la aplicación de un método que ofrece grandes posibilidades entre las que inicialmente hay que destacar la capacidad de utilizar mallas irregulares directamente, evitando artificios en la realización de algunos modelos híbridos. Se desea hacer hincapié en la solidez metodológica puesto que no solo se presenta la formulación en GFD para la aplicación indicada, sino que se estudia de forma rigurosa su estabilidad y se analiza detalladamente la dispersión.

==2. Esquema en diferencias finitas generalizadas para el problema de propagación de ondas sísmicas en dominios homogéneos==

===2.1. Planteamiento del problema===

El problema de propagación de ondas sísmicas en 2-D viene caracterizado por el sistema de ecuaciones diferenciales en derivadas parciales

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}
\frac{{\partial }^2U(x,y,t)}{\partial t^2}={\alpha }^2\frac{{\partial }^2U(x,y,t)}{\partial x^2}+{\beta }^2\frac{{\partial }^2U(x,y,t)}{\partial y^2}+({\alpha }^2-{\beta }^2)\frac{{\partial }^2V(x,y,t)}{\partial x\partial y}\\
\frac{{\partial }^2V(x,y,t)}{\partial t^2}={\beta }^2\frac{{\partial }^2V(x,y,t)}{\partial x^2}+{\alpha }^2\frac{{\partial }^2V(x,y,t)}{\partial y^2}+({\alpha }^2-{\beta }^2)\frac{{\partial }^2U(x,y,t)}{\partial x\partial y}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

con las condiciones iniciales

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>U(x,y,0)=f_1(x,y);V(x,y,0)=f_2(x,y)</math>
|-
|<math>\frac{\partial U(x,y,0)}{\partial t}=f_3(x,y);\frac{\partial V(x,y,0)}{\partial t}=</math><math>f_4(x,y)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

y la condición de contorno

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}
a_1U(x_0,y_0,t)+b_1\frac{\partial U(x_0,y_0,t)}{\partial n}=g_1(t)\\
a_2V(x_0,y_0,t)+b_2\frac{\partial V(x_0,y_0,t)}{\partial n}=g_2(t)
\end{array}\quad en\quad \Gamma </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

siendo ''f''<sub>1</sub> (''x'' , ''y'' ), ''f''<sub>2</sub> (''x'' , ''y'' ), ''f''<sub>3</sub> (''x'' , ''y'' ), ''f''<sub>4</sub> (''x'' , ''y'' ), ''g''<sub>1</sub> (''t'' ) y ''g''<sub>2</sub> (''t'' ) funciones conocidas,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\alpha =\sqrt{\frac{\lambda +2\mu }{\rho }},\quad \beta =</math><math>\sqrt{\frac{\mu }{\rho }}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

''ρ''  es la densidad, ''λ''  cortante de Lamé y ''μ''  el módulo de rigidez. Las funciones ''U'' (''x'' , ''y'' , ''t'' ), ''V'' (''x'' , ''y'' , ''t'' ) son al menos dos veces diferenciables con continuidad en el dominio Ω con frontera Γ.      

===2.2. Esquema en diferencias finitas generalizadas===

Para la obtención de las fórmulas explícitas en diferencias finitas de las derivadas espaciales, una vez discretizado el dominio Ω ∪ Γ, a cada nodo del mismo se le asigna un conjunto de nodos a su alrededor que se denomina estrella. En la [[#fig0005|figura 1]]  se muestra una estrella irregular con ocho nodos seleccionados con el criterio de los cuatro cuadrantes (dos nodos por cuadrante, que son los más próximos al nodo central, considerando este como origen).

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_402164503-1-s2.0-S0213131511000150-gr1.jpg|center|241px|Estrella irregular de 9 nodos.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Estrella irregular de 9 nodos.

</span>
|}

La derivada segunda respecto del tiempo en el nodo central de la estrella se aproxima por:

<span id='eq0025'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}
\frac{{\partial }^2U(x_0,y_0,n\vartriangle t)}{\partial t^2}=\frac{u_0^{n+1}-2u_0^n+u_0^{n-1}}{(\vartriangle t)^2}\\
\frac{{\partial }^2V(x_0,y_0,n\vartriangle t)}{\partial t^2}=\frac{v_0^{n+1}-2v_0^n+v_0^{n-1}}{(\vartriangle t)^2}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

donde <math display="inline">u_0^n</math> , <math display="inline">v_0^n</math> , <math display="inline">u_0^{n+1}</math>  y <math display="inline">v_0^{n+1}</math>  son los valores aproximados de las funciones ''U'' (''x'' , ''y'' , ''t'' ), ''V'' (''x'' , ''y'' , ''t'' ) en el nodo central de coordenadas espaciales (''x''<sub>0</sub> , ''y''<sub>0</sub> ) para los tiempos ''n''  ▵ ''t''  y (''n''  + 1) ▵ ''t''  respectivamente.      

Las derivadas espaciales se aproximan por las expresiones explícitas en diferencias finitas generalizadas, obtenidas en [[#bib0180|[36]]]  y [[#bib0190|[38]]] :

<span id='eq0030'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\partial }^2U(x_0,y_0,n\vartriangle t)}{\partial x^2}=</math><math>-m_0u_0^n+\sum_{j=1}^Nm_ju_j^n;\quad \frac{{\partial }^2V(x_0,y_0,n\vartriangle t)}{\partial x^2}=</math><math>-m_0v_0^n+\sum_{j=1}^Nm_jv_j^n</math>
|-
|<math>\frac{{\partial }^2U(x_0,y_0,n\vartriangle t)}{\partial y^2}=</math><math>-{\eta }_0u_0^n+\sum_{j=1}^N{\eta }_ju_j^n;\quad \frac{{\partial }^2V(x_0,y_0,n\vartriangle t)}{\partial y^2}=</math><math>-{\eta }_0v_0^n+\sum_{j=1}^N{\eta }_jv_j^n</math>
|-
|<math>\frac{{\partial }^2U(x_0,y_0,n\vartriangle t)}{\partial x\partial y}=</math><math>-{\zeta }_0u_0^n+\sum_{j=1}^N{\zeta }_ju_j^n;\quad \frac{{\partial }^2V(x_0,y_0,n\vartriangle t)}{\partial x\partial y}=</math><math>-{\zeta }_0v_0^n+\sum_{j=1}^N{\zeta }_jv_j^n</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

donde ''N''  es el número de nodos de la estrella cuyo nodo central tiene de coordenadas (''x''<sub>0</sub> , ''y''<sub>0</sub> ) (en este trabajo ''N''  = 8 y son seleccionados mediante el criterio de los cuatro cuadrantes), ''m''<sub>0</sub> , ''η''<sub>0</sub>  y ''ζ''<sub>0</sub>  son los coeficientes que multiplican a los valores aproximados de las funciones ''U''  y ''V''  en el nodo central para el tiempo ''n''  ▵ ''t  ''  (<math display="inline">u_0^n</math>  y <math display="inline">v_0^n</math>  respectivamente) de las expresiones explícitas en diferencias finitas generalizadas para las derivadas espaciales <math display="inline">\frac{{\partial }^2}{\partial x^2}</math> , <math display="inline">\frac{{\partial }^2}{\partial y^2}</math>  y <math display="inline">\frac{{\partial }^2}{\partial x\partial y}</math>  respectivamente, ''m''<sub>''j''</sub> , ''η''<sub>''j''</sub>  y ''ζ''<sub>''j''</sub>  son los coeficientes que multiplican a los valores aproximados de las funciones ''U''  y ''V''  en el resto de los nodos de la estrella para el tiempo ''n''  ▵ ''t  ''  (<math display="inline">u_j^n</math>  y <math display="inline">v_j^n</math>  respectivamente) de las expresiones explícitas en diferencias finitas generalizadas para las derivadas espaciales <math display="inline">\frac{{\partial }^2}{\partial x^2}</math> , <math display="inline">\frac{{\partial }^2}{\partial y^2}</math>  y <math display="inline">\frac{{\partial }^2}{\partial x\partial y}</math>  respectivamente ([[#bib0160|[32]]] , [[#bib0165|[33]]] , [[#bib0170|[34]]] , [[#bib0185|[37]]] ).      

Si en el sistema de ecuaciones [[#eq0005|1]]  se sustituyen las derivadas parciales por las expresiones [[#eq0025|4]]  y [[#eq0030|5]]  se obtiene:

<span id='eq0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}
u_0^{n+1}=2u_0^n-u_0^{n-1}+(\vartriangle t)^2[{\alpha }^2(-m_0u_0^n+\sum_1^Nm_ju_j^n)+{\beta }^2(-{\eta }_0u_0^n+\sum_1^N{\eta }_ju_j^n)\\
+({\alpha }^2-{\beta }^2)(-{\zeta }_0v_0^n+\sum_1^N{\zeta }_jv_j^n)]\\
v_0^{n+1}=2v_0^n-v_0^{n-1}+(\vartriangle t)^2[{\beta }^2(-m_0v_0^n+\sum_1^Nm_jv_j^n)+{\alpha }^2(-{\eta }_0v_0^n+\sum_1^N{\eta }_jv_j^n)\\
+({\alpha }^2-{\beta }^2)(-{\zeta }_0u_0^n+\sum_1^N{\zeta }_ju_j^n)]
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

que constituyen el esquema en diferencias finitas generalizadas para el pro-blema planteado.      

==3. Criterio de estabilidad==

Para estudiar la estabilidad se utiliza el análisis de von Neumann. Para lo cual se realiza una descomposición armónica de la solución aproximada de la forma

<span id='eq0040'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u_0^n=A{\xi }^ne^{i\boldsymbol{k}^T\boldsymbol{x}_\boldsymbol{0}};\quad u_j^n=</math><math>A{\xi }^ne^{i\boldsymbol{k}^T\boldsymbol{x}_\boldsymbol{j}};\quad v_0^n=</math><math>B{\xi }^ne^{i\boldsymbol{k}^T\boldsymbol{x}_\boldsymbol{0}};\quad v_j^n=</math><math>B{\xi }^ne^{i\boldsymbol{k}^T\boldsymbol{x}_\boldsymbol{j}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

donde '''''x'''''<sub>'''0'''</sub>  es vector de posición del nodo central de la estrella, '''''x'''''<sub>'''j'''</sub> , ''j''  = 1, …, ''N''  son los vectores de posición del resto de los nodos de la estrella y '''h'''<sub>'''j'''</sub>  son los vectores de las posiciones relativas de los nodos de la estrella con respecto al nodo central cuyas coordenadas son ''h''<sub>''jx''</sub>  = ''x''<sub>''j''</sub>  − ''x''<sub>0</sub> ,  ''h''<sub>''jy''</sub>  = ''y''<sub>''j''</sub>  − ''y''<sub>0</sub> .      

''ξ''  = ''e''<sup>−''iw'' ▵''t''</sup>  es el llamado factor de amplificación, cuyo valor determinará la condición de estabilidad, <math display="inline">w</math>  es la frecuencia angular y '''''k'''''  es el vector del número de onda (ver [[#fig0010|figura 2]] )

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_402164503-1-s2.0-S0213131511000150-gr2.jpg|center|226px|Número de onda k→.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Número de onda <math display="inline">\boldsymbol{\overset{\rightarrow}{\mbox{k}}}</math> .                  

</span>
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{k}=\left\{\begin{array}{c}
{\kappa }_x\\
{\kappa }_y
\end{array}\right\}=\kappa \left\{\begin{array}{c}
cos\varphi \\
sin\varphi 
\end{array}\right\}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Al sustituir las expresiones [[#eq0040|7]]  en la ecuación [[#eq0035|6]] , y simplificar, se obtiene.

<span id='eq0050'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>A\xi =2A-\frac{A}{\xi }+(\vartriangle t)^2[{\alpha }^2(-</math><math>{Am}_0+A\sum_1^Nm_je^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>{\beta }^2(-A{\eta }_0+A\sum_1^N{\eta }_je^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>({\alpha }^2-{\beta }^2)(-B{\zeta }_0+B\sum_1^N{\zeta }_je^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})]</math>
|-
|<math>B\xi =2B-\frac{B}{\xi }+(\vartriangle t)^2[{\beta }^2(-</math><math>{Bm}_0+B\sum_1^Nm_je^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>{\alpha }^2(-B{\eta }_0+B\sum_1^N{\eta }_je^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>({\alpha }^2-{\beta }^2)(-A{\zeta }_0+A\sum_1^N{\zeta }_je^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

Teniendo en cuenta que los coeficientes ''m''<sub>''j''</sub> , ''η''<sub>''j''</sub> , ''ζ''<sub>''j''</sub> , son particiones de la unidad, se tiene

<span id='eq0055'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>m_0=\sum_1^Nm_j;\quad {\eta }_0=\sum_1^N{\eta }_j;\quad {\zeta }_0=</math><math>\sum_1^N{\zeta }_j</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

al sustituir las expresiones [[#eq0055|9]]  en la ecuación [[#eq0050|8]] , se obtiene el sistema de dos ecuaciones con dos incógnitas, ''A''  y ''B'' , siguiente

<span id='eq0060'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>A[\xi -2+\frac{1}{\xi }+(\vartriangle t)^2{\alpha }^2\sum_1^Nm_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>(\vartriangle t)^2{\beta }^2\sum_1^N{\eta }_j(1-e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})]</math>
|-
|<math>+B(\vartriangle t)^2({\alpha }^2-{\beta }^2)\sum_1^N{\zeta }_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})=</math><math>0</math>
|-
|<math>A(\vartriangle t)^2({\alpha }^2-{\beta }^2)\sum_1^N{\zeta }_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>B[\xi -2+\frac{1}{\xi }+(\vartriangle t)^2{\beta }^2\sum_1^Nm_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})</math>
|-
|<math>+(\vartriangle t)^2{\alpha }^2\sum_1^N{\eta }_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})]=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

en el que al despejar ''B''  de la segunda ecuación del sistema  [[#eq0060|10]] , y sustituir en la primera ecuación, y tras operar, se obtiene

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>[2cosw\vartriangle t-2+(\vartriangle t)^2({\alpha }^2\sum_1^Nm_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>{\beta }^2\sum_1^N{\eta }_j(1-e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}}))]\times </math>
|-
|<math>[2cosw\vartriangle t-2+(\vartriangle t)^2({\beta }^2\sum_1^Nm_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})+</math><math>{\alpha }^2\sum_1^N{\eta }_j(1-e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}}))]</math>
|-
|<math>=(\vartriangle t)^4({\alpha }^2-{\beta }^2)^2[\sum_1^N{\zeta }_j(1-</math><math>e^{i\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}})]^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

Tras operar, se obtienen las siguientes condiciones:

==Parte real==

<span id='eq0070'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>(1-cosw\vartriangle t)^2-2(1-cosw\vartriangle t)\frac{(\vartriangle t)^2}{4}({\alpha }^2+</math><math>{\beta }^2)\sum_1^N(m_j+{\eta }_j)(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})+</math><math></math>
|-
|<math>\frac{(\vartriangle t)^4}{4}[({\alpha }^2\sum_1^Nm_j(1-</math><math>cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})+</math><math>{\beta }^2\sum_1^N{\eta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}))({\beta }^2\sum_1^Nm_j(1-</math><math>cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})</math>
|-
|<math>+{\alpha }^2\sum_1^N{\eta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}))-</math><math>({\alpha }^2\sum_1^Nm_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}+</math><math>{\beta }^2\sum_1^N{\eta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})({\beta }^2\sum_1^Nm_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}</math>
|-
|<math>+{\alpha }^2\sum_1^N{\eta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})-</math><math>({\alpha }^2-{\beta }^2)^2[(\sum_1^N{\zeta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}))^2-</math><math>(\sum_1^N{\zeta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})^2]]=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)
|}

==Parte imaginaria==

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>2(1-cosw\vartriangle t)({\alpha }^2+{\beta }^2)\sum_1^N(m_j+</math><math>{\eta }_j)sin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}-</math><math>(\vartriangle t)^2[({\alpha }^2\sum_1^Nm_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}</math>
|-
|<math>+{\beta }^2\sum_1^N{\eta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})({\beta }^2\sum_1^Nm_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}+</math><math>{\alpha }^2\sum_1^N{\eta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})+</math><math></math>
|-
|<math>({\alpha }^2\sum_1^Nm_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}+</math><math>{\beta }^2\sum_1^N{\eta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})({\beta }^2\sum_1^Nm_j(1-</math><math>cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})+</math><math></math>
|-
|<math>{\alpha }^2\sum_1^N{\eta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}))+</math><math>2({\alpha }^2-{\beta }^2)^2(\sum_1^N{\zeta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})\sum_1^N{\zeta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}]=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

Al operar con estas ecuaciones, simplificando mediante el uso de criterios conservadores y eligiendo entre ambas la condición más restrictiva, se obtiene

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\vartriangle t<\sqrt{\frac{4}{\left({\alpha }^2+{\beta }^2)[(\vert m_0\vert +\vert {\eta }_0\vert )+\sqrt{(m_0+{\eta }_0)^2+{\zeta }_0^2}\right]}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

que es la condición de estabilidad de la estrella. Por tanto, el paso de tiempo (▵''t'' ) deberá ser el mínimo de los obtenidos para cada una de las estrellas de la malla.      

==4. Dispersión en la estrella==

===4.1. Dispersión de las ondas P===

De la ecuación [[#eq0070|12]]  se tiene

<span id='eq0085'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\omega =\frac{1}{\vartriangle t}arccos\Phi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Phi =1-\frac{(\vartriangle t)^2}{4}(({\alpha }^2+</math><math>{\beta }^2)(a_1+a_3)+(({\alpha }^2+{\beta }^2)^2(a_1+</math><math>a_3)^2+</math>
|-
|<math>4[({\alpha }^2-{\beta }^2)^2(a_5^2-a_6^2)+({\alpha }^2a_2+</math><math>{\beta }^2a_4)({\beta }^2a_2+{\alpha }^2a_4)-({\alpha }^2a_1+</math><math>{\beta }^2a_3)({\beta }^2a_1+{\alpha }^2a_3)])^{\frac{1}{2}})</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a_1=\sum_1^Nm_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})\Rightarrow \frac{\partial a_1}{\partial k}=</math><math>a_{1,k}=\sum_1^Nm_jdsinkd</math>
|-
|<math>a_2=\sum_1^Nm_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}\Rightarrow \frac{\partial a_2}{\partial k}=</math><math>a_{2,k}=\sum_1^Nm_jdcoskd</math>
|-
|<math>a_3=\sum_1^N{\eta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})\Rightarrow \frac{\partial a_3}{\partial k}=</math><math>a_{3,k}=\sum_1^N{\eta }_jdsinkd</math>
|-
|<math>a_4=\sum_1^N{\eta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}\Rightarrow \frac{\partial a_4}{\partial k}=</math><math>a_{4,k}=\sum_1^N{\eta }_jdcoskd</math>
|-
|<math>a_5=\sum_1^N{\zeta }_j(1-cos\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j})\Rightarrow \frac{\partial a_5}{\partial k}=</math><math>a_{5,k}=\sum_1^N{\zeta }_jdsinkd</math>
|-
|<math>a_6=\sum_1^N{\zeta }_jsin\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}\Rightarrow \frac{\partial a_6}{\partial k}=</math><math>a_{6,k}=\sum_1^N{\zeta }_jdcoskd</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 17)
|}

y

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{k}^T\boldsymbol{h}_\boldsymbol{j}=k(h_{jx}cos\varphi +</math><math>h_{jy}sin\varphi )=kd</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Se sabe que

<span id='eq0105'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\omega =2\pi \frac{c^{grid}}{{\lambda }^{grid}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

donde ''c''<sup>''grid''</sup>  y ''λ''<sup>''grid''</sup>  son la velocidad de fase (''α''<sup>''grid''</sup>  o ''β''<sup>''grid''</sup> ) y la longitud de onda (<math display="inline">{\lambda }_P^{grid}</math>  y <math display="inline">{\lambda }_S^{grid}</math> ) en la estrella respectivamente ([[#bib0020|[4]]] , [[#bib0025|[5]]] , [[#bib0030|[6]]] , [[#bib0035|[7]]] , [[#bib0060|[12]]] , [[#bib0065|[13]]] , [[#bib0070|[14]]] ).      

Si se definen las relaciones:

<span id='eq0110'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s=\frac{2}{{\lambda }_S^{grid}\sqrt{\left(r^2+1)[(\vert m_0\vert +\vert {\eta }_0\vert )+\sqrt{(m_0+{\eta }_0)^2+{\zeta }_0^2}\right]}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

<span id='eq0115'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s_P=\frac{2}{{\lambda }_P^{grid}\sqrt{\left(r^2+1)[(\vert m_0\vert +\vert {\eta }_0\vert )+\sqrt{(m_0+{\eta }_0)^2+{\zeta }_0^2}\right]}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

<span id='eq0120'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>p=\frac{\beta \vartriangle t\sqrt{\left(r^2+1)[(\vert m_0\vert +\vert {\eta }_0\vert )+\sqrt{(m_0+{\eta }_0)^2+{\zeta }_0^2}\right]}}{2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

La razón de velocidades es

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>r=\frac{\alpha }{\beta }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)
|}

y por tanto

<span id='eq0130'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s_P=\frac{s}{r}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 23)
|}

Si se sutituyen las ecuaciones [[#eq0085|(15)]]  and [[#eq0115|(20)]] , [[#eq0120|21]]  y [[#eq0130|23]]  en la ecuación [[#eq0105|18]]  particularizada para las ondas P, tras operar, se obtiene la relación de dispersión para las ondas ''P''

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\alpha }^{grid}}{\alpha }=\frac{arccos\Phi }{2\pi sp}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)
|}

===4.2. Dispersión de las ondas S===

Se sabe que

<span id='eq0140'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\omega =2\pi \frac{{\beta }^{grid}}{{\lambda }_S^{grid}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)
|}

Si se sustituyen las ecuaciones [[#eq0085|(15)]]  and [[#eq0110|(19)]]  y [[#eq0120|21]]  en la ecuación [[#eq0140|25]] , tras operar, se obtiene

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\beta }^{grid}}{\beta }=\frac{arccos\Phi }{2\pi sp}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 26)
|}

===4.3. Dispersión de la velocidad de grupo ondas P===

Por definición la velocidad de grupo para las ondas ''P  '' , se obtiene al derivar <math display="inline">w</math> , dada por la expresión [[#eq0085|15]] , respecto de ''k'' , por tanto

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\alpha }_{group}^{grid}=\frac{\partial w}{\partial k}=</math><math>\frac{\vartriangle t}{4}\frac{{\beta }^2\Upsilon }{\sqrt{1-{\Phi }^2}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 27)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Upsilon =(r^2+1)(a_{1,k}+a_{3,k})+\frac{1}{2}[2(r^2+</math><math>1)^2(a_1+a_3)(a_{1,k}+a_{3,k})+</math>
|-
|<math>4[2(r^2-1)^2(a_5a_{5,k}-a_6a_{6,k})+(r^2a_{2,k}+a_{4,k})(a_2+</math><math>r^2a_4)+</math>
|-
|<math>(r^2a_2+a_4)(a_{2,k}+r^2a_{4,k})-(r^2a_{1,k}+a_{3,k})(a_1+</math><math>r^2a_3)-</math>
|-
|<math>(r^2a_1+a_3)(a_{1,k}+r^2a_{3,k})]]\times [(r^2+1)^2(a_1+</math><math>a_3)^2+</math>
|-
|<math>4[(r^2-1)^2(a_5^2-a_6^2)+(r^2a_2+a_4)(a_2+r^2a_4)-</math><math>(r^2a_1+a_3)(a_1+r^2a_3)]]^{-\frac{1}{2}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 28)
|}

designando

<span id='eq0160'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F=(r^2+1)(a_1+a_3)+[(r^2+1)^2(a_1+a_3)^2+</math>
|-
|<math>4[(r^2-1)^2(a_5^2-a_6^2)+(r^2a_2+a_4)(a_2+r^2a_4)-</math><math>(r^2a_1+a_3)(a_1+r^2a_3)]]^{\frac{1}{2}}]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)
|}

Teniendo en cuenta las ecuación [[#eq0160|29]]  y [[#eq0120|21]] , se tiene

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\alpha }_{group}^{grid}}{\alpha }=\frac{1}{2\sqrt{2}r}\frac{\Upsilon }{\sqrt{F-(\frac{pF}{\sqrt{\left(r^2+1)[(\vert m_0\vert +\vert {\eta }_0\vert )+\sqrt{(m_0+{\eta }_0)^2+{\zeta }_0^2}\right]}\sqrt{2}})^2}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 30)
|}

===4.4. Dispersión de la velocidad de grupo ondas S===

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\beta }_{group}^{grid}}{\beta }=\frac{1}{2\sqrt{2}}\frac{\Upsilon }{\sqrt{F-(\frac{pF}{\sqrt{\left(r^2+1)[(\vert m_0\vert +\vert {\eta }_0\vert )+\sqrt{(m_0+{\eta }_0)^2+{\zeta }_0^2}\right]}\sqrt{2}})^2}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)
|}

===4.5. Irregularidad de la estrella y dispersión===

Los coeficientes, ''m''<sub>0</sub> , ''η''<sub>0</sub> , ''ζ''<sub>0</sub> , que aparecen en el criterio de estabilidad dependen del número de nodos en la estrella, de las coordenadas de cada nodo respecto del nodo central de la estrella y de la función de ponderación elegida ([[#bib0170|[34]]] ).      

Fijados el número de nodos por estrella, en nuestro caso 9 (''N''  = 8), y la función de ponderación, <math display="inline">\Omega (h_{jx},h_{jy})=\frac{1}{(\sqrt{h_{jx}^2+h_{jy}^2})^3}</math> , la expresión

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{1}{\sqrt{\left(r^2+1)[(\vert m_0\vert +\vert {\eta }_0\vert )+\sqrt{(m_0+{\eta }_0)^2+{\zeta }_0^2}\right]}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 32)
|}

depende de la razón de velocidades y de las coordenadas relativas de cada nodo de la estrella respecto del nodo central de la misma.      

Al disminuir las coordenadas relativas de los nodos de la estrella, aumentan los coeficientes ''m''<sub>0</sub> , ''η''<sub>0</sub> , ''ζ''<sub>0</sub> , lo cual ocurre al aumentar el número de nodos en la malla. Para un mismo número de nodos en la malla el mínimo valor de los coeficientes mencionados se obtiene para la malla regular ([[#bib0180|[36]]] ).      

Sea ''τ''<sub>''l''</sub>  la media de las distancias de los nodos de la estrella ''l''  a su nodo central. Sea ''τ''  la media de todas las ''τ''<sub>''l''</sub>  de las estrellas de la malla, una vez discretizado nuestro dominio. Se tiene que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{h}_\boldsymbol{j}=\tau \left\{\begin{array}{c}
\overline{h_{jx}}\\
\overline{h_{jy}}
\end{array}\right\}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overline{m_0}=m_0{\tau }^2;\quad \overline{{\eta }_0}=</math><math>{\eta }_0{\tau }^2;\quad \overline{{\eta }_0}={\zeta }_0{\tau }^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 34)
|}

en las que la barra superior indica valores medios.      

La condición de estabilidad puede reescribirse

<span id='eq0190'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\vartriangle t<\frac{2\tau }{\beta \sqrt{\left(r^2+1\right)}\sqrt{(\vert \overline{m_0}\vert +\vert \overline{{\eta }_0}\vert )+\sqrt{(\overline{m_0}+\overline{{\eta }_0})^2+{\overline{{\zeta }_0}}^2}}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)
|}

En el caso de malla regular de ocho nodos la expresión [[#eq0190|35]]  resulta ser

<span id='eq0195'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\vartriangle t<\frac{\tau }{\beta \sqrt{r^2+1}}\frac{2(\sqrt{2}-1)\sqrt{3}}{\sqrt{5}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)
|}

Si se multiplica el segundo miembro de la expresión [[#eq0195|36]]  por el factor

<span id='eq0200'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{IIS}_0=\frac{\sqrt{5}(\sqrt{2}+1)}{\sqrt{3(\vert \overline{m_0}\vert +\vert \overline{{\eta }_0}\vert +\sqrt{(\overline{m_0}+\overline{{\eta }_0})^2+{\overline{{\zeta }_0}}^2})}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)
|}

se obtiene la expresión [[#eq0190|35]] .      

El factor obtenido en [[#eq0200|37]] , que para una estrella regular (esquema de ocho nodos más el nodo central) es la unidad, resulta ser un indicador de la irregularidad de la estrella, ''IIS''<sub>0</sub> , cuyo nodo central tiene de coordenadas (''x''<sub>0</sub> , ''y''<sub>0</sub> ). Se adopta como indicador de irregularidad de la malla (IIC) el valor mínimo de los indicadores de irregularidad obtenidos para cada una de las estrellas de la malla.      

El control sobre la irregularidad de la malla es importante, ya que al disminuir el ''IIC  '' , disminuye el paso del tiempo. Sin embargo, mediante un proceso adaptativo de la malla (a<math display="inline">\tilde{n}</math> adiendo nodos de forma selectiva [[#bib0175|[35]]] ) se pueden conseguir índices de irregularidad más próximos a la unidad y, por tanto, pasos de tiempo mayores y menores relaciones de dispersión.      

==5. Resultados numéricos==

===5.1. Irregularidad y estabilidad===

En esta subsection se muestran los resultados obtenidos al resolver el problema definido en 2.1 en Ω = [0, 1] × [0, 1] ⊂ '''R'''<sup>2</sup> , con las condiciones de contorno Dirichlet y las condiciones iniciales

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>U(x,y,0)=sinxsiny;V(x,y,0)=cosxcosy</math>
|-
|<math>\frac{\partial U(x,y,0)}{\partial t}=0;\frac{\partial V(x,y,0)}{\partial t}=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)
|}

utilizando irregulares (ver figuras [[#fig0020|4]]  y [[#fig0025|5]] ) de 121 nodos. La solución analítica

<span id='eq0210'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>U(x,y,t)=cos(\sqrt{2}\beta t)sinxsiny;\quad V(x,y,t)=</math><math>cos(\sqrt{2}\beta t)cosxcosy</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)
|}

La función de ponderación utilizada ha sido

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Omega (h_{jx},h_{jy})=\frac{1}{\sqrt{(h_{jx}^2+h_{jy}^2)^3}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)
|}

y el criterio de selección de los nodos es el del cuadrante. El error global ha sido calculado para cada paso de tiempo utilizando la siguiente norma

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Error\quad global=\frac{\sqrt{\frac{{\sum }_{j=1}^{NT}(sol(j)-exac(j))^2}{NT}}}{\vert {exac}_{max}\vert }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 41)
|}

donde ''sol'' (''j'' ) es el valor de la solución aproximada en el nodo ''j'' , ''exac'' (''j'' ) es el valor de la solución exacta en el nodo ''j'' , ''exac''<sub>''max''</sub>  es el máximo valor de la solución exacta en los nodos interiores de de la malla considerada y ''NT''  es el número de nodos del interior.

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_402164503-1-s2.0-S0213131511000150-gr4.jpg|center|469px|Malla irregular (IIC=0, 65).]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Malla irregular (''IIC''  = 0, 65).                  

</span>
|}

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_402164503-1-s2.0-S0213131511000150-gr5.jpg|center|494px|Malla irregular (IIC=0, 89).]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

Malla irregular (''IIC''  = 0, 89).                  

</span>
|}

Los valores de los errores globales, para ''n''  = 500, al variar ▵''t'' , utilizando la malla irregular de 121 nodos de la  [[#fig0020|figura 4]] , con ''IIC''  = 0, 6524, se muestran en la [[#tbl0005|tabla 1]] .

<span id='tbl0005'></span>

Tabla 1.

Influencia del ▵''t''  en el error global con ''α''  = 1 ; ''β''  = 0, 5 ; ''IIC''  = 0, 6524 ; ''n''  = 500  

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"

|-

! Error Global ''V''
! ▵''t''
! Error Global ''U''
|-

| 0,0316
| 2, 078 × 10<sup>−3</sup>
| 9, 875 × 10<sup>−4</sup>
|-

| 0,0223
| 7, 423 × 10<sup>−4</sup>
| 3, 123 × 10<sup>−4</sup>
|-

| 0,01
| 2, 514 × 10<sup>−4</sup>
| 1, 217 × 10<sup>−4</sup>
|-

| 0,007
| 1, 423 × 10<sup>−4</sup>
| 1, 012 × 10<sup>−4</sup>
|}

Los valores de los errores globales, para ''n''  = 500, al variar ▵''t'' , utilizando la malla irregular de 121 nodos de la  [[#fig0025|figura 5]] , con ''IIC''  = 0, 8944, se muestran en la [[#tbl0010|tabla 2]] .

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 2.

Influencia del ▵''t''  en el error global con ''α''  = 1 ; ''β''  = 0, 5 ; ''IIC''  = 0, 8944 ; ''n''  = 500                  

|-

! ▵''t''
! Error Global ''U''
! Error Global ''V''
|-

| 0,0316
| 7, 331 × 10<sup>−4</sup>
| 5, 167 × 10<sup>−4</sup>
|-

| 0,0223
| 2, 685 × 10<sup>−4</sup>
| 1, 950 × 10<sup>−4</sup>
|-

| 0,01
| 1, 102 × 10<sup>−4</sup>
| 8, 428 × 10<sup>−5</sup>
|-

| 0,007
| 9, 575 × 10<sup>−5</sup>
| 7, 340 × 10<sup>−5</sup>
|}

===5.2. Resultados para distintas longitudes de onda===

La [[#tbl0015|tabla 3]]  muestra los resultados de los errores globales de ''U''  y ''V'' , para ''n''  = 500 y ▵''t''  = 0, 01, al resolver el problema definido en 2.1, utilizando una malla regular de 121 nodos ([[#fig0015|figura 3]] ), para funciones cuyas soluciones analíticas son: [[#eq0210|39]]  y

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 3.

Errores globales para distintas soluciones analíticas del problema 2.1, ''n''  = 50 ; ▵ ''t''  = 0, 01                  

|-

! Ecuación (sol. anal.)
! Error global U
! Error global V
|-

| 42
| 1, 646 × 10<sup>−6</sup>
| 1, 778 × 10<sup>−6</sup>
|-

| 39
| 1, 232 × 10<sup>−5</sup>
| 2, 443 × 10<sup>−5</sup>
|-

| 43
| 4, 081 × 10<sup>−4</sup>
| 2, 001 × 10<sup>−4</sup>
|-

| 44
| 9, 188 × 10<sup>−2</sup>
| 8, 647 × 10<sup>−2</sup>
|-

| 45
| 3, 035 × 10<sup>−1</sup>
| 3, 275 × 10<sup>−1</sup>
|-

| 46
| 4, 942 × 10<sup>−1</sup>
| 4, 999 × 10<sup>−1</sup>
|}

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_402164503-1-s2.0-S0213131511000150-gr3.jpg|center|489px|Malla regular.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Malla regular.

</span>
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}
U(x,y,t)=cos(0,5\sqrt{2}\beta t)sin(0,5x)sin(0,5y)\\
V(x,y,t)=cos(0,5\sqrt{2}\beta t)cos(0,5x)cos(0,5y)
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 42)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}
U(x,y,t)=cos(2\sqrt{2}\beta t)sin(2x)sin(2y)\\
V(x,y,t)=cos(2\sqrt{2}\beta t)cos(2x)cos(2y)
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 43)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}
U(x,y,t)=cos(4\pi \beta t)sin(2\sqrt{2}\pi x)sin(2\sqrt{2}\pi y)\\
V(x,y,t)=cos(4\pi \beta t)cos(2\sqrt{2}\pi x)cos(2\sqrt{2}\pi y)
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 44)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}
U(x,y,t)=cos(8\pi \beta t)sin(4\sqrt{2}\pi x)sin(4\sqrt{2}\pi y)\\
V(x,y,t)=cos(8\pi \beta t)cos(4\sqrt{2}\pi x)cos(4\sqrt{2}\pi y)
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 45)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}
U(x,y,t)=cos(16\pi \beta t)sin(8\sqrt{2}\pi x)sin(8\sqrt{2}\pi y)\\
V(x,y,t)=cos(8\pi \beta t)cos(8\sqrt{2}\pi x)cos(8\sqrt{2}\pi y)
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 46)
|}

y las correspondientes condiciones de contorno Dirichlet e iniciales      

===5.3. Resultados de la dispersión===

Los valores obtenidos de la dispersión en cada una de las estrellas de la malla irregular con ''IIC''  = 0, 8944 de las ondas ''P'' , para tres valores distintos de del ángulo ''φ''  = 0, 45, 90, se han representado sobre los nodos en la [[#fig0030|figura 6]] . Puesto que los valores en los tres casos únicamente difieren a partir de la tercera cifra decimal la [[#fig0030|figura 6]]  solo muestra el resultado para el caso ''φ''  = 45

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_402164503-1-s2.0-S0213131511000150-gr6.jpg|center|362px|Dispersión de las ondas P (IIC=0, 89).]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Dispersión de las ondas P (''IIC''  = 0, 89).                  

</span>
|}

==6. Conclusiones==

En este artículo se muestran las expresiones explíctas, en diferencias finitas generalizadas, para la propogación de ondas sísmicas en 2-D. Se estudia la estabilidad, obteniéndose un criterio en función de la razón de velocidades y los coeficientes de las fórmulas explícitas.

Se analiza la dispersión y se relaciona con la irregularidad de la estrella utilizando el índice de irregularidad de la malla. La utilización de mallas irregulares, según la geometría del problema, puede suponer dispersiones elevadas en algunas estrellas, lo cual va asociado con valores peque<math display="inline">\tilde{n}</math> os del ''IIC''  y un menor paso de tiempo. En este caso se puede aumentar el ''IIC'' , y por tanto, el paso del tiempo y disminuir las dispersiones, realizando una redefinición de la malla, para lo que se puede utilizar un proceso adaptativo como el presentado por los autores en  [[#bib0175|[35]]] .      

Por último es fácil ver como el error aumenta al disminuir el número de puntos de la discretización por longitud de onda, apreciándose que hasta el orden de diez puntos por longitud de onda se obtienen errores generalmente admisibles inferiores al 10%.

==Agradecimientos==

Los autores agradecen la ayuda recibida del Ministerio de Ciencia e Innovación de Espa<math display="inline">\tilde{n}</math> a en el proyecto TISMANCA, Ref.: CGL2008-01757/CLI.

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] Z. Alterman, F.C. Karal; Propagation of elastic waves in layered media by finite-difference methods, Bulll; Seism. Society American, 58 (1968), pp. 367–398</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] L.R. Kelly, R.W. Ward, S. Freitel, R.M. Alford; Synthetic seismograms: a finite-difference approach; Geophysics, 41 (1976), pp. 2–27</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] R. Madariaga; Dynamics of an expanding circular fault, Bull; Seism. Soc. Am., 67 (1976), pp. 163–182</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] J. Virieux; SH-wave propagation in heterogeneus media: velocity-stress finite difference method; Geophysics, 49 (1984), pp. 1933–1957</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] J. Virieux; P-SV wave propagation in heterogeneus media: velocity-stress finite difference method; Geophysics, 51 (1986), pp. 889–901</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] A.R. Levander; Fourth-order finite-difference P-SV seismograms; Geophysics, 53 (1988), pp. 1425–1436</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] R.W. Graves; Simulating seimisc wave propagation in 3-D elastic media using staggered-grid finite differences, Bull; Seism. Soc. Am., 93 (1996), pp. 283–300</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] R.W. Graves, S.M. Day; Stability and accuracy analysis of course-grain viscoelastic simulations, Bull; Seism. Soc. Am., 87 (2003), pp. 494–515</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] T. Ohnimato, B.A. Chovet; A free-surface boundary condition for including 3-D topography in finite difference method, Bull; Seism. Soc. Am., 87 (1997), pp. 494–515</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] J. Kristek, P. Moczo, P.J. Archuleta; Efficient method to simulate planar free surface in the 3-D 4<sup>''th''</sup> order staggered-grid finite difference schemes, Studia Geophys                                        ; Geod., 46 (2002), pp. 355–381</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] J. Kristek, P. Moczo; Seismic wavw propagation in viscoelastic media with material discontinuities- a 3-D 4<sup>''th''</sup> order staggered grid finite-difference modeling, Bull                                        ; Seism. Soc. Am., 93 (2003), pp. 2273–2280</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] E. Crase, Ch. Wideman, N. Noble, A. Tarantula; Non linear elastic waveform inversion of land seismic reflection data; J. Geophys. Rev., 97 (1992), pp. 4685–4703</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] P. Moczo, J. Kristek, L. Halada; 3-D 4<sup>''th''</sup> order staggered-grid finite-difference schemes: stability and grid dispersion, Bull                                        ; Seism. Soc. Am., 90 (2000), pp. 587–603</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] P. Moczo, J. Kristek and L. Halada, The finite-difference method for siemologists. An introduction, Comenius University, Bratislava.</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] D.J. Andrews; A numerical study of tectonic stress relase by underground explosions, Bulll; Seism. Society American, 63 (1973), pp. 1375–1391</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] J. Zhang; Quadruple-grid velocity-stress finite-difference method for elastic wave-propagation simulation; Geophysics J. Int., 131 (1997), pp. 127–134</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] V.M. Cruz-Atienza, J. Virieux; Dynamic rupture simulation of non-planar faultnwith finite-difference aproach, Geophys; J. Int., 158 (2004), pp. 939–954</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] D.M. Boore; Love waves in non uniform waveguides: finite-difference calculations; J. Geophys. Rev, 75 (1970), pp. 1512–1527</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] T. Mikumo,  Miyatake; B.A. Bott (Ed.), Numerical modeling of realistic fault rupture process. Seismic motion synthetics, Academic Press (1987), pp. 91–151</li>
<li><span id='bib0100'></span>
[[#bib0100|[20]]] A. Pitarka; 3-D elastic finite-difference modeling of seismic motion using staggered grids with non uniform spacing, Bull; Seism. Soc. Am., 89 (1999), pp. 54–68</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] Y. Wang, J. Xu and G.T. Schuster, Viscoelastic wave simulation in basins by a variable-grid finite-difference method, Bull. Seism. Soc. Am. 91, 1641-1749 (2001Bull. Seism. Soc. Am. 90, 587-603 (2000).).</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] L. Lucy; A numerical approach to testing the fission hypotesis; App. Jour., 82 (1977), pp. 1013–1024</li>
<li><span id='bib0115'></span>
[[#bib0115|[23]]] J.J. Monaghan; Why particle methods work?, SIAM Scientific; Journal of stat. Comp., 3 (1982), pp. 422–433</li>
<li><span id='bib0120'></span>
[[#bib0120|[24]]] J.J. Monaghan; And introduction to SPH; Comput. Phys. Comun., 48 (1988), pp. 89–96</li>
<li><span id='bib0125'></span>
[[#bib0125|[25]]] P. Lancaster y, K. Salkauskas; Surfaces generated by moving least square methods; Mathematics of computation, 37 (1981), pp. 141–158</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] B. Nayroles, G. Touzpt y, P. Villon; Generalizing the finite element method: Diffuse approximation and diffuse elements; Computational mechanics, Springer-Verlag, 30 (1992), pp. 7–318</li>
<li><span id='bib0135'></span>
[[#bib0135|[27]]] T. Belytschko, Y.Y. Lu y, L. Gu; Element-free Galerkin methods; International journal for numerical methods in engineering, 37 (1994), pp. 229–256</li>
<li><span id='bib0140'></span>
[[#bib0140|[28]]] C.A. Duarte y, J.T. Oden; A meshless method to solve boundary-value problems; TICAM report, 9 (1995), pp. 5–105</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] I. Babuska and J.M. Melenk, The partition of unity method, International journal for numerical methods in engineering, vol 40, 727-758. Jhon Wiley & Sons, 1995.</li>
<li><span id='bib0150'></span>
[[#bib0150|[30]]] P.S. Jensen; Finite difference techniques for variable grids; Computers & Structures, 2 (1972), pp. 17–29</li>
<li><span id='bib0155'></span>
[[#bib0155|[31]]] N. Perrone, R. Kao; A general finite difference method for arbitrary meshes; Computers & Structures, 5 (1975), pp. 45–58</li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] T. Liszka y J. Orkisz, The finite difference method at arbitrary irregular grids and its application in applied mechanics, Computers & Structures 11, 83-95. Pergamon Press 1980.</li>
<li><span id='bib0165'></span>
[[#bib0165|[33]]] J. Orkisz, Meshless finite difference method II. Adaptative approach, Computation Mechanics, CIMNE, 1998.</li>
<li><span id='bib0170'></span>
[[#bib0170|[34]]] J.J. Benito, F. Ureñy, L. Gavete; Influence several factors in the generalized finite difference method; Applied Mathematical Modeling pp., 25 (2001), pp. 1039–1053</li>
<li><span id='bib0175'></span>
[[#bib0175|[35]]] J.J. Benito, F. Ureñ, L. Gavete y, R. Alvarez; An h-adaptive method in the generalized finite difference; Comput. Methods Appl. Mech. Eng., 192 (2003), pp. 735–759</li>
<li><span id='bib0180'></span>
[[#bib0180|[36]]] J.J. Benito, F. Ureñ, L. Gavete y, B. Alonso; Application of the Generalized Finite Difference Method to improve the approximated solution of pdes; Computer Modelling in Engineering & Sciences, 38 (2009), pp. 39–58</li>
<li><span id='bib0185'></span>
[[#bib0185|[37]]] J.J. Benito, F. Ureñy, L. Gavete; Leading-Edge Applied Mathematical Modelling Research (chapter 7); Nova Science Publishers, New York (2008)</li>
<li><span id='bib0190'></span>
[[#bib0190|[38]]] J.J. Benito, F. Ureñ, L. Gavete y, B. Alonso; Solving parabolic and hyperbolic equations by Generalized Finite Difference Method; Journal of Computational and Applied Mathematics, 209 (2) (2007), pp. 208–233</li>
<li><span id='bib0195'></span>
[[#bib0195|[39]]] F. Ureñ, J.J. Benito, L. Gavete, R. Alvarez, Resolución de ecuaciones diferenciales en derivadas parciales dependientes del tiempo de segundo orden utilizando Diferencias Finitas Generalizadas. Revista Internacional de Métodos Numéricos para Cálculo y Diseñen Ingeniería. Vol. 19, 3, 331-340 (2003).</li>
<li><span id='bib0200'></span>
[[#bib0200|[40]]] P. Moczo, J. Kristek, M. Galis, P. Pazak, M. Balazokjech; The finite difference and finite element modelling of seismic wave propagation and earthquake moltion; Acta physica slovaca, 57 (2) (2007), pp. 177–406</li>
</ol>