Presentamos un ejemplo sencillo en el que la técnica variacional de los elementos finitos (caso particular de modelo pseudoespectral) ha sido aplicada a la resolución de la ecuación unidimensional de las
ondas de gravedad

García Rendo, Jesús

Hortal, Mariano

Spanish

Hortal Reymundo, Mariano

Agencia Estatal de Meteorología

- 69 -TRATAMIENTO NUMERICO POR TECNICAS PSEUDOESPECTRALES (ELEMENTOS Fl NITOS) DE LA ECUACION UNIDIMENSIONAL DE LAS ONDAS DE GRAVEDAD RESUMEN Tras la introducción que hicimos de los modelos pseudoespectrales en el anterior número de esta revista (por cierto que lamen-tamos la involuntaria substitución de Legen-dre por Lagrange ) presentamos un ejemplo sencillo en el que la técnica variacional de los elementos finitos (caso particular de modelo pseudoespectral ) ha sido aplicada a la resolu-ción de la ecuación unidimensional de las ondas de gravedad. En este ejemplo se utilizan diferentes funciones básicas para obtener una solución aproximada, con objeto de apreciar la influen-cia de la regularidad de las mismas en el grado de aproximación y estabilidad de la solución. La solución exacta pertenece, en cada instante tk, al espacio de Sobolev (Hilbert ) H 1 JO ,L[ y por tanto haremos la compara-ción entre la solución exacta y las soluciones aproximadas mediante la norma de H 1 aplica-da en diferentes instan tes del tiempo. El resultado es que, utilizando como ba-ses las funciones "chapeau" o las "exponen-ciales", se alcanza mayor estabilidad que utili-zando las funciones " características" debido a que aquéllas son más regulares que éstas. por J. G. Rendo y M. Hortal Equipo de Predicción Numérica deii.N.M. 1.- Modelo La propagación en una dimensión (x) de una onda de gravedad hasta el instante t = T viene descrita en el dominio [O ,L] x [O ,T] = (O ~ x ~ L; O ~ t ~ T ) por el sistema lineal de ecuaciones en derivadas parciales: (la) a u ah --+g - . =0 a t a x ( 1 b) ah a u -+H--=0 at a x ( le) u (x,O ) = IP(x) ; h(x,O) = H'l' (x) (condiciones iniciales ) (ld) u(L,t) = u(O,t) =O (condiciones de contorno de pared rígida) donde u es la velocidad de movimiento del fluido en la dirección x, h su altura, g la ace-leración de la gravedad y H una constante que coincide con la media espacial de h. Las funciones lf?(x ) y 'l'(x ) son conoci-das y de cuadrado integrable sobre el dominio ] O, L [. - 70 -2.- Espacio de definición de u y h Se denomina L 2 ( ]O ,s [ ) al espacio de las funciones numéricas definidas sobre el intervalo abierto ] O,s [ cuyo cuadrado es integrable sobre este intervalo. Este espacio es de Hi.lbert para el producido escalar ·S (f,g) ~jo f(x)g( x)dx ; f,g ,L2 ( ]O,s[ ) Sea si a H 1 ( ]O,s[ ) le dotamos de producto escalar ((f,g)) ~ ¡;t(x)g(x) + ~~ i;! ) dx resulta un espacio de Hi.lbert Sea ahora (2) E= L 2 [O,T; H 1 ( ]O,L[ )] el espacio de las funciones de fin idas sobre ] O ,T[ de cuadrado integrable sobre este inter-valo y de valores en H 1 (]O,L[). Definimos en E el producto escalar (( f,g)J. ~[~( f,g) )dt con lo que E resulta un espacio de Hilbert separable y por tanto admite bases libres . En este espacio tomaremos las funciones u y h del modelo. 3 .- Téc nicas de aproximac ió n 3.1.- Espacios aproximantes de E Introducimos como primeros aproximan-tes de E los espacios producto senso rial que son subespacios prehilbertianos de E Una función ~(x,t) que pertenezca a E 1 o a E2 pertenece a E y puede escribirse como ( 3) donde ( 4) es una base de E 1 o de E2 y Aik es una su-cesión de números reales que cumple la condición ¿ A2 < oo j, k J k Si truncamos · el desarrollo ( 3) para O ~ j ~ N 1 ; O ~ K ~ N 2 (N 1 puede ser igual a N2 y ambos son enteros positivos ) tenemos una estimación~ N de ~ dada por N 2 N 1 (4b ) ~N =k~o ( i~o Ajkwi (x))ek(t) El espacio de las combinaciones lineales de las funciones ( 4) es un aproximan te de E 1 (o de E2 ) y por tanto un segundo aproxi-mante de E. Representamos este espacio vectorial por V N ( O,T ; O,L ) 3.2.- Método variacional de los elemen-tos finitos 3.2.1 Aproximación en E 1 • a) Elegimos como base de V N ( O,T;O,L) las funciones (S) - 71 -(es decir, en (4) se hace e k= wk) donde las w son las funciones "chapeau" Para constituir dichas funciones hacemos sendas particiones regulares de los intervalos [O,L] y [O,T] con longitudes (distancia inter-nodal)a = xi- xi-t y b = tk - tk-I respectiva-mente y definimos (6) (x-xi)/a ; xi-l ~ x ~ xj wi(x) = (xi+l -x)/a ; xi~ x =< xi+I O, en el resto del intervalo [O,L] (t - tk-1 )/b; tk_ l {;; t =< tk (7 ) wk (t ) = ( tk+ 1 - t)/b;tk<t~tk+ t O, en el resto del intervalo [O,T ] Si escribimos las ecuaciones (la) y (lb) en forma matricial , tenemos (8 ) o u ~ +AU= O o t "" donde }¿ = {u,h}EExE 5x y o A =(H: &x Multiplicando ( 8) escalarmente en L 2 por g _§___) '!!_E E x E se obtiene la ecuación variacional del modelo que tiene por aproximan te la e cuación {f,g} EVN( O,T ;O,L) X VN (O,T;O,L) donde uN y hN vienen dados por ( 4b) Eligiendo como W N expresiones de la forma{wik wm n} que constituyen . una base del espacio producto cartesiano V N [ O,T ;O,L] X VN [O,T ;O,LJ; (wik Y wmn vienen dados por (S) Llamando en ( 4b) pik a ~k para ~N =uN Qik para ~N= hN La ecuación (10) tiene por expresión J.TJL N () w. 1 (11) [( L P¡ 1 --1- w.k + o o i ,1 o t J ~ O W¡¡ · .. N O W ¡¡ Q¡¡ -8-- wmn ) + (g L Q¡¡-- w.k + i,l t i ,l ox J N . ~ W·¡ H .L P¡ 1- 5-1- wmn ]dtdx = O 1,! X 1 ~ m,n <N 1 ~ j,k =<N que es uh sistema lineal de ecuaciones alge-bráicas en los coeficientes P y Q (N= N 1 = N2 ) denominado esquema de Galerkin. Las condiciones iniciales (le) se adaptan al problema aproximante (11), obligando a c umplir las condiciones de norma mínima a las diferencias ( 12) y las condiciones de contorn o ( Id) se trans-forman en - 72 -Las 2N ecuaciones en las que (n=O; k=O) re -lacionan las cantidades P¡ 0 ,Q¡ 0 ,P¡¡ yQ¡¡ porque la integral de la función producto de 8wil 8w¡ 1 -- w y w . Dt . jO ~ JO es distinta de cero solamente para 1 = O, 1 Estas 2N ecuaciúnes, juntamente con las (12) forman un sistema completo para las 2(N+l) incógnitas P¡ 1 y Qil. Las 2N ecuaciones en las que n= k =v relacionan los coeficientes P. Q. p 'Q p Q l(V-1)' l(V-1)' iV• iV• i(V+l)• i(V+l) porque las únicas integrales distintas de cero son aquellas en que 1 = v-1,v,v+ 1 ' Estas ecuaciones, junto con las ( 12 ), una vez conocidos P¡(v- 1 ), Qi(V-1)' P¡v Y Q¡v . forman un sistema completo para obtener las 2(N + 1) incógnitas Pi(V+l) y Qi(V+t) b) Elegimos como base de VN(O,T;O,L) las funciones (14) wjk(x,t)=ej(x) ~ ek(t)=ej(x)ek(t) donde exp {-(x-xj)2 /[(x-xj) 2 - a2 ]} si xj- l ~ x ~ xH 1 O en el resto del intervalo [O,L] (16) ek (t)= exp {-(t-tk) 2 /[(t-tk) 2 - a2 ]} si t k -l~ t ~ tk+l O en el resto del inte rvalo [O,T] Llevando estas funciones a ( 3) y sustituyendo en (11), (12) y (13) resulta un esquema análogo al del apartado a) e) Siguiendo en el espacio E 1 , definimos que es una base combinada de las dos anterio-res con funciones base "chapeau" en el espa-cio y "exponenciales" en el tiempo, resultan-do un esquema similar a los anteriores. 3.2. 2. Aproximación en E2 Situémonos ahora en el espacio E2 donde definimos wjk (x,t) = wj(x)8 k (t) donde las wj(x) son de nuevo las funciones "chapeau" y8k(t)son las funciones "caracte-rísticas" definidas por (17) 8k(t)J1 si tk-1~ t~ tk+l lo en el resto del intervalo [O,T] la derivación respecto del tiempo para estas funciones se toma en el sentido de las distri-buciones y viene dada por de k ( t) -(e k -vz-e k -1/2) /b dt haciendo las substituciones de estas bases en (3) se obtiene el esquema equivalente al (11), (12), (13) en este espacio. 4.- Discusión y resultados Comparamos la solucion aproximada ob-tenida para(u) y (h ) con la solución exacta en instantes tk (donde K es múltiplo de 5) me-diante la n orma en H 1 ]O ,L[ p ara las fun ciones diferencia (u N-u) y (hN-h). Como en el ins-tante tk los valores de cada fun ción base con índice de tiempo dist'n t o de k son cero, y el - 73 -de la de índice k es 1, el valor de la solución aproximada es simplemente N uN(x,t=tk) = ~ Pikwi(x) J debe observarse, sin embargo, que en los in-tantes tk (nodos de la variable temporal) es donde la solución aproximada difiere más de la solución exacta en el caso de las funciones "chapeau" temporales. Consideremos el caso L = 10500 km y T 30 horas, con las condiciones iniciales u(x,O) = U0 sen (2rrx/L) h(x,O) = H = cte. con U0 = 54,6 m/s y H= 9184 m que son los que toman Wang et al con funcio-nes base cúbicas. La solución analítica en este caso es u(x,t) = U0 sen (2rrx/L)cos(2rrtc/L) h(x,t) = H-(HU0 /e )e os( 2rrx/L )sen( 2rrct/L) donde e= gHf 12 El periodo de esta onda es L/c que en nuestro caso es 35000 s (9,7 horas). La longitud de malla a = xi-xi_ 1 se toma en todos los casos igual a O,l.L. Para el intervalo de tiempos b=tk -tk-l se toman tres valores diferentes : 3600 s (1 ho-ra), 1800 s (30 min) y 600 s (10 min). Para que el número resultan te de la com-paración de uN y u (hN y h) dé una idea intuitiva de la discrepancia relativa entre ambas funciones definimos el error funcional (E .F.) de la siguiente forma 2 /LU2o )t/2 E.F.(u) = ( 11 uN -u 1[ H 1 h h 2 1 LH2 )t/2 E.F.(h)= (11 N - II Hl Las tablas 1,2 y 3 dan los valores del E.F cada 5 horas para los intervalos de 3600 s, 1800 s y 600 s respectivamente. Como puede verse en las tablas , el ajuste inicial en la aproximación de El es mucho mejor con las funciones "chapeau" que con las " exponenciales" porque, con estas últimas como bases, las derivadas en los puntos xi de la función aproximada son necesariamente nu-las; este ajuste puede mejorarse, a expensas de un aumento de cálculo, introduciendo, junta-mente con las "exponenciales" otras funcio-nes base cuya derivada en los puntos xi sea la unidad, tal como hacen Wang et al. ( 197 3); como nuestro propósito fue principalmente estudiar las ventajas de algunos elementos finitos como base y la estabilidad de la solu-ción en un espacio de Hilbert o de Sobolev, no introducimos este mejor ajuste y utiliza-mos principalmente las funciones "chapeau'' en la parte espacial. Las funciones " exponenciales" como ba-se en la dimensión temporal dan, en general, mejor estabilidad en ambos u y h que las otras bases, debido a su mayor regularidad. No obs-tante, para periodos cortos, la mayor simpli-cidad de las funciones "chapeau" las hace pre-feribles a las "exponenciales" por las mismas razones expuestas para la dimensión espacial. En cualquier caso, las integrales pueden calcularse de una vez por todas al comienzo de los cálculos con lo que, en cada intervalo de tiempo, lo único que hay que hacer es re-solver un sistema lineal de ecuaciones alge-bráicas. Utilizando las funciones características 8 k como base de tiempo, no se puede utilizar un paso de tiempo largo (en este caso 1 hora) porque aparece inestabilidad semejan te a la presen tada por las diferencias finit as en esque-mas no Hl trados. No obstan te , un paso de tiempo demasiado pequeño origina acumula-ción de errores de cálculo (redondeo ) y por lo tanto, el tiempo ideal debe estar compren-dido entre esos dos extremos. - 74 -Cuando la aproximación en la dimensión temporal es median te funciones "chapeau" o "exponenciales", el método de Galerkin apli-cado en la deducción del esquema numérico es estable aún para pasos largos de tiempo. En el caso de las funciones "chapeau", la dis-minución del intervalo de tiempo produce acumulación de errores de redondeo, lo que no sucede con las "exponenciales" dada su mayor regularidad. ABST RAC T In this paper we present, as a simple example of pseudospectral models, the application of the finite elements variationa! technique to the solution of the unidimen-sional gravity wave equation. We use different basis functions to obtain an approximate solution in order to estímate the influence of their regularíty on the stability and degree of approximation of this solution. The exact solution belongs, in every instant tk, to the Sobolev (Hilbert) space H1 ]O,L[ and therefore the comparison bet -ween the exact ant the app roximate solutions is made by means of the H 1 norm at different time points. The result obtained is that , using as basis functions piecewise linear or "exponential" functions , a greater stahility is obtained than using ''characteristic" functions due to their greater regulatity . BIBLIOG RA FIA Bourke, W., ''A muti-level spectral model. l. Formulation and hemispheric integrations Mo. Wea. Rev. 102 pp. 687-701. (1974). Cea, J., "Approximation variationnelle des problemes aux limites Ann ". Ins. Fourier 14,345-444. (1964) Ciarlet, D.G. and Raviart, P.A. "General Lagrange and Hermite interpolation en Rn with applications to finite element methds ". Arch. Rational Mech. Anal 46, 177-199. (1972). Khoan, YO-K. "Distributions Analyse de Fourier Operateurs aux limites derivées parcielles". tome 1,1·289 et tome 2,1 -332. Vuibert ; París. (1970). Lions, J.L. "Equations differentielles opéra-tionnelles ". Springe-Verlag, 1-291 (1961). Lions, J .L. et Maagenes, E. "Problemes aux limites non homogeneos et applications ". Dunod; Paris, 1-372. (1968). Mikhlin, S.G. "Numerical realizations of the variational methds ". Traslated from the Russian by DR.R.S. Andersen; Computer Centre, the Australian National University, Camberra, 1-373. (1966). Raviart, P.A. ''Sur l 'approximation de cer-teines équations d'evolution lineaires". J our. de Math. pures et appliqués, 46, 11-107. (1967). Schwartz, L. "Théorie des Distributions ". Hermann, París; 1-413. (1966). Teman, R. "Analyse Numerique ". Presses Universitaires de France; Paris, 1-11.(1970) Treves, F. ' Topological vector spaces, distri-butions and Kemels" Pure and applied Mathematics; 96, l-565. (1967). Wang, H.H.; Halpern , P., Douglas Jr., J., and Dupont, T. "Numerical solutions of the Une-Dimensional Primitive Equations using Galerk in App roximations with localized basis fu nc tions ". Mon. Wea. Rev. 100, 738-746.(1973). Y osida, K. "Fun ctional Analysis ". Springer-Verlag. Secod Edition , 1-465. (1968). - 75 -Funciones básicas Tiempo/ h o 5 10 15 20 25 30 F.E. wi (x)wk (t) u 1 ,4. 1 o -2 2.1 o -2 1 ,4.1 o-2 2.1 o-2 6.1 o-2 1,2.10-1 1,5.10-1 h 2,6.10-3 2,8.1 o -2 6.1 o -2 8 .1 o -2 1.1 o-1 1,2.10-l 1,5.10-1 u 8.1 o-2 1.1 o-1 1,8.10-1 2,1. 1Q-l 1,4 .10-l 7. 1 o-2 8.1 o-2 ei(x)ek (t) h 1 ,3.1 o -l 3,7.10-1 1,3.1 0-1 3,7 .1 o-1 1,4.10-1 3,8.1 o-2 1,4.10-1 u 1,4.10-2 5,6.1 o-2 1,5.10-1 2.0.1 o-1 1,2.10-1 1 ,5.1 o-2 a,o.1 o-2 wl (x)~ (t) h 2,6.1 o-3 3,5. 1 o-l 2,6. 1 o -2 3,5 .1 o-1 4,0 .1 o -3 3,5 .1 o-1 3,0.1 Q-2 u 1 ,4. 1 o -2 3,2 .1 o -1 4,2.101 9,3.1 03 1 ,6 .1 06 3,4.1 08 6,3. 1030 wi (x)ek (t) h 2,6.10 -3 6,7 .1 o -2 8,9.10° 1,5.103 3,2.1 05 5,8.1 07 1 ,2 .1 020 TA BLA 1 Error funcional de la solución aproximada para difere ntes instantes de tiempo y diferentes funciones base para un paso de tiempo de 3. 600 seg. Funclon~ Tiempo/ h o 5 10 15 20 25 30 básicas F . E . u 1,4.10-2 3.1 o-2 6.1 o -2 1, 1.1 o-1 1,3.1Q-l 1,6.10-1 2.1 o-1 wi(x)wk (t) h 2,6.10 -3 5.1 o-2 9,7 .1 o -2 1 ,4.1 o-1 1,8.10-1 2.2.1 o-1 2,6.10-1 u 8.1 o-2 1,9.10-1 2.1 o-1 8.1 o-2 1,6.10-1 1,8.10-1 8,4, 1 o -2 ei (x)~ (t) h 1,3.1Q-l 1,4.10-1 1,4.10 -1 1,4 .1Q-l 1 ,4.1 o-1 1,4.10-1 1,4.10-l u 1,4. 1 o -2 1 ,8. 1 o -z 1,7.10-1 2,4.1 o -2 1 ,5.1 o-1 1,2.10-1 7 ,5.1 o -2 w1 (x)~ (t) h 2,6.10 - 3 1 ,3.1 o -2 2,9.1 o-2 1 ,6.1 o -2 3,5.1 o-3 4,6.1 o-2 2,8.1 o -2 u 1 ,4.1 o -2 3,7 .1 o -2 5,5.10-2 6,3.1 o -2 5, 1.1 o-2 2,9.1 o-2 2,9.1 o-2 w 1 (x )e k ( t) h 2 ,6.1 o -3 5,3.10 -2 1, 1.1 o-1 1,7 .1 o -1 2,3.10-1 2,9.1 o-1 3,6.10-1 TABLA 2 Error funcional de la solución aproximada para diferentes instantes de tiempo y diferentes funciones base para un paso de tiempo 1.800 seg. - 76 -Funciones Tiempo/h o 5 básicas 10 15 20 25 30 F .E . u 1,4.10-2 9,5.10-2 1,8 .1Q-l 2,5 .1 o-1 3. 1o-l 3,4.1 o-1 3,8.1 o-~ w¡(x)wk(t ) h 2,6.1 o-3 1 ,4.1 o-1 2,5.1 o-1 3,5.1 o-1 4,4.1 o-~ 5,2.1 o-1 5,9.10-1 u 8.1 o-2 2.10-1 2.1 o-1 8.1 o-2 1,7.10-l 1,7.10-l 8,5.1 o-2 e1 (x)~ (t) h 1 ,3.1 o-l 1,4.10-l 1,4.10-l 1,4.10-l 1,4.10-1 1 ,4.1 o-~ 1,4.10-l u 1,4.10-l 1,9.10-1 1,7.10-1 2,4.1 o-2 1 ,6.1 o-~ 1,2.10-1 7 ,5.1 o-2 w 1 (x)~ (t) h 2,6.1 o-3 1 ,4.1 o-2 3.1 o-2 1 ,5.1 o-2 3,6. 1 o-j 4,7 .1 o-2 2,7 .1 o-2 u 1,4.1 o-2 1,1. 1Q-l 2,3.1 o-l 3,6. 10-1 5.1 o-1 6,4 .1 o-1 7,7.1Q-l w¡ (x)8k (t) 2,6.1 o-3 h 1,5.10-l 3,3.1 o-~ 5,4 .1 o -1 7,8 .1Q-I 1,1 1,4 TABLA 3 Error funcional de la solución aproximada para diferentes instantes de tiempo y diferentes funciones base para un paso de tiempo de 600 seg. 

Tratamiento numérico por técnicas pseudoespectrales (elementos finitos) de la ecuación unidimensional de las ondas de gravedad

https://repositorio.aemet.es/bitstream/20.500.11765/1153/1/TyC_1983_1_08.pdf

Tratamiento numérico por técnicas pseudoespectrales (elementos finitos) de la ecuación unidimensional de las ondas de gravedad

Abstract

Similar works

Full text

Available Versions

Agencia Estatal de Meteorología