36
Prelininares.
1.4.
El problema de Cauchy. Resoluci´ on num´ erica de E.D.: generalidades.
1.4.1.
El problema de Cauchy. Existencia y unicidad.
Aqu´ı, enumeramos los resultados principales de existencia y unicidad sde soluci´on del problema de valores iniciales (p.v.i.) siguiente: y 0 = f (x, y) y (x0 ) = µ
(1.29)
donde f : D ⊆ R2 → R e y es una funci´on desconocida de x. Se dice que y(x) es una soluci´on de (1.29) en un intervalo I tal que x0 ∈ I si y(x) es funci´on derivable en I, y 0 (x) = f (x, y (x)) ∀x ∈ I e y (x0 ) = µ. Ahora bien, ¿est´a asegurada la existencia de soluci´on?, ¿y la unicidad? La respuesta a sendas cuestiones se explicita en el teorema siguiente:
Teorema 1.36 (existencia y unicidad) Si f : D ⊆ R2 → R es continua en D = {(x, y) : a ≤ x ≤ b, y ∈ R} y satisface la condici´on de Lipschitz: |f (x, y1 ) − f (x, y2 )| < L |y1 − y2 |
(1.30)
(L es la cte de Lipschitz respecto de y). Entonces, el problema (1.29) admite una u ´nica soluci´on y(x) en [a, b] : x0 ∈ [a, b]. Ahora bien, en ocasiones, cuando una funci´on f es suficientemente suave la condici´on (1.30) puede sobre D. As´ı, si ∂f es continua y acotada en D, entonces el teorema evaluarse desde el an´alisis de ∂f ∂y ∂y sigue siendo v´alido. EJERCICIO 1 Busque en la bibliograf´ıa de E.D.O. ejemplos de p.v.i. en los que falle la existencia y/o unicidad de soluci´on justificando el motivo. EJERCICIO 2 Pruebe que el p.v.i. y 0 =
p |y|, y (0) = 0 admite soluci´on sobre R. ¿Es u ´nica?
EJERCICIO 3 Demuestre que el p.v.i. x0 = t2 + ex , x (0) = 0 admite soluci´on u ´nica en el intervalo |t| ≤ 0,351
Apuntes de Jer´onimo Lorente
37
Teorema 1.37 (existencia y unicidad) Si f : R ⊆ R2 → R es continua en el rect´angulo R = {(x, y) / |x − x0 | ≤ a, |y − µ| ≤ b} y |f (x, y1 ) − f (x, y2 )| < L |y1 − y2 | Entonces, el problema (1.29) admite una u ´nica soluci´on y(x) en [x0 − h, x0 + h] donde h = m´ın a, Mb y M = cota de |f (x, y)| sobre R. Este teorema se puede extender, f´acilmente, para el caso de sistemas de ecuaciones diferenciales ordinarias de primer orden sin m´as que modificar en (1.30) los valores absolutos por normas vectoriales/matriciales seg´ un el caso.
1.4.2.
Generalidades sobre m´ etodos de resoluci´ on num´ erica de E.D.
A lo largo de los cap´ıtulos siguientes, analizaremos distintos m´etodos para resolver num´ericamente tanto el p.v.i. (1.29) como otros que precisaremos en su momento. Dado que se utilizan con frecuencia ecuaciones en diferencias finitas, damos aqu´ı los aspectos m´as b´asicos de este tipo de ecuaciones. Ecuaciones en diferencias. Sea V un espacio vectorial de sucesiones {yn } sobre R. Este es un espacio de dimensi´on infinita. Sea L : V → V un operador lineal y sea el operador desplazamiento E:V →V y 7→ Ey/ (Ey)n = yn+1 n = 1, 2, . . . Entonces, todo operador sobre V dado como combinaci´on lineal de potencias del operador desplazamiento se llama operador en diferencias lineales finitas; es decir, L=
m X
ci E i
(1.31)
i=0
donde, el operador E i es tal que: (E i y)n = yn+i i = 0, 1, . . .; ci ∈ R. m P i El subconjunto H = L : V → V : L = L = ci E es un subespacio propio de V con dimensi´on m+1 cuya base es {I, E, . . . , E m }.
i=0
El operador L dado en (1.31) tiene asociado un polinomio de grado m: p (λ) =
m X
ci λi
i=0
llamado polinomio caracter´ıstico asociado a L. As´ı, L = p (E).
38
Prelininares. Diremos que una ecuaci´on es en diferencias finitas (E.D.F.), homog´enea, si es de la forma: Ly = 0
(1.32)
donde y es la inc´ognita de la ecuaci´on. Resolver (1.32) ser´a encontrar todas las sucesiones y = {yn } de V que satisfacen la ecuaci´on. Dado que L es un operador lineal sobre V resolver (1.32) equivale a conocer el subespacio ker (L). Ahora bien, en la pr´actica, (1.32) puede aparecer en muy diversas formas; por ejemplo: Analice la ecuaci´on en diferencias finitas: a) yn+2 − 3yn+1 + 2yn = 0, n = 1, 2, . . . ¿{yn }?; ´o, b) (E 2 − 3E + 2E 0 ) y = 0 ¿y?; ´o, c) p (E) y = 0 con p (λ) = λ2 − 3λ + 2; ¿y? ¿C´omo calcular las soluciones de a), b) ´o c)? Para ello, ensayamos con tipos especiales de candidatas a soluciones (de modo similar al tratamiento de E.D.O. lineales de orden superior y coeficientes ctes). Por ejemplo, supongamos que tomamos y = {yn = rn }, entonces se ha de cumplir, en a), ( rn+2 − 3rn+1 + 2rn = 0 ∀n ⇔ (r2 − 3r + 2) rn = 0 ⇔
r = 0 (soluci´on trivial) r2 − 3r + 2 = 0 ⇔ r = 1 o´ 2
esto conduce a dos soluciones linealmente independientes: y1 = {1, 1, 1, . . . 1, . . .} e y2 = {2, 4, 8, . . . , 2n , . . .} y desde ´estas podemos generar otras mediante combinaciones lineales: y = α 1 y1 + α 2 y2 ¿Se puede expresar toda soluci´on de a) en la forma anterior? La respuesta la damos un poco m´as adelante. Del ejemplo se extrae la idea de que para resolver la E.D.F. (1.32) hemos de buscar los ceros del polinomio caracter´ıstico asociado a la ecuaci´on. Todo esto se aclara en el teorema siguiente: Teorema 1.38 1. Si p (λ) es el polinomio caracter´ıstico y r es un cero suyo, entonces y = {rn } es una soluci´on de p (E) y = 0. 2. Si todos los ceros de p (λ) son simples y no nulos (p.e. r1 , r2 , . . . , rm ), entonces toda soluci´on de p (E) y = 0 es combinaci´on lineal de y1 , y2 , . . . , ym donde yj = rjn j = 1, . . . , m; es decir, n y = α1 r1n + α2 r2n + · · · + αm rm
(1.33)
Apuntes de Jer´onimo Lorente
39
3. Si p (λ) tiene un cero, r, de multiplicidad k > 1 y p (0) 6= 0, entonces: ∂ j−1 y yj = j−1 j = 1, . . . , k ∂r
(1.34)
donde y = {rn } son soluciones linealmente independientes de p (E) y = 0. ¿Es una condici´on severa la hip´otesis de que p (0) 6= 0? ESTABILIDAD DE UNA E.D.F. Se dice que una E.D.F. p (E) y = 0 es estable si toda soluci´on est´a acotada; es decir, ∃C ∈ R tal que |yn | ≤ C ∀n. Teorema 1.39 (estabilidad) Sea p (λ) con p (0) 6= 0 el polinomio caracter´ıstico de (1.31), entonces son equivalentes: 1. La ecuaci´on (1.32) es estable; 2. Si r es un cero de p (λ), entonces: |r| < 1; ´o, |r| = 1 y r es simple EJERCICIO. Analice la estabilidad de la ecuaci´on en diferencias: 4yn + 7yn−1 + 2yn−2 − yn−3 = 0. Pasamos, ahora, a describir algunos aspectos generales sobre los m´etodos de resoluci´on num´erica de E.D.O. Un m´etodo de resoluci´on num´erica de un P.V.I. trata, fundamentalmente, 1. tomar una partici´on (habitualmente uniforme) {xn }N a definida n=0 , del intervalo [a,b] donde est´ la soluci´on, y (x), del problema (1.29). on. 2. obtener una sucesi´on o conjunto de valores {yn }N n=0 asociados a los nodos de la partici´ 3. objetivo: tal sucesi´on reproduzca la soluci´on exacta en cada punto del intervalo; es decir, que yn ' y (xn ) es una aproximaci´on tanto mejor cuanto mayor sea N. Pero, ¿c´omo se obtienen los valores yn ? Son varias las t´ecnicas (integraci´on y derivaci´on num´ericas, desarrollo de Taylor, interpolaci´on, etc..) usadas para el c´alculo de tales valores que conducen a m´etodos que se suelen clasificar en dos grupos: ´ METODOS DE UN PASO.Estos son los que consiguen el valor yn+1 a partir de los valores yn y xn (p.e. m´etodos de Euler, Runge-Kutta, etc...)
40
Prelininares.
´ METODOS MULTIPASO o ´ k-PASOS (k > 1).Estos son los que consiguen el valor yn+k a partir de los valores yn+k−1 , . . . , yn y xn+k−1 , . . . , xn (p.e. m´etodos de Adams-Bashfort, Milne-Simpson, etc.....) Una segunda clasificaci´on de m´etodos, en cada grupo, distingue entre: EXPL´ICITOS.- son aquellos en los que el valor yn+k se obtiene a parir de la evaluaci´on de una determinada funci´on que depende u ´nicamente de xn+k−1 , . . . , xn e yn+k−1 , . . . , yn ; es decir, yn+k = g (yn+k−1 , . . . , yn ; xn+k−1 , . . . , xn ; h, f ) (h =
b−a N
es el paso del m´etodo cuando la partici´on es uniforme)
IMPL´ICITOS.- son aquellos en los que el valor yn+k se obtiene a partir de la resoluci´on de una ecuaci´on del tipo: G (yn+k , yn+k−1 , . . . , yn ; xn+k , xn+k−1 , . . . , xn ; h, f ) = 0 Con esto se puede decir que un m´etodo para obtener una soluci´on num´erica de (1.29) tiene la forma (en la mayor´ıa de los casos) siguiente: yn+k = αk−1 yn+k−1 + · · · + α0 yn + hΦ (xn ; yn , . . . , yn+k−1 , yn+k ; h, f )
(1.35)
donde supondremos que la funci´on Φ cumple ∀h < h0 :
si f ≡ 0 ⇒ Φ ≡ 0 ∗ |Φ x; yn∗ , . . . , yn+k ; h, f − Φ (x; yn , . . . , yn+k ; h, f ) | 6 M
(1.36) k X
∗ kyn+j − yn+j k
(1.37)
j=0
Definici´ on 1.12 Diremos que (1.35) es convergente si se satisfacen las dos condiciones siguientes: l´ım
h→0 x=a+nh
yn = y (x) para x ∈ [a, b];
l´ım yr (h) = µ r = 0, 1, . . . , k − 1
h→0
Una definici´on alternativa de la convergencia es: l´ım m´ax |yn − y (xn )| = 0
n→∞ 0≤n≤N
Definici´ on 1.13 Diremos que (1.35) es consistente con el p.v.i. si se satisface la condici´on siguiente: l´ım
h→0 x=a+nh
Rn+k =0 h
(1.38)
Apuntes de Jer´onimo Lorente
41
donde Rn+k = y (xn+k ) −
k−1 X
αj y (xn+j ) − hΦ (xn ; y (xn ) , . . . , y (xn+k−1 ) , y (xn+k ) ; h)
j=0
se conoce como error de truncatura local. Teorema 1.40 k−1 P El m´etodo es consistente con el p.v.i. sii p (λ) = λk − αj λj (llamado primer polinomio caracter´ıstij=0
co) cumple: p (1) = 0 Φ (xn ; y (xn ) , . . . , y (xn ) , y (xn ) ; 0) = p0 (1) f (xn , y (xn )) Definici´ on 1.14 Sean {δn } y {δn∗ } sendas perturbaciones en el m´etodo num´erico y {zn } y {zn∗ } las soluciones num´ericas respectivas de (1.39) y (1.40). Entonces el m´etodo es cero-estable sii existen ctes M, h0 tales que ∀h ∈ ]0, h0 ] se cumple: |zn − zn∗ | ≤ M ε 0 ≤ n ≤ N siempre que |δn − δn∗ | ≤ ε ∀n. zn+k = αk−1 zn+k−1 + · · · + α0 zn + h [Φ (xn ; zn , . . . , zn+k−1 , zn+k ; h, f ) + δn+k ] zj = µr (h) + δr
j = 0, . . . , k − 1
∗ ∗ ∗ ∗ ∗ zn+k = αk−1 zn+k−1 + · · · + α0 zn∗ + h Φ xn ; zn∗ , . . . , zn+k−1 , zn+k ; h, f + δn+k zj∗ = µr (h) + δj∗
j = 0, . . . , k − 1
(1.39)
(1.40)
Teorema 1.41 El m´etodo (1.35) es estable sii el primer polinomio caracter´ıstico cumple: todos sus ceros est´an en el disco unidad (|r| ≤ 1∀r/p (r) = 0); si un cero, r, est´a en la frontera del disco unidad (|r| = 1), entonces es simple.
Teorema 1.42 El m´etodo (1.35) es convergente si, y s´olo si, es consistente y cero-estable.
Por u ´ltimo, dos m´etodos para resolver un mismo problema pueden compararse respecto a su precisi´on usando el concepto de orden.
42
Prelininares.
Definici´ on 1.15 Diremos que (1.35) es de orden p > 1 si se satisface la condici´on siguiente: l´ım
h→0 x=a+nh
Rn+k = C (cte. respecto de h) hp+1
donde Rn+k es el error de truncatura local del m´etodo. En la pr´actica, para la determinaci´on del orden de un m´etodo, se intentar´a buscar el valor de p para el que se tiene: Rn+k = y (xn+k ) −
k−1 X
αj y (xn+j ) − hΦ (xn ; y (xn ) , . . . , y (xn+k ) ; h) = O hp+1
j=0
La t´ecnica, habitual, para este prop´osito es el uso de desarrollos de Taylor apropiados.