Aproximación numérica de integrales de convolución mediante el método de cuadratura de convolución

Trabajo Fin de Grado Aproximación numérica de integrales de convolución mediante el método de cuadratura de convolución Miguel Picos Maiztegui (2019-2020) UNIVERSIDAD DE SANTIAGO DE COMPOSTELA GRADO DE MATEMÁTICAS Trabajo Fin de Grado Aproximación numérica de integrales de convolución mediante el método de cuadratura de convolución Miguel Picos Maiztegui (2019-2020) UNIVERSIDAD DE SANTIAGO DE COMPOSTELA iii iv Trabajo propuesto Área de Conocimiento: Matemática Aplicada Título: Aproximación numérica de integrales de convolución mediante el método de cuadratura de convolución Tutor: Jerónimo Rodríguez Cotutor: Alfredo Ríos Breve descrición del contenido: En múltiples aplicaciones se necesita calcular integrales de convolución de dos funciones f∗g(x) := Zx 0 f(t−x)g(t) d t. El objetivo de este trabajo es la comprensión y programación del método de cuadratura de convolución de cara a aproximar numéricamente este tipo de integrales. La aproximación de la convolución entre dos funciones f y g sobre una malla se obtiene mediante la convolución discreta con los valores de g sobre la misma malla. Los pesos de cuadratura se determinan mediante la transformada de Laplace de la función f (función llamada con frecuencia el núcleo de convolución), un integrador de Runge-Kutta y la fórmula integral de Cauchy. Una vez se haya comprendido y programado el método para ejemplos sencillos, se tratará de aplicar a la resolución de EDO. Asimismo se estudiará la convergencia del tal aproximación. Breve Planicación 1Recordatorio de conceptos básicos que se usarán. 2Comprensión del método de cuadratura de convolución. 3Programación del mismo en Matlab. 4Aplicación del método al cálculo de integrales de convolución y resolución de EDO sencillas. 5Redacción de la memoria. Índice general Resumen vii Introducción ix 1. Integradores de tipo Runge-Kutta 1 1.1. ProblemadeCauchy............................... 1 1.2. Integradores de tipo Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3. Propiedades.................................... 3 1.3.1. Caracter Bien Planteado . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3.2. Convergencia ............................... 3 1.3.3. A-estabilidad ............................... 4 1.3.4. Rígidamente Preciso . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2. Desarrollo del Método 7 2.1. Reescritura de la Integral de Convolución . . . . . . . . . . . . . . . . . . . 7 2.2. Introducción del método Runge-Kutta . . . . . . . . . . . . . . . . . . . . . 9 2.3. Transformada Z de {xn}n∈N ........................... 10 2.4. Descomposición Lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5. TeoremadeCauchy................................ 17 2.6. Cálculo de la serie de potencias de ˆ f∆(z) h .................. 21 2.7. No rígidamente preciso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3. Resultados Numéricos 27 3.1. Consideraciones Previas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2. TestAcadémico.................................. 30 3.2.1. FunciónEscalón ............................. 30 3.2.2. Polinomios ................................ 32 3.3. EDOsLineales .................................. 37 v vi ÍNDICE GENERAL 3.3.1. EDOs de Primer Orden . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3.2. EDOsdeordenn............................. 40 3.4. EcuacionesIntegrales............................... 47 3.4.1. Ejemplo.................................. 49 Códigos 51 Bibliografía 57 xiv INTRODUCCIÓN Capítulo 1 Integradores de tipo Runge-Kutta Durante este capítulo introduciremos los métodos Runge-Kutta que, como vimos en la introducción, serán una parte principal del desarrollo del método de cuadratura de convolución. Nos centraremos en las propiedades que pueden tener estos integradores, especialmente en aquellas que serán más importantes para el método. Para el estudio teórico se utilizarán como principales referencias los libros [5] y [6] , que exponen de forma extensa todos los componentes de este capítulo. 1.1. Problema de Cauchy Los métodos Runge-Kutta son una familia de integradores que aproximan numéricamente la solución de un problema de Cauchy. Llamamos problema de Cauchy al problema cuya incognita es la solución de una ecuación diferencial que cumpla una cierta condición inicial. En este trabajo solo afrontaremos ecuaciones diferenciales de primer orden en una dimensión espacial. Denición 1.1 (Problema de Cauchy) . Sean f(t, x) una función real de dos variables reales, t0, x0∈R . (x0(t) = f(t, x(t)), x(t0) = x0. donde x(t) es una función derivable denida en un intervalo de R . Debemos poder garantizar la existencia y unicidad de este problema, para ello la función f(t, x) deberá cumplir unas ciertas propiedades. El siguiente resultado nos permitirá armar que los problemas que consideremos tengan solución única. 1 2 CAPÍTULO 1. INTEGRADORES DE TIPO RUNGE-KUTTA Teorema 1.1. Sea f: [a, b]×R→R una función real de dos variables reales tal que f∈C1([a, b]×R,R) . Entonces, el problema de Cauchy denido en 1.1 tiene solución única denida en el intervalo [a, b] Observación 1.1. El desarrollo anterior puede ser generalizado para dimensiones mayores, tal que f: [a, b]×Rn→Rm , pero para este trabajo no será necesario. 1.2. Integradores de tipo Runge-Kutta Los métodos numéricos para la resolución de problemas de Cauchy son útiles para calcular aproximaciones de soluciones que no pueden ser obtenidad analíticamente. En el caso de los métodos Runge-Kutta, se discretiza el dominio del problema mediante una malla uniforme con un paso de discretización h . Esta malla estará formada por N+ 1 nodos de discretización de modo que estos puntos dividan el intervalo [a, b] en N intervalos de longitud h . Estos métodos se caracterizan por la obtención de un mayor orden de aproximación gracias a la evaluación de la función f en varios puntos intermedios entre los nodos. Dados un intervalo [a, b]⊂R y {tn}0≤n≤N los nodos equidistantes de la malla sobre el dominio de la función f del problema de Cauchy (1.1). Se dene un método Runge-Kutta de m etapas como: xn,i =xn+h m X j=1 aijf(tn,j, xn,j) , i= 1,··· , m xn+1 =xn+h m X i=1 bif(tn,i, xn,i) (1.1) donde tn,j =tn+hcjj= 1,··· , m y los coecientes aij bi y cj están dados. Los coecientes son los que denen un método particular. Estos suelen expresarse en una tabla de Butcher: c A bT donde c , b∈Rm y A= (aij)m i,j=1 ∈Mmxm(R) . Si denimos Xn:= (xn,1,··· , xn,m)T , fn:= (f(tn,1, xn,1),··· , f(tn,m, xn,m)) y −→ 1 := (1,··· ,1)T , podemos reformular (1.1) del siguiente modo: Xn=xn−→ 1 + hAfn xn+1 =xn+hbTfn (1.2) Observación 1.2. Esta última expresión nos será útil durante el trabajo ya que la función de la ecuación diferencial que vamos a considerar es f(t, x) = sx +g(t) de modo que fn=sXn+gn , siendo g(tn+cih) . 1.3. PROPIEDADES 3 1.3. Propiedades Estudiaremos ahora las propiedades que pueden tener los integradores descritos en la sección anterior. Comenzaremos por propiedades básicas como que el carácter bien planteado del esquema numérico o sea convergente y tratando también otras propiedades más particulares que serán especialmente importantes para el método de cuadratura de convolución como puede ser la A-estabilidad. 1.3.1. Caracter Bien Planteado Para comprobar que el método esta bien denido debemos comprobar que el sistema de ecuaciones que tenemos que resolver en (1.1) para obtener el vector Xn tiene solución. Para un h sucientemente pequeño podremos asegurar que el sistema tenga solución. Proposición 1.1. Si h < 1 LM , con M=maxi=1,···,m Pm j=1 |aij| y L constante de Lipchitz de f con respecto de x, entonces el sistema de (1.1) que dene Xn tiene solución única. 1.3.2. Convergencia Para que un método sea útil, tiene que ser convergente, es decir, las aproximaciones de la solución del Problema de Cauchy con el que estemos trabajando deben acercarse más a la solución exacta a medida que disminuye el parámetro de discretización del problema. Los métodos Runge-Kutta son métodos de un paso, por lo que adoptan la forma: (xn+1 =xn+hΦ(xn+1, xn, tn;h) x0 dado (1.3) Podremos asegurar la convergencia siempre que los métodos sean estables y consistentes. Estabilidad El polinomio de estabilidad de los métodos Runge-Kutta es ρ(r) = r−1 por ser un método de un paso de la forma mostrada en (1.3). En consecuencia, como su única raíz es 1 es un polinomio de Dahlquist (todas sus raíces son de módulo menor o igual que 1 y las de modulo 1 son simples). Por esta razón, todos los métodos Runge-Kutta son estables. Para entrar en profundidad en este resultado se puede consultar [5]. Consistencia Como podemos observar en la expresión (1.1), para el cálculo de la aproximación en el instante tn+1 realizamos una combinación lineal de los valores f(tn,i, xn,i) de modo que 4 CAPÍTULO 1. INTEGRADORES DE TIPO RUNGE-KUTTA querremos que esta combinación sea una media ponderada. Así, a medida que el parámetro de discretización tienda a 0, la combinación lineal tenderá al valor de f(t, x) . Podemos ver, intuitivamente, la condición para que un método Runge-Kutta sea consistente. Proposición 1.2. Un método Runge-Kutta es consistente si y solo si: m X i=1 bi= 1 Uno de los motivos principales por los que utilizamos integradores Runge-Kutta es porque podremos alcanzar órdenes de convergencia mayores que 2 en casos donde necesitamos métodos A-estables. Por esta razón necesitaremos unos criterios que nos permitan dictaminar cuando un método tiene un cierto orden de convergencia. Este estudio es bastante complejo, se puede ver con detalle en [5]. Generalmente se pide la condición de la para los métodos Runge-Kutta, ya que esta propiedad relaja las condiciones orden. Denición 1.2 (Condición de la) . Denimos la condición de la como: X j aij =ci para todo i∈ {1,··· , m} Esto se puede entender como que las etapas intermedias del método sean una aproximación de orden 1 de la solución de la EDO a aproximar. 1.3.3. A-estabilidad Podemos entender la A-estabilidad de un método como la capacidad de este de aproximar una EDO de la forma: x0(t) = λx(t) para todo λ∈C tal que Re(λ)≤0 . La solución exacta de la EDO anterior es x(t) = eλt . Para valores de λ con parte real negativa la solución decae a 0 rápidamente. Consideraremos que un cierto valor λ pertenece a la región de estabilidad de un método si las aproximación de la solución de la EDO se mantiene acotada a medida que aumentamos las interaciones. Denición 1.3 (A-estabilidad) . Diremos que un método Runge-Kutta es A-estable si la región de estabilidad del método contiene al conjunto {z∈C tal que Re(z)<0} . Esta característica va a ser muy importante para el método de cuadratura de convolución, pues la ecuación diferencial que aproximaremos es una EDO de esta forma. Por esta razón se introducen los métodos Runge-Kutta, ya que, con métodos lineales multipaso 1.3. PROPIEDADES 5 solo podemos obtener orden de aproximación 2 manteniendo la A-estabilidad, mientras que con los métodos Runge-Kutta podemos mantener la A-estabilidad con métodos de orden arbitrariamente grande. Para conocer la región de estabilidad de un método se calcula su función de estabilidad, que en el caso de los integradores de tipo Runge-Kutta tiene la forma: R(z) = 1 + zbt(I−zA)−1−→ 1 donde −→ 1 := (1,··· ,1)T . Esta función corresponde al crecimiento de la aproximación de una etapa a la siguiente, ya que: yn+1 =R(z)yn Con el uso de esta función podemos obtener el siguiente resultado, que nos permitirá clasicar los métodos según sean A-estable o no. Proposición 1.3. Un método Runge-Kutta es A-estable si: 1. |R(z)| ≤ 1 para todo z∈C tal que Re(z)≤0 . 2. I−zA es no singular para todo z∈C tal que Re(z)≤0 . 1.3.4. Rígidamente Preciso Introducimos ahora esta propiedad, ya que, aunque no sea estrictamente necesaria, será una de las condiciones que busquemos para los métodos Runge-Kutta que utilizaremos en el método de cuadratura de convolución. Si exigimos la condición de la, podemos considerar el vector Xn como las aproximaciones de la solución de la EDO en los tiempos tn+cih . Con esta consideración podemos entender la propiedad de un método de ser rígidamente preciso como que la última componente del vector Xn sea la aproximación de la solución en el siguiente nodo de discretización. Es decir, que la última coordenada de Xn es igual que xn+1 . Esto conlleva que la última la de la matriz A tiene que ser igual al vector b . Denición 1.4 (Rígidamente Preciso) . Un método Runge-Kutta es rígidamente preciso si el vector b coincide con la última la de la matriz A o, lo que es lo mismo, que se cumpla la condición: bTA−1= (0,··· ,0,1) 6 CAPÍTULO 1. INTEGRADORES DE TIPO RUNGE-KUTTA Capítulo 2 Desarrollo del Método de Cuadratura de Convolución En este capítulo explicaremos la justicación teórica del método de cuadratura de convolución. Afrontamos la aproximación de una integral de convolución, es decir, queremos aproximar la función y(t) : y(t)=[f∗g](t) = Zt 0 f(t−τ)g(τ)dτ (2.1) donde la función g(t) y ˆ f(s) , Transformada de Laplace de f(t) , son conocidas. Para resolver este problema, utilizaremos las relaciones entre la Transformada de Laplace y su inversa. Después, identicaremos una expresión dentro de la integral con la solución de una ecuación diferencial conocida. Posteriormente, realizaremos la discretización del problema aproximando la solución de la ecuación diferencial anterior con un método Runge-Kutta. Finalmente, utilizaremos la transformada Z para, una vez calculadas las integrales utilizando el Teorema de Cauchy, podamos igualar términos de las series de potencias resultantes resolviendo asi el problema. El desarrollo de este capítulo sigue el esquema del artículo de referencia [4]. 2.1. Reescritura de la Integral de Convolución En esta sección, buscaremos reescribir la expresión de la integral de convolución (2.1), para llegar a un problema en función de datos conocidos. Para ello comenzaremos sustituyendo la función f(t) , que puede ser conocida o no, por su expresión en función de su Transformada de Laplace. Podemos obtener el valor de la función f(t) mediante las relaciones entre la transformada de Laplace y su transformada 7 8 CAPÍTULO 2. DESARROLLO DEL MÉTODO inversa, en este trabajo no se profundizará en el uso de la Transformada de Laplace, puede obtenerse mas información en la referencia [7]. f(t) = 1 2πiZc+i∞ c−i∞ est ˆ f(s)ds (2.2) siendo c∈R tal que todos los polos de la función ˆ f(s) tienen parte real menor que c. Por tanto podemos sustituir en la integral inicial, la función f(t−τ) por su correspondiente expresión en función de ˆ f(s) dada por la fórmula de la Transformada Inversa de Laplace (2.2). y(t) = Zt 0 g(τ)1 2πiZc+i∞ c−i∞ es(t−τ)ˆ f(s)dsdτ Necesitaremos realizar suposiciones sobre las funciones ˆ f(s) y g(t) para utilizar el método de cuadratura de convolución. Para este trabajo, asumiremos las siguientes propiedades: 1. La función g(t) es continua. 2. La función ˆ f(s) es analítica. 3. La función ˆ f(s) es la Transfromada de Laplace de una función f(t) . Con estas suposiciones tenemos que las funciónes ˆ f(s) y g(t) son continuas por lo que la función g(τ)es(t−τ)ˆ f(s) es continua y por tanto integrable en compactos. Por otra parte, por ser ˆ f(s) una Transformada de Laplace tendremos que la función g(τ)es(t−τ)ˆ f(s) será integrable en cada sección jando un valor de τ . En conclusión, podremos aplicar el teorema de Fubini, intercambiando el orden de las integrales: y(t) = 1 2πiZc+i∞ c−i∞ ˆ f(s)est Zt 0 e−sτ g(τ)dτds (2.3) Consideremos ahora el siguiente problema de Cauchy: (x0(t) = sx(t) + g(t) x(0) = 0 (2.4) Podemos comprobar que x(t)=est Rt 0e−sτ g(τ)dτ es la solución de este problema. En efecto, (x0(t) = sest Rt 0e−sτ g(τ)+este−stg(t) = sx(t) + g(t) x(0) = 0 En conclusión x(t) es la única solución del problema de Cauchy (2.4). Por lo tanto podemos introducir la solución en la expresión (2.3), reescribiendo la expresión de la integral: y(t) = 1 2πiZc+i∞ c−i∞ ˆ f(s)x(t)ds (2.5) 2.2. INTRODUCCIÓN DEL MÉTODO RUNGE-KUTTA 9 Para el cálculo de esta integral, realizaremos la aproximación de la solución x(t) mediante un método numérico Runge-Kutta. Para, de este modo, poder evaluar la función y(t) . 2.2. Introducción del método Runge-Kutta Como la obtención analítica de soluciones no es siempre posible, para continuar con el método realizaremos una discretización del problema (2.5). En este caso utilizaremos un mallado uniforme del dominio. La función y(t) , que queremos aproximar, esta denida en [0,∞) . Sea h∈R+ el parámetro de discretización elegido para el método. Realizamos un mallado uniforme del dominio de la función y(t) , de modo que {t∈[0,∞) tal que t=ti=ih con i∈N} es el conjunto de nodos de discretización de la malla. De este modo tenemos que, evaluando en los nodos la expresión (2.5), obtenemos: y(tn) = 1 2πiZc+i∞ c−i∞ ˆ f(s)x(tn)ds para todo n∈N (2.6) Sean yn≈y(tn) y xn≈x(tn) con n∈N las aproximaciones de las funciones y(t) y x(t) en los nodos de discretización. Los valores yn son los que buscamos aproximar, mientras que los valores xn son aproximaciones obtenidas utilizando un método Runge-Kutta para resolver el problema de Cauchy (2.4). Con estas aproximaciones llegamos a: yn=1 2πiZc+i∞ c−i∞ ˆ f(s)xnds para todo n∈N (2.7) Utilizaremos un método Runge-Kutta de m etapas, denido por su tabla de Butcher: c A bT donde c , b∈Rm y A= (ai,j)m i,j=1 ∈Mmxm(R) . Necesitaremos que este método cumpla unas ciertas propiedades para que el método de cuadratura de convolución sea convergente: 1. El método Runge-Kutta es A-estable. 2. La matriz de coecientes de método, A , es no singular. 3. El método cumple que bTA−1−→ 1 = 1 . Observación 2.1. Una condición suciente para obtener la propiedad de bTA−1−→ 1 = 1 , es que el método Runge-Kutta sea rígidamente preciso. De este modo tendremos que bTA−1= 16 CAPÍTULO 2. DESARROLLO DEL MÉTODO donde Jλi es el i-ésimo bloque de una forma canónica de Jordan J . Para el cálculo de la inversa de un bloque de Jordan, de tamaño k, utilizaremos el método de Gauss. Jλ=       λ1 ...... λ1 λ               λ1 1 ......... λ1 1 λ1        1 λ∗Jλ −−−→       11 λ 1 λ ......... 11 λ 1 λ 11 λ        Fk−1−1 λFk . . . F1−1 λF2 −−−−−−−−−−−→       11 λ−1 λ2 ......... 11 λ−1 λ2 11 λ        De modo que la inversa de un bloque de Jordan es: J−1 λ=       1 λ−1 λ2 ...... 1 λ−1 λ2 1 λ        En conclusión, calculando el caso particular con el que estamos tratando, obtenemos: (Jz,h −sI)−1=    J−1 λ1−s ... J−1 λl−s     (2.22) donde {λ1,··· , λl} son los autovalores de la matriz ∆(z) h , asociados a los diferentes bloques que forman Jz,h . Como vimos para el caso general, tenemos que para r∈ {1,2,··· , l} : J−1 λr−s=       1 λr−s−1 (λr−s)2 ...... 1 λr−s−1 (λr−s)2 1 λr−s        (2.23) 2.5. TEOREMA DE CAUCHY 17 Finalmente podemos armar, por (2.22) y (2.23) que : (Jz,h −sI)−1kj =       1 λr−s si j=k −1 (λr−s)2 si j=k+ 1 y [J]kk+1 = 1 0 en otro caso (2.24) Siendo en ambos casos λr∈ {λ1,··· , λl} el autovalor asociado con la caja en la que se encuentra la posición kj en la matriz Jz,h . Por lo tanto, utilizando la expresión de (2.21) en la igualdad (2.18), obtenemos: ∞ X n=0 yn+1zn=∞ X n=0 1 2πiZc+i∞ c−i∞ ˆ f(s)bTA−1∆(z) h−sI−1 gnds!zn =∞ X n=0   1 2πiZc+i∞ c−i∞ ˆ f(s)  m X k=0 m X j=0 βkαjh(J−sI)−1ikj ds zn =∞ X n=0 m X k=0 m X j=0 βkαj1 2πiZc+i∞ c−i∞ ˆ f(s)h(J−sI)−1ikj dszn (2.25) Observación 2.8. En término general, la matriz ∆(z) h no tendrá autovalores multiples y por tanto será diagonalizable, de modo que la matriz Jz,h será una matriz diagonal. Por lo tanto obtendremos que: h(J−sI)−1ikj =(1 λk−s Si k=j 0 Si k6=j Para el cálculo de la integrales de la expresión (2.25) buscaremos utilizar el Teorema de Cauchy. De este modo llegaremos a una expresión donde no tendremos que calcular ninguna integral. 2.5. Teorema de Cauchy Buscamos calcular el valor de integrales del tipo: 1 2πiZc+i∞ c−i∞ ˆ f(s) s−λds 1 2πiZc+i∞ c−i∞ ˆ f(s) (s−λ)2ds que son el tipo de integrales que tenemos en la expresión (2.25). Para ello utilizaremos las siguientes curvas, para cada n∈N sean: φn,1:t∈[−n, n]→C , tal que φn,1(t) = c+ it 18 CAPÍTULO 2. DESARROLLO DEL MÉTODO φn,2:t∈h−π 2,π 2i→C , tal que φn,2(t) = c+neit φn la curva cerrada formada por las curvas −φn,1 y φn,2 . Observación 2.9. Para que las curvas φn tengan orientación positiva debemos invertir la orientación de la curva φn,1 . En la gura 2.5 podemos ver como cualquier número complejo con parte real mayor que c acabará encontrándose dentro de la sucesión de curvas considerada. En la gura se muestran las 5 primeras curvas de la sucesión en el caso de c= 10 . Figura 2.1: En esta gura podemos ver las curvas φn para n entre 1 y 5 y c= 10 . Utilizando el teorema de Cauchy sobre las curvas cerradas φn obtenemos: ˆ f(λ) = 1 2πiZφn ˆ f(s) s−λds ˆ f0(λ) = 1 2πiZφn ˆ f(s) (s−λ)2ds (2.26) para todo λ∈C que este en el interior de la curva cerrada φn . Observación 2.10. En efecto las funciones integradas cumplen las condiciones del Teorema de Cauchy pues ˆ f(s) es diferenciable y, en el interior de φn , ˆ f(s) no tiene polos por denición de la transformada inversa de Laplace. Además estamos integrando sobre una región convexa {x∈C tal que Re(x)≥c} ⊂ {z∈C tal que Re(z)>0} abierto convexo. Para más detalle sobre este tema se puede ver la referencia [8]. 2.5. TEOREMA DE CAUCHY 19 Podemos dividir las integrales sobre φn en dos integrales, sobre φn,1 y φn,2 . De este modo tenemos que: l´ım n→∞ 1 2πiZφn,1 ˆ f(s) s−λds= l´ım n→∞ 1 2πiZc+in c−in ˆ f(s) s−λds=1 2πiZc+i∞ c−i∞ ˆ f(s) s−λds l´ım n→∞ 1 2πiZφn,1 ˆ f(s) (s−λ)2ds= l´ım n→∞ 1 2πiZc+in (c−in)2 ˆ f(s) s−λds=1 2πiZc+i∞ c−i∞ ˆ f(s) (s−λ)2ds (2.27) Calcularemos ahora el valor del límite de las integrales sobre φn,2(t) . l´ım n→∞Zφn,2 ˆ f(s) s−λds= l´ım n→∞Zπ 2 −π 2 ˆ f(c+neit) c+neit−λdt l´ım n→∞Zφn,2 ˆ f(s) (s−λ)2ds= l´ım n→∞Zπ 2 −π 2 ˆ f(c+neit) (c+neit−λ)2dt Por ser ˆ f(s) transformada de Laplace se tiene que l´ımRe(s)→∞ ˆ f(s) = 0 y por tanto existe M∈R tal que |ˆ f(s)|<=M para todo s tal que Re(s)> c . Con lo que obtenemos: Zπ 2 −π 2 ˆ f(c+neit) c+neit−λdt≤Zπ 2 −π 2 |ˆ f(c+neit)| |c+neit−λ|dt≤Zπ 2 −π 2 M ndt=πM n Zπ 2 −π 2 ˆ f(c+neit) (c+neit−λ)2dt≤Zπ 2 −π 2 |ˆ f(c+neit)| |c+neit−λ|2dt≤Zπ 2 −π 2 M n2dt=πM n2 Observación 2.11. Para el cálculo anterior hemos supuesto que |ˆ f(s)|<=M para todo s tal que Re(s)> c . Esto es posible porque en la denición de la Transformada Inversa de Laplace, tomamos ˆc∈R tal que todos los polos de la función ˆ f(s) tengan parte real menor que ˆc . Dado que la transformada de Laplace tiene límite 0 cuando la parte real de s tiende a innito, existirá un real ˜c∈R tal que |ˆ f(s)|<=M para todo s tal que Re(s)>˜c . De modo que podemos tomar c=max{ˆc, ˜c} , de modo que se cumplen las suposiciones que tomamos en los cálculos anteriores. Con lo cual 0≤l´ım n→∞Zφn,2 ˆ f(s) s−λds≤l´ım n→∞ πM n= 0 0≤l´ım n→∞Zφn,2 ˆ f(s) (s−λ)2ds≤l´ım n→∞ πM n2= 0 Y por tanto: l´ım n→∞Zφn,2 ˆ f(s) s−λds= 0 l´ım n→∞Zφn,2 ˆ f(s) (s−λ)2ds= 0 (2.28) 20 CAPÍTULO 2. DESARROLLO DEL MÉTODO En conclusión, de lo obtenido en (2.26),(2.27) y (2.28), podemos calcular el valor de las siguientes integrales: 1 2πiZc+i∞ c−i∞ ˆ f(s) λ−sds= l´ım n→∞−1 2πiZφn,2 ˆ f(s) s−λds+1 2πiZφn,1 ˆ f(s) s−λds = l´ım n→∞ 1 2πiZφn ˆ f(s) s−λds=ˆ f(λ) 1 2πiZc+i∞ c−i∞ ˆ f(s) (λ−s)2ds= l´ım n→∞ 1 2πiZφn,2 ˆ f(s) (s−λ)2ds−1 2πiZφn,1 ˆ f(s) (s−λ)2ds = l´ım n→∞ 1 2πiZφn ˆ f(s) (s−λ)2ds=−ˆ f0(λ) (2.29) Observación 2.12. Para el desarrollo del cálculo de las integrales de (2.29) debemos suponer que Re(λ)> c . Pues de ser esto así, existirá n∈N tal que λ se encuentre en el interior de φN para todo N∈N tal que N > n . En conclusión, podemos calcular las integrales de la expresión (2.25) utilizando lo obtenido en (2.24) y (2.29). 1 2πiZc+i∞ c−i∞ ˆ f(s)h(J−sI)−1ikj ds= =         1 2πiRc+i∞ c−i∞ ˆ f(s) λr−sds=ˆ f(λr) Si h(J−sI)−1ikj =1 λr−s −1 2πiRc+i∞ c−i∞ ˆ f(s) (λr−s)2ds=ˆ f0(λr) Si h(J−sI)−1ikj =−1 (λr−s)2 0 en otro caso (2.30) Observación 2.13. Los cálculos que acabamos de realizar son posibles ya que todos los autovalores de la matriz ∆(z) h tienen parte real mayor que c y por tanto podemos aplicar el Teorema de Cauchy en (2.26). Esto se debe a que los autovalores de ∆(z) tienen parte real positiva y por tanto podemos escoger h sucientemente pequeño como para que se cumpla la condición anterior. Los autovalores de A−1 son de parte real positiva porque son los inversos de los autovalores de A que tienen parte real positiva por ser el método RungeKutta A-estable. Los autovalores de I−z−→ 1bTA−1 son 1 con multiplicidad 1−m y 1−z que tendrá parte real positiva para todo z∈C tal que la parte real de z sea menor que 1 (esto se cumple pues trabajamos con |z|<1 ). Como ∆(z) = A−1(I−z−→ 1bTA−1) , entonces todos sus autovalores son positivos. Puesto que todos los autovalores de ∆(z) son todos de parte real positiva, entonces existirá h∈R tal que los autovalores de ∆(z) h tengan parte real mayor que c , cumpliendo así las condiciones del Teorema de Cauchy. 2.6. CÁLCULO DE LA SERIE DE POTENCIAS DE ˆ F∆(Z) H 21 Puesto que estamos trabajando con matrices, nos será de utilidad poder aplicar la función ˆ f(s) a matrices. Para ello denimos la extensión de la función ˆ f:s∈C→C al conjunto de matrices complejas Mmxm(C) de la siguiente forma: Sea A∈Mmxm(C) con su descomposición canónica de Jordan A=LAJAL−1 A siendo JA la forma canónica de Jordan de A y LA la matriz de cambio de base asociada. Denimos ˆ JA como: hˆ JAikj :=        ˆ f[JA]kj si j=k ˆ f0(λr) si j=k+ 1 y [JA]kk+1 = 1 0 en otro caso Denimos ˆ f:A∈Mmxm(C)→Mmxm(C) tal que ˆ f(A) := LAˆ JAL−1 A . Partiendo de (2.25) y utilizando (2.30), (2.20) y la nueva denición de ˆ f(A) obtenemos: ∞ X n=0 yn+1zn=∞ X n=0 m X k=0 m X j=0 βkαj1 2πiZc+i∞ c−i∞ ˆ f(s)h(J−sI)−1ikj dszn =∞ X n=0 m X k=0 m X j=0 βkαjhˆ J∆(z) hikj zn =∞ X n=0 m X k=0 βkek!T ˆ J∆(z) h  m X j=0 αjej zn =∞ X n=0 bTA−1L∆(z) h ˆ J∆(z) h L−1 ∆(z) h gnzn =∞ X n=0 bTA−1ˆ f∆(z) hgnzn (2.31) De este modo hemos conseguido obtener los valores de las integrales de la expresión (2.18) en función de términos que podemos calcular. Para continuar, descompondremos la función ˆ f∆(z) h en su serie de potencias, de modo que podremos transformar la expresión (2.31) en una serie de potencias. 2.6. Cálculo de la serie de potencias de ˆ f∆(z) h Buscamos expresar la función ˆ f∆(z) h en serie de potencias, de modo que podamos operar. ∆(z) h es una matriz que depende únicamente de z , ya que h está jado. Por lo que para cada i, j ∈ {1,2,··· , m} , ˆ f∆(z) hkj :z∈C→C 22 CAPÍTULO 2. DESARROLLO DEL MÉTODO es una función compleja de variable compleja que, debido a la denición de ˆ f∆(z) h , es analítica pues ˆ f(s) es analítica por hipótesis. En conclusión, admite una representación en serie de potencias. ˆ f∆(z) hkj =∞ X n=0 hWh n(ˆ f)ikj zn Por tanto podemos expresar ˆ f∆(z) h como serie de potencias con las matrices Wh n(ˆ f) como coecientes: ˆ f∆(z) h=∞ X n=0 Wh n(ˆ f)zn (2.32) Observación 2.14. La matriz ∆(z) tiene la forma A−1−zA−1−→ 1bTA−1 , como vimos en (2.14). Por lo tanto es una función analítica pues cada componente de la matriz es una combinación lineal de funciones analíticas. De este modo la función ˆ f∆(z) h es analítica (pues ˆ f(s) es analítica). Para resolver el problema necesitaremos conocer los coecientes, Wh n(ˆ f) , de la serie de potencias de ˆ f∆(z) h . Como hˆ f∆(z) hikj es una función analítica y estamos trabajando en una región convexa (un entorno de 0), utilizando el Teorema de Cauchy, tenemos que para todo k, j ∈ {1,2,··· , m} podemos calcular los coecientes de modo que: hWh n(ˆ f)ikj =1 2πiZ|z|=Rhˆ f∆(z) hikj zn+1 dz (2.33) Aproximaremos la integral de (2.33) mediante un método número de integración. Tomando la curva φR:t∈[0,2π]→C tal que φR(t) = Reit y utilizando el método de los trapecios, obtenemos: 1 2πiZ|z|=Rhˆ f∆(z) hikj zn+1 dz=1 2πiZ2π 0hˆ f(∆(Reit) h)ikj (Reit)ndt=R−n 2πiZ2π 0ˆ f(∆(Reit) h)kj e−intdt ≈R−n L L−1 X l=0 ˆ f∆(Reit) hkj Cnldt siendo C=ei2π L Por tanto tomando la misma discretización, que depende de R y L , para cada elemento de la matriz Wh n(ˆ f) podemos aproximar los coecientes de la serie de potencias por: Wh n(ˆ f)≈R−n L L−1 X l=0 ˆ f∆(Reit) hCnl siendo C=ei2π L (2.34) 2.7. NO RÍGIDAMENTE PRECISO 23 Observación 2.15. Utilizamos la regla de los trapecios para esta aproximación por sus buenas propiedades de convergencia. Estamos integrando sobre una curva cerrada regular una función analítica. Por esta razón, el problema se reduce a integrar una función regular y periódica sobre su dominio. En este ámbito, la aproximación tiene una convergencia exponencial, por lo que el orden de convergencia del método será el que aporte el método Runge-Kutta y no el que aporta esta aproximación. Para más información sobre el tema mirar [9]. Finalmente, introduciendo la expresión de ˆ f∆(z) h como serie de potencias, (2.32), en la expresión (2.31) y utilizando la formula de Cauchy para el producto de series innitas, obtenemos: ∞ X n=0 yn+1zn=∞ X n=0 bTA−1ˆ f∆(z) hgnzn =∞ X n=0 bTA−1 ∞ X n=0 Wh n(ˆ f)zn!gnzn =∞ X n=0 bTA−1 n X k=0 Wh n−k(ˆ f)gk!zn En virtud de la unicidad de expresión en serie de potencias, podemos armar que: yn+1 =bTA−1 n X k=0 Wh n−k(ˆ f)gk Wh n(ˆ f)≈R−n L L−1 X l=0 ˆ f∆(Reit) hCnldt siendo C=ei2π L (2.35) Quedando asi completamente denido el método para el cálculo de la integral (2.1). 2.7. No rígidamente preciso En esta sección afrontaremos el mismo problema (2.1) eliminando la condición de que el método Runge-Kutta sea rígidamente preciso. Nos encontraremos con dos casos. El primero es que se cumpla la condición bTA−1−→ 1 = 1 , en cuyo caso el método funcionará igualmente ya que, como menciona la observación (2.2), pedir que el método sea rígidamente preciso es más fuerte que pedir bTA−1−→ 1 = 1 que es realmente la condición necesaria para que el método expuesto anteriormente funcione. El segundo caso es que K=bTA−1−→ 16= 1 . En este caso la expresión (2.17) no sería válida. Para obtener una expresión similar utilizaremos los resultados de las observaciones 24 CAPÍTULO 2. DESARROLLO DEL MÉTODO (2.5) y (2.6) y de la expresión (2.15). ∞ X n=0 xn+1zn=1 KbTA−1∞ X n=0 Xnzn+K−1 KhbT(s∞ X n=0 Xnzn+∞ X n=0 gnzn) =1 KbTA−1∆(z) h−sI−1∞ X n=0 gnzn +K−1 KhbT s∆(z) h−sI−1∞ X n=0 gnzn+∞ X n=0 gnzn! donde ∆(z) = A−1−z 1+(K−1)zA−1−→ 1bTA−1 . Análogamente al caso rígidamente preciso podemos introducir esta expresión en (2.8) obteniendo: ∞ X n=0 yn+1zn=1 K ∞ X n=0 1 2πiZc+i∞ c−i∞ ˆ f(s)bTA−1∆(z) h−sI−1 gn!zn +hK−1 K ∞ X n=0 1 2πiZc+i∞ c−i∞ sˆ f(s)bT∆(z) h−sI−1 gnds!zn +hK−1 K ∞ X n=0 1 2πiZc+i∞ c−i∞ ˆ f(s)dsbTgnzn (2.36) El primer término de la expresión (2.36) es un múltiplo de la expresión a la que se llega si el método es rígidamente preciso. Por este motivo podremos aproximarlo con las fórmulas calculadas anteriormente. Para el segundo sumando, necesitaremos hacer un desarrollo similar al realizado a lo largo de este capítulo, pero con la función sˆ f(s) . De modo que necesitamos que las acotaciones que realizamos en (2.28) también funcionasen para esta nueva función. Para que estas acotaciones sean válidas debemos pedir que el comportamiento asintótico de la función sˆ f(s) sea que este acotada. De este modo podremos realizar un desarrollo similar al caso anterior. Finalmente, el tercer sumando de esta expresión contiene la integral: 1 2πiZc+i∞ c−i∞ ˆ f(s)ds (2.37) Por lo tanto necesitamos que esta integral sea convergente. Esta expresión es l´ımt→0+f(t) pues es la fórmula de inversión de la transformada de Laplace para t= 0 . Para que sea convergente necesitamos que exista el límite de la función f(t) en 0 . Si se cumpliesen todas las condiciones anteriores, podríamos continuar con el método. 2.7. NO RÍGIDAMENTE PRECISO 25 El primer término seria: 1 K ∞ X n=0 bTA−1 n X k=0 Wh n−k(ˆ f)gk! Wh n(ˆ f)≈R−n L L−1 X l=0 ˆ f∆(Reit) hCnldt siendo C=ei2π L El tercer término podemos calcularlo directamente: hf(0)K−1 KbT∞ X n=0 gnzn Para el segundo término, supongamos que se cumple que, utilizando las curvas denidas para (2.29): 1 2πiZc+i∞ c−i∞ sˆ f(s) λ−sds= l´ım n→∞−1 2πiZφn,2 sˆ f(s) s−λds+1 2πiZφn,1 sˆ f(s) s−λds = l´ım n→∞ 1 2πiZφn sˆ f(s) s−λds=λˆ f(λ) para todo λ∈C tal que su parte real sea mayor que c. Consideraremos que ∆(z) h es siempre diagonalizable. Se podría hacer los mismos cálculos, como en el apartado de rígidamente preciso, si la matriz no fuese diagonalizable pero los términos serián más extensos. Sea A∈Mmxm(C) con su diagonalización A=LADAL−1 A siendo DA la matriz diagonal asociada a A y LA la matriz de cambio de base. Denimos ˜ DA como: h˜ DAikj :=    [DA]kj ˆ f[DA]kj si j=k 0 si k6=j Denimos ˜ f:A∈Mm×m(C)→Mmxm(C) tal que ˜ f(A) := LA˜ DAL−1 A . Entonces podemos expresar el segundo sumando como: hK−1 K ∞ X n=0 bT n X k=0 ˜ Wh n−k(˜ f)gk! ˜ Wh n(ˆ f)≈R−n L L−1 X l=0 ˆ f∆(Reit) hCnldt siendo C=ei2π L donde ˜ Wh n−k son los coecientes de la descomposición en serie de potencias de ˜ f(∆(z) h) . 32 CAPÍTULO 3. RESULTADOS NUMÉRICOS Figura 3.2: Resultado utilizando métodos rígidamente precisos caso en el que x= 1 + √2 2 con el mismo método Runge-Kutta. 3.2.2. Polinomios Realizaremos ahora otro test académico en el que la solución del problema sea una función polinómica. Afrontaremos el problema: t2∗12t=Zt 0 12τ(t−τ)2dτ=t4 (3.2) En efecto, t4 es la solución de esta integral de convolución pues por el teorema de convolución tenemos que la Transformada de Laplace de una convolución es el producto de las Transformadas de sus factores y como la transformada de tn es Ltn(s) = n! sn+1 , podemos ver: L[t2∗12t] = L[t2]L[12t] = 2 s3 12 s2=4! s5=L[t4] Por lo tanto, tomando como g(t) = 12t y como ˆ f(s) = 2 s3 , podemos resolver la integral de convolución (3.2) con el método. Como conocemos el resultado exacto podremos realizar un estudio del error de aproximación del método. El las tablas 3.4-3.8 podemos ver los 3.2. TEST ACADÉMICO 33 errores de aproximación de los métodos rígidamente precisos para el problema. Puesto que es necesario tomar límites en el dominio de deción del problema, lo consideramos denido en [0,3] . Tabla 3.4: Error, Euler Implícito, ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 4,9500e+ 00 2,4524e+ 00 1,2206e+ 00 4,8689e−01 2,4320e−01 1,0055 Tabla 3.5: Error, Radau IIA de 2 etapas, ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 3,4754e−04 9,3393e−05 2,5393e−05 4,9489e−06 4,3653e−05 1,8474 Tabla 3.6: Error, Radau IIA de 3 etapas, ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 3,7454e−04 9,6768e−05 2,5815e−05 4,9759e−06 4,3688e−05 1,8767 Podemos ver que, dependiendo del método Runge-Kutta que utilicemos, el error del método de cuadratura de convolución varia. Para los métodos Radau IIA podemos ver que las aproximaciones con un paso de h= 3/1000 son más exactas que con un paso de h= 3/2000 . Una de las posibles razones de esto es que el error de aproximación del método es menor que el error de aproximación de la máquina al realizar las operaciones, pues en alguno de los cálculos internos del método los valores se podrían estar acercando al error de la máquina. En cuanto al orden de convergencia de los métodos podemos ver como como el método de Euler Implícito da un orden de aproximadamente 1, como sería de esperar por el orden del método Runge-Kutta. En general, para el resto de métodos el orden es mayor que 1 ya que estamos utilizando métodos Runge-Kutta con órdenes más altos. En las tablas 3.9-3.14 podemos ver los resultados de la resolución del mismo problema, 34 CAPÍTULO 3. RESULTADOS NUMÉRICOS Tabla 3.7: Error, Lobatto IIIC de 2 etapas, ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 3,2268e−02 8,0336e−03 2,0030e−03 3,1927e−04 1,1095e−04 1,9172 Tabla 3.8: Error, Lobatto IIIC de 3 etapas, ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 4,0233e−01 1,9958e−01 9,9397e−02 3,9664e−02 1,9791e−02 1,0051 utilizando métodos no rígidamente precisos. Podemos observar una disparidad entre los errores de aproximación de los diferentes métodos. Así como para algunos métodos, como el método de Crouzeix, tienen errores similares a los de los métodos rígidamente precisos. Por otro lado, algunos métodos, como el método diagonal implícito con x= 1 + √2/2 , resultan en unas aproximaciones con errores muy altos. Cada método será más o menos preciso en función del problema que vayamos a aproximar. En general, los métodos rígidamente precisos tendrán errores menores que los métodos no rígidamente precisos. Aunque puede que para casos particulares, algún método no rígidamente preciso pueda dar resultados más exactos que un método rígidamente preciso especíco. Para el resto de las aplicaciones que afrontaremos en este capítulo, nos centraremos en los métodos Runge-Kutta rígidamente precisos únicamente. De este modo, tendremos, en general, errores de aproximación menores. Además, no tendremos que preocuparnos sobre las condiciones necesarias que necesitamos para el uso del método de cuadratura de convolución en el caso no rígidamente preciso. 3.2. TEST ACADÉMICO 35 Tabla 3.9: Error, Pareschi y Russo de 2 etapas con x= 1 + √2/2 , ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 2,7494e−01 6,7700e−02 1,6795e−02 2,6739e−03 6,9336e−04 1,9993 Tabla 3.10: Error, Pareschi y Russo de 2 etapas con x= 0,3 , ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 8,7907e−03 2,2018e−03 5,5219e−04 8,9208e−05 3,6842e−05 1,8626 Figura 3.3: Resultados aplicando métodos no rígidamente precisos. 36 CAPÍTULO 3. RESULTADOS NUMÉRICOS Tabla 3.11: Error, Diagonalmente Implícito de 2 etapas con x= 1 + √2/2 , ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 8,5019e+ 00 4,1996e+ 00 2,0869e−00 8,3170e−01 4,1531e−01 1,0073 Tabla 3.12: Error, Diagonalmente Implícito de 2 etapas con x= 0,8 , ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 3,9496e+ 00 1,9594e+ 00 9,7584e−01 3,8941e−01 1,9453e−01 1,0047 Tabla 3.13: Error, Crouzeix de 2 etapas, ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 2,0001e−04 7,4952e−05 2,3088e−05 4,8014e−06 4,3684e−05 1,6306 Tabla 3.14: Error, Nørsett de 3 etapas, ec. polinómica h= 3/100 h= 3/200 h= 3/400 h= 3/1000 h= 3/2000 Orden Aprox. 3,7430e−04 9,6708e−05 2,5800e−05 4,9735e−06 4,3664e−05 1,8766 3.3. EDOS LINEALES 37 3.3. EDOs Lineales En esta sección afrontaremos la resolución numérica de ecuaciones diferenciales ordinarias lineales no homogéneas. Es decir, ecuaciones del tipo: xn)(t) + an−1xn−1)(t) + ···+a1x0(t) + a0x(t) = g(t) (3.3) La solución de este tipo de ecuaciones, en el caso en el que g(t) = 0 , se puede obtener analíticamente. Pero en el momento en el que se introduce una función g(t) no nula, la solución de la EDO se puede volver complicada y puede no ser posible obtener analíticamente. Buscaremos transformar la resolución de un problema de Cauchy en un problema que podamos resolver utilizando integrales de convolución. Comenzaremos desarrollando el método para EDOs de orden 1, resolviendo algún ejemplo particular, y nalmente desarrollaremos el método para EDOs de orden n. Esta aplicación del método de integrales de convolución viene inspirada por la misma aplicación que se realiza en el trabajo [3]. 3.3.1. EDOs de Primer Orden Afrontaremos la resolución numérica de EDOs lineales de primer orden no homogéneas. Particularizando en la expresión (3.3) para n= 1 e introduciendo una condición inicial, obtenemos la expresión general de los problemas de Cauchy con los que vamos a trabajar: (x0(t) + ax(t) = g(t) x(0) = x0 (3.4) Introducimos la condición inicial x(0) = 0 para obtener el problema de Cauchy (3.5). Con esta condición inicial estamos restringiendo este estudio a únicamente un problema de Cauchy particular. Esto no es ningún inconveniente, ya que, las diferentes soluciones de una EDO lineal no homogénea son de la forma: xp(t) + xh(t) donde xp(t) es una solución particular de la EDO no homogénea y xh(t) es una solución de la EDO homogénea, que son fáciles de calcular. De este modo tomando la solución de la EDO homogénea xh(t) tal que xh(0) = x0 y sumándosela a la solución obtenida para (3.5) obtendremos la solución de nuestro problema (3.4). (x0(t) + ax(t) = g(t) x(0) = 0 (3.5) Para resolver este problema particular, introducimos el siguiente problema dependiente de un parámetro k que consideraremos real positivo ( k∈[0,∞) ): (df dt(t;k) + af(t;k) = δ(t−k) f(0; k) = 0 (3.6) 38 CAPÍTULO 3. RESULTADOS NUMÉRICOS De este modo podemos ver que si denimos ˜x(t) como: ˜x(t) := Z∞ 0 f(t;k)g(k)dk (3.7) podemos comprobar como ˜x(t) es, en efecto, la solución del problema (3.5). ˜x0(t) + a˜x(t) = Z∞ 0df dt(t;k) + af(t;k)g(k)dk=Z∞ 0 δ(t−k)g(t)dk=g(t) ˜x(0) = Z∞ 0 f(0; k)g(k)dk= 0 Por tanto, por unicidad de solución, ˜x(t) es la única solución del problema de Cauchy (3.5). Una vez tenemos esto, solo es necesario calcular la solución del problema (3.6). Para esto debemos darnos cuenta que δ(t−k)=0 para todo k∈R+ , k6=t . Con lo que, para cada k∈R+ , en los intervalos (0, k) y (k, ∞) , f(t;k) es solución de la EDO homogénea. Conocemos las soluciones de la EDO homogénea, son de la forma ce−at con c∈R . Por tanto conocemos la expresión de la función f(t;k) : f(t;k) = (c1e−at si t∈(0, k) c2e−at si t∈(k, ∞) Puesto que tenemos como condición inicial del problema (3.6), f(0; k) = 0 , entonces podemos asegurar que c1= 0 . Por otro lado, para que se verique el problema (3.6), tiene que cumplirse que: l´ım t→k+f(t;k)−l´ım t→k− f(t;k)=1 Observación 3.2. En efecto, como la función δ(t−k) podemos entenderla como la derivada de la función heavyside H(t−k) denida como 1 si t≥k y 0 si t<k , entonces la función f(t;k) , jando k, tiene que tener un salto de longitud 1 en t=k , es decir, l´ımt→k+f(t;k)− l´ımt→k−f(t;k) = 1 . De este modo en todo punto tal que t6=k la función f(t;k) será la solución de la EDO homogénea y en el punto t=k tendremos que df dt(k;k) = δ(k−k) , vericando la ecuación de (3.6) para todos los puntos. con lo que podemos calcular el valor de c2 : c2e−ak = 1 c2=eak De este modo podemos expresar la función f(t;k) como: f(t;k) = F(t−k) := H(t−k)e−a(t−k) 3.3. EDOS LINEALES 39 donde H(t−k) es la función que toma el valor de 0 si t<k y toma el valor de 1 si t≥k . De este modo introduciendo esta expresión en la integral de (3.7) podemos obtener una fórmula con la que poder calcular la solución del ploblema de Cauchy (3.5). ˜x(t) = Z∞ 0 f(t;k)g(k)dk=Z∞ 0 F(t−k)g(k)dk =Zt 0 F(t−k)g(k)dk= [F∗g](t) Por lo que podemos calcular la solución buscada, calculando la convolución de las funciones F y g . La función F(t) = H(t)e−at podemos considerarla como si fuese F=e−at ya que todas las funciones con las que estamos trabajando durante el trabajo estamos suponiendo que están denidas únicamente sobre R+ . Por lo tanto, conocemos su Transformada de Laplace ˆ F(s) = 1 s+a puesto que es la transfromada de la función exponencial. De este modo hemos reducido el cálculo de la solución del problema de Cauchy (3.5) al cálculo de una integral de convolución. Veamos ahora los resultados de aplicar este método a la resolución de un ejemplo concreto. Consideremos la siguiente EDO lineal de primer orden no homogénea: (x0(t)+3x(t) = cos(t)−sin(t) x(0) = 3 (3.8) Calculando analíticamente la solución del problema podemos ver que e−3t es solución de la EDO homogénea y que (−2 sin(t) + 4 cos(t))/10 es una solución particular de la EDO no homogénea. Por lo tanto, las soluciones de la EDO no homogénea son del tipo (−2 sin(t) + 4 cos(t))/10 + Ce−3t con C∈R . Entonces particularizando para la condición inicial x(0) = 3 tenemos que x(t)=(−2 sin(t) + 4 cos(t) + 26e−3t)/10 es la solución exacta del problema de valor inicial (3.8). Para aproximar la solución de este problema, primero buscamos la solución de la EDO homogénea que cumple la condición inicial x(0) = 3 , esta es x(t)=3e−3t . Después calculamos la solución aproximada de la EDO no homogénea con condición inicial x(0) = 0 realizando la convolución explicada anteriormente y nalmente sumamos las dos funciones para obtener la solución de nuestro problema de Cauchy (3.8). Tabla 3.15: Error, Euler implícito, EDO lineal orden 1 h= 4/100 h= 4/200 h= 4/400 h= 4/1000 h= 4/2000 Orden Aprox. 9,8804e−03 5,0498e−03 2,5534e−03 1,0285e−03 2,4194e−03 0,9831 40 CAPÍTULO 3. RESULTADOS NUMÉRICOS Tabla 3.16: Error, Radau IIA de 2 etapas, EDO lineal orden 1 h= 4/100 h= 4/200 h= 4/400 h= 4/1000 h= 4/2000 Orden Aprox. 3,6503e−06 4,6305e−07 5,8337e−08 5,9370e−09 2,4351e−03 2,7976 Tabla 3.17: Error, Radau IIA de 3 etapas, EDO lineal orden 1 h= 4/100 h= 4/200 h= 4/400 h= 4/1000 h= 4/2000 Orden Aprox. 4,8580e−10 1,5331e−11 2,1276e−12 5,50943 −09 2,4346e−03 3,9175 En las tablas 3.15-3.19 podemos ver los errores de aproximación de la solución del problema (3.8), en el intervalo [0,4] , utilizando el método desarrollado en esta sección. Se muestran los resultados utilizando los cuatro métodos Runge-Kutta rígidamente precisos de la tabla 3.1 y Euler implicíto. Podemos observar como el error disminuye a medida que tomamos un parámetro de discretización menor hasta que llegamos a un punto donde los errores de aproximación aumentan de nuevo. Esto podría ser debido a que los errores de las operaciones realizadas por el ordenador podrían ser mayores que el error de aproximación del método, por lo que el error total aumenta en lugar de disminuir. En este caso, en mayor grado que en la sección anterior, podemos ver como el orden del método aumenta bastante a medida que lo hace el orden del integrador Runge-Kutta que utilicemos. Esto último se ve, sobre todo, en los métodos Radau IIA. 3.3.2. EDOs de orden n Afrontaremos ahora el problema de Cauchy asociado con la EDO lineal de orden n no homogénea (3.3). Introduciendo condiciones iniciales podemos obtener el problema de Cauchy que trataremos de resolver:                  xn)(t) + an−1xn−1)(t) + ···+a1x0(t) + a0x(t) = g(t) x(0) = x0 x0(0) = x1 . . . xn−1)(0) = xn−1 (3.9) Como en el caso de orden 1, resolveremos el problema de Cauchy para la EDO no homogénea con las condiciones inicales iguales a 0. Del mismo modo que en el caso anterior, podremos obtener fácilmente una solución de la EDO homogénea que cumpla las condiciones iniciales de (3.9) y sumársela a la solución del problema con condiciones iniciales nulas para obtener la solución del problema con condiciones iniciales no nulas. 3.3. EDOS LINEALES 41 Tabla 3.18: Error, Lobatto IIIC de 2 etapas, EDO lineal orden 1 h= 4/100 h= 4/200 h= 4/400 h= 4/1000 h= 4/2000 Orden Aprox. 4,0061e−04 1,0440e−04 2,6656e−05 4,3196e−06 2,4349e−03 1,9682 Tabla 3.19: Error, Lobatto IIIC de 3 etapas, EDO lineal orden 1 h= 4/100 h= 4/200 h= 4/400 h= 4/1000 h= 4/2000 Orden Aprox. 9,1414e−04 4,5556e−04 2,2742e−04 9,0873e−05 1,0264e−04 1,0025 Resolveremos entonces el problema de valor inicial:                  xn)(t) + an−1xn−1)(t) + ···+a1x0(t) + a0x(t) = g(t) x(0) = 0 x0(0) = 0 . . . xn−1)(0) = 0 (3.10) para ello introduciremos el problema dependiente del parámetro k∈R+ :                  dnf dtn(t;k) + an−1dn−1f dtn−1(t;k) + ···+a1df dt(t;k) + a0f(t;k) = δ(t−k) f(0; k)=0 df dt(0; k)=0 . . . dn−1f dtn−1(0; k)=0 (3.11) Análogamente al caso de orden 1, podemos denir la función ˜x(t) ˜x(t) := Z∞ 0 f(t;k)g(k)dk (3.12) de modo que será la única solución del problema de Cauchy (3.10). En efecto podemos ver como ˜xn)(t) + an−1˜xn−1)(t) + ···+a1˜x0(t) + a0˜x(t) = Z∞ 0dnf dtn(t;k) + an−1 dn−1f dtn−1(t;k) + ···+a1 df dt(t;k) + a0f(t;k)g(k)dk= Z∞ 0 δ(t−k)g(k)dk=g(t) y por tanto ˜x(t) es la solución de (3.10). Por tanto, debemos únicamente calcular la solución del problema (3.11). Para ello consideramos los intervalos, para cada k∈R+ , (0, k) y (k, ∞) , donde f(t;k) es una de las soluciones de la EDO homogénea. Puesto que el conjunto 48 CAPÍTULO 3. RESULTADOS NUMÉRICOS de este modo obtenemos una expresión análoga a la expresión (2.16) para cada etapa del método Runge-Kutta. Observación 3.6. En efecto, como Ai es la i-ésima la de la matriz A, AiA−1 es el vector la con un 1 en la i-ésima posición y ceros en las demás. La la Ai hace el papel que hace el vector bT en el desarrollo del método en el segundo capítulo. En efecto, si el método es rígidamente preciso se cumple que bT=Am y por tanto xn,m =xn . Utilizando la expresión (2.15) podemos llegar a: ∞ X n=0 xn,izn=AiA−1 ∞ X n=0 Xnzn!=AiA−1∆(z) h−sI−1∞ X n=0 gnzn Por tanto aproximando la expresión (2.5) en los tiempos tn+cih con i∈ {1,··· , m} . Llegamos a la expresión: yn, i =1 2πiZc+i∞ c−i∞ ˆ f(s)xn,ids de donde, realizando una transformada Z e introduciendo la expresión anterior llegamos a una expresión similar a la expresión (2.18): ∞ X n=0 yn,izn=∞ X n=0 1 2πiZc+i∞ c−i∞ ˆ f(s)AiA−1∆(z) h−sI−1 gnds!zn Podemos realizar los mismos procedimientos que en las secciones 2.4 y 2.5 con el vector la Ai haciendo el papel del vector bT . De modo que llegamos a la expresión nal: ∞ X n=0 yn,izn=AiA−1∞ X n=0 ˆ f∆(z) hgnzn De donde podemos obtener, agrupando cada componente del vector Yn e introduciendo la serie de potencias de la función ˆ f(∆(z) h) una expresión similar a la expresión (2.35): Yn= n X k=0 Wh n−k(ˆ f)gk Wh n(ˆ f)≈R−n L L−1 X l=0 ˆ f∆(Reit) hCnldt siendo C=ei2π L (3.18) De este modo hemos obtenido un sistema de ecuaciones determinado con el que podremos resolver la ecuación integral planteada en el problema. Se puede ver un desarrollo más riguroso en el artículo [11]. Observación 3.7. Para la resolución numérica del sistema determinado por la expresión (3.18), no será necesario resolver un gran sistema de ecuaciones, sino que se puede resolver un sistema mas pequeño en cada paso de tiempo, puesto que para calcular el los términos del vector gn no son necesarios los vectores posteriores, gn+1, gn+2 ··· . 3.4. ECUACIONES INTEGRALES 49 3.4.1. Ejemplo Nos centraremos ahora en un caso particular. Resolveremos la siguiente ecuación integral: exp(t)∗x(t) = Zt 0 et−τx(τ)dτ=sen(t) (3.19) Cuya incógnita x(t) es uno de los factores de la convolución. Podemos resolver este problema analíticamente, puesto que estas funciones particulares nos lo permiten. Las Transformadas de Laplace de las funciones seno y exponencial son conocidas, por lo tanto utilizando el Teorema de Convolución llegamos a: L[exp(t)∗x(t)] = L[exp(t)]L[x(t)] = 1 s−1L[x(t)] = L[sin(t)] = 1 s2+ 1 L[x(t)] = s−1 s2+ 1 =L[cos(t)] −L[sen(t)] Por lo tanto la solución exacta de la ecuación (3.19) es x(t) = cos(t)−sen(t) . Podemos utilizar el método desarrollado en esta sección para aproximar esta solución. Utilizaremos los métodos rígidamente precisos Radau IIA de 2 y 3 etapas, Lobatto IIIC de 2 etapas y el método de Euler Implícito. Tabla 3.25: Error, Euler Implícito, ec. integral h= 1/200 h= 1/400 h= 1/800 h= 1/1200 h= 1/1600 Orden Aprox. 5,0465e−03 2,5013e−03 1,2503e−03 8,3345e−04 6,2507e−04 1,0038 Tabla 3.26: Error, Radau IIA de 2 etapas, ec. integral h= 1/200 h= 1/400 h= 1/800 h= 1/1200 h= 1/1600 Orden Aprox. 4,9280e−03 2,5018e−03 1,2505e−03 8,3358e−04 6,2514e−04 0,9939 Tabla 3.27: Error, Radau IIA de 3 etapas, ec. integral h= 1/200 h= 1/400 h= 1/800 h= 1/1200 h= 1/1600 Orden Aprox. 5,1350e−03 2,5039e−03 1,2508e−03 8,3371e−04 6,2521e−04 1,0111 En las tablas 3.25-3.28 podemos ver los errores de aproximación de cada uno de los métodos utilizados para resolver la ecuación (3.19). Podemos observar como, independientemente del método, el orden de la aproximación se encuentra en torno a 1. Puesto que 50 CAPÍTULO 3. RESULTADOS NUMÉRICOS Tabla 3.28: Error, Lobatto IIIC de 2 etapas, ec. integral h= 1/200 h= 1/400 h= 1/800 h= 1/1200 h= 1/1600 Orden Aprox. 5,0041e−03 2,5010e−03 1,2503e−03 8,3345e−04 6,2507e−04 1,0003 no hemos realizado el estudio de convergencia del método no podemos armar nada, pero en vista de los resultados podríamos conjeturar que el método tiene orden 1 pues el orden de las aproximaciones en las etapas intermedias de los métodos que estamos utilizando es de orden 1. Por esta razón para obtener mayor orden quizá necesitásemos métodos que tuvieran ordenes superiores para las etapas intermedias. Observación 3.8. En efecto, si alguna de las componentes de los vectores gn tuviese orden de aproximación de 1, como utilizamos los vectores anteriores para calcular los siguientes, el orden de las aproximaciones siguientes se vería reducido a orden 1. Códigos En este apéndice veremos los códigos de Matlab que fueron utilizados para la resolución de los problemas que se presentaron a lo largo del trabajo. Estas dos primeras funciones permiten calcular la aproximación de una integral de convolución con el método desarrollado en este trabajo. El programa funciona para métodos Runge-Kutta rígidamente precisos. Toma como argumentos de entrada los datos del integrador Runge-Kutta A,b,c , la función ˆ f(s) Transformada de Laplace de la función f(t) , la función g(t) , la longitud del intervalo en el que calculamos la integral de convolución y los parámetros de discretización npas, L y R. El parámetro npas es el número de pasos del método, que es equivalente a dar el parámetro h. La salida del programa será un vector y que contendrá las aproximaciones de la función y(t) en los nodos de discretización. function [ y]=convolucion_nucleo (A, b, c , f , g , npas ,T,L,R) h=T/npas ; t =[0 ,(1: npas) * h ] ; y= zeros (1 , npas+1); dim= length (b ); A_1= inv (A); gn= zeros (dim , npas ); for i =1:npas gn (: , i)=g( t ( i)+h * c ); end [W]=calc_nucleo (A_1, f , npas = 1,h ,L,R); for n=0:npas = 1 for k=0:n y(n+2)=y(n+2)+W(n = k+1 ,:) * gn (: , k+1); end end end 51 52 CÓDIGOS function [W]=calc_nucleo (A_1, f ,N,h ,L,R) dim= length (A_1); W= zeros (N+1,dim ); C= exp (2 * pi * 1 i /L); e= zeros (1 ,dim ); e(dim)=1; for l =0:L = 1 delta=A_1 * ( eye (dim) = R * C^( = l ) * sparse (1: dim , . . . dim * ones (1 ,dim) , ones (1 ,dim) ,dim ,dim )); delta=delta./h; [V,D]= eig ( delta ); for j =1:dim D( j , j)=f (D( j , j )); end V_1= inv (V); for n=0:N W(n+1,:)=W(n+1,:)+e * V * D * V_1 * C^(n * l ); end end for n=0:N W(n+1 ,:)=((R^( = n))/L). * W(n+1 ,:); end end Para el caso no rígidamente preciso, podemos utilizar las siguientes 3 funciones. En contrapartida con el caso anterior, para este es necesario introducir el valor de l´ımx→0+= f0 y en lugar de introducir R, utilizaremos epsilon para calcular R internamente en el programa, ya que en el caso no rígidamente preciso debemos variar el parámetro R. function [ y]=convolucion_NRP(A,b , c , f , g , f_0 , npas ,T,L, epsilon ) h=T/npas ; t =[0 ,(1: npas) * h ] ; A_1= inv (A); unos=ones ( length (b) ,1); K=b' * A_1 * unos ; [y_1]=convolucion_primertermino (A,b , c , f ,g , npas ,T,L, epsilon ); [y_2]=convolucion_segundotermino (A,b, c , f , g , npas ,T,L, epsilon ); gn= zeros ( length (b) , npas ); 53 for i =1:npas gn (: , i)=g( t ( i)+h * c ); end y_3=[0 ,b' * gn ] ; y=K^( = 1) * y_1+h * (1 = K^( = 1)) * (y_2+f_0 * y_3); end function [W]=calc_nucleo_primertermino (b ,A_1, f ,N,h,L, epsilon ) dim= length (A_1); W= zeros (N+1,dim ); C= exp (2 * pi * 1 i /L); unos=ones ( length (b ) ,1); K=b' * A_1 * unos ; e=b' * A_1; R= min ([1 , abs (1/(1 = K))]) = epsilon ; for l =0:L = 1 delta=A_1 = R * C^( = l )/(1+R * C^( = l ) * (K = 1)) * A_1 * unos * e ; delta=delta./h; [V,D]= eig ( delta ); for j =1:dim D( j , j)=f (D( j , j )); end for n=0:N W(n+1,:)=W(n+1,:)+e * V * D * inv (V) * C^(n * l ); end end for n=0:N W(n+1 ,:)=((R^( = n))/L). * W(n+1 ,:); end end function [W]=calc_nucleo_segundotermino (b ,A_1, f ,N,h ,L, epsilon ) dim= length (A_1); W= zeros (N+1,dim ); C= exp (2 * pi * 1 i /L); unos=ones ( length (b ) ,1); K=b' * A_1 * unos ; 54 CÓDIGOS R= min ([1 , abs (1/(1 = K))]) = epsilon ; for l =0:L = 1 delta=A_1 = R * C^( = l )/(1+R * C^( = l ) * (K = 1)) * A_1 * unos * b' * A_1; delta=delta./h; [V,D]= eig ( delta ); %delta * V=V * D for j =1:dim D( j , j)=D( j , j ) * f (D( j , j ) ); end for n=0:N W(n+1,:)=W(n+1,:)+b' * V * D * inv (V) * C^(n * l ); end end for n=0:N W(n+1 ,:)=((R^( = n))/L). * W(n+1 ,:); end end En el caso de la resolución de ecuaciones integrales es necesario utilizar otros programas diferentes. En estos programas la función y(t) es un dato mientras que la función g(t) será la salida del programa. function [ g]=ec_int (A, c , f ,y , npas ,T,L,R) h=T/npas ; t =(0:npas ) * h; dim= length (c ); A_1= inv (A); [W]=calc_nucleo_ecint (A_1, f , npas = 1,h ,L,R); W= real (W); g= zeros (1 , npas+1); G= zeros (dim , npas ); W_0=W(: , : ,1 ); for n=0:npas = 1 Y=y( t (n+1)+c * h ); Wg= zeros (dim , 1 ); for k=0:n = 1 Wg=Wg+W( : ,: , n = k+1) * G(: , k+1); end g_n1=W_0\(Y = Wg); 55 G(: ,n+1)=g_n1; end g (2: end )=G(dim , : ) ; g(1)=g (2); end function [W]=calc_nucleo_ecint (A_1, f ,N,h ,L,R) dim= length (A_1); W= zeros (dim ,dim ,N+1); C= exp (2 * pi * 1 i /L); for l =0:L = 1 delta=A_1 * ( eye (dim) = R * C^( = l ) * sparse (1: dim , dim * ones (1 ,dim) , ones (1 ,dim) ,dim ,dim )); delta=delta./h; [V,D]= eig ( delta ); %delta * V=V * D for j =1:dim D( j , j)=f (D( j , j )); end V_1= inv (V); for n=0:N W(: , : , n+1)=W( : ,: , n+1)+V * D * V_1 * C^(n * l ); end end for n=0:N W(: , : , n+1)=((R^( = n))/L). * W(: , : , n+1); end end 56 CÓDIGOS Bibliografía [1] C.Lubich, Convolution Quadrature and Discretized Operational Calculus I. , Numer. Math. 52,129-145 (1988) [2] C.Lubich, Convolution Quadrature and Discretized Operational Calculus II. , Numer. Math. 52, 413-425 (1988) [3] A. Ríos, Aproximación numérica de integrales de convolución mediante el método de cuadratura de convolución Santiago de Compostela (Trabajo n de grado), Universidad de Santiago de Compostela [4] Lehel Banjai, Matthias Messner, Martin Schanz, Runge-Kutta convolution quadrature for the Boundary Element Method , Comput. Methods Appl. Mech. Engrg. 245-246 (2012) 90-101 [5] J.C Butcher, Numerical Methods for Ordinary Dierential Equations , 2nd ed., The University of Auckland, New Zealand [6] Hairer, Nørsett, Wanner, Solving Ordinary Dierential Equations I Nonsti Problems , (1993), Springer-Verlag Berlin Heidelberg. [7] Dyke, P.P.G., An Introduction to Laplace transforms and Fourier series ,London [etc.] : Springer, cop. 2000 [8] I. Márquez, J.J. Nieto, Variable Compleja , Santiago de Compostela (2017), Universidad de Santiago de Compostela [9] Lloyd N. Trefethen and J. A. C. Weideman, The Exponentially Convergent Trapezoidal Rule , SIAM Rev., 56(3), 385458. (74 pages) [10] Wikipedia,(5/06/2020), List of RungeKutta methods , recuperado de https://en. wikipedia.org/wiki/List_of_Runge-Kutta_methods 57