MÉTODOS MATEMÁTICOS (Curso 2011-2012) Tercer Curso de Ingeniería Aeronáutica Departamento de Matemática Aplicada II. Universidad de Sevilla

LECCIÓN 8: OTROS MÉTODOS NUMÉRICOS PARA ECUACIONES DIFERENCIALES ORDINARIAS. PROBLEMAS DE CONTORNO. Índice 1. Otros métodos: Métodos lineales multipaso 1.1. Métodos de Adams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2. Fórmulas de diferenciación regresiva (BDF) . . . . . . . . . . . . . . . . . . . . . . 1.3. Otros métodos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1 3 7

2. Problemas stiff

8

3. Otras cuestiones de orden práctico

12

4. El método del disparo para problemas de contorno 4.1. Cuestiones de repaso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. El método del disparo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13 13 16

5. Cuestiones y problemas

17

1.

Otros métodos: Métodos lineales multipaso

De los llamados métodos lineales multipaso, sólamente haremos unos breves comentarios, pues por un lado, el desarrollo en las últimas décadas de los pares encajados de métodos Runge-Kutta han hecho que este tipo de métodos pierdan mucha competitividad y, por otro, un estudio más detallado aporta poco en el orden práctico a lo que se haya podido aprender en la secciones anteriores. A su vez, estos métodos están muy bien documentados en multitud de libros de texto. Debe saber sin embargo que para algunos problemas muy determinados y concretos, presentan ventajas que, a día de hoy, ningún método Runge-Kutta ha conseguido igualar. Lo particular de estas situaciones no justifica, en todo caso, que vayamos más alla de la descripción somera que haremos en esta sección.

1.1.

Métodos de Adams

Son anteriores a los métodos Runge-Kutta y se deben a los británicos Adams, Bashforth (1883, 1885) y Moulton (1926). Recordemos que dados el PVI  y ′ (t) = f (t, y(t)), , t ∈ [t0 , tf ]. (1) y(t0) = y0 , 1

y una partición t0 < t1 < . . . < tN = tf , del intervalo [t0 , tf ], la solución y satisface la igualdad Z Z tn+1 ′ y (t) dt = y(tn+1 ) = y(tn ) +

tn+1

(2)

f (t, y(t)) dt.

tn

tn

Los métodos Runge-Kutta aproximan esta última integral mediante una fórmula de cuadratura con nodos en el intervalo [t0 , tf ] (y aproximan los valores del integrando con las etapas ki ). Los métodos de Adams aproximan dicha integral con una fórmula de cuadratura basada en los nodos de cuadratura tn−k+1, . . . , tn

o tn−k+1 , . . . , tn , tn+1 .

donde k es el número de pasos del método. Para ello, obviamente se requiere haber obtenido (con ∗ otro método) esas k aproximaciones anteriores. Si Pn,k el polinomio de grado k − 1 que interpola a los k valores del campo f (t, y) en las aproximaciones numéricas (tn , yn ), . . . , (tn−k+1 , yn−k+1), esto es,  j−1 k−1 Y X ∗ Pn,k (t) = (t − tn−l ) f [tn , . . . , tn−j ], j=0

l=0

donde f [tn , . . . , tn−j ] son las diferencias divididas definidas recurrentemente como f [tn−j ] = f (tn−j , yn−j ), j = 0, . . . , k − 1, f [tn , . . . , tn−j+1 ] − f [tn−1 , . . . , tn−j ] f [tn , . . . , tn−j ] = , tn − tn−j

j = 1, 2, . . . .

el correspondiente método de k-pasos es entonces Z tn+1 ∗ Pn,k (t) dt. yn+1 = yn +

(3)

tn

y se conoce como método de Adams explícito de k pasos (más conocido también como método Adams-Bashforth de k pasos, aunque no se deba a Bashforth). También podemos considerar el polinomio de grado k que interpola k valores del campo f (t, y) en las aproximaciones numéricas (tn , yn ), . . . , (tn−k+1, yn−k+1) anteriores, junto con el el valor de f en (tn+1 , yn+1 ), aunque no lo conozcamos, Pn,k (t) =

∗ Pn,k

 Y k−1 (t − tn−l ) f [tn+1 , tn , . . . , tn−k+1 ]. + l=0

El correspondiente método numérico será entonces Z tn+1 ∗ Pn,k (t) dt yn+1 − hn βn,k f (tn+1 , yn+1 ) = yn + tn

+f

(0)

[tn+1 , tn , . . . , tn−k+1]

Z

tn+1 k−1 Y

tn

2

l=0

(t − tn−l ) dt,

(4)

donde βn,k

1 = hn

Z

(0)

tn+1 k−1 Y

tn

l=0

t − tn−l dt, tn+1 − tn−l

y f [tn+1 , tn , . . . , tn−k+1] es la diferencia dividida que se obtiene cuando se cambia yn+1 por 0. Observe que en (4), la incógnita yn+1 aparece a la izquierda del signo igual, y que lo que hay a la derecha sólamente depende de los valores conocidos yn , . . . , yn−k+1. El correspondiente método se conoce como método de Adams implícito de k pasos, o también de Adams-Moulton de k pasos. Observe que para obtener el valor de yn+1 hay que resolver la ecuación (o sistema de ecuaciones) dado por la igualdad en (4). En la práctica esto se hace por iteración de punto fijo, procedimiento que convergerá siempre que hn sea suficientemente pequeño. Debe saber lo siguiente. El método de Adams explícito de un paso es el método de Euler yn+1 − yn = hn f (tn , yn ). El método de Adams implícito de un paso se conoce como la regla de los trapecios (iqual que la fórmula de cuadratura) y es de la fomra yn+1 − yn =


hn f (tn , yn ) + f (tn+1 , yn+1) . 2

El método de Adams explícito de k pasos es convergente de orden k, y el implícito de k pasos de orden k + 1 (para recordarlo, puede Vd. pensar en los casos k = 1). En la práctica se suelen emplear ambos métodos juntos, en lo que se conoce como pares ∗ predictor-corrector. Se utiliza el método explícito para obtener un valor de yn+1 que se toma como iterante inicial en la iteración de punto fijo necesaria para encontrar el el valor de yn+1 del método implícito. Existen procedimientos en los pares predictor-corrector para estimar el error local, con lo que generalmente se emplean con paso variable. De hecho, se pueden estimar los errores locales correspondientes a los métodos de k − 1 pasos y de k + 1 pasos, por lo que los métodos de Adams se emplean en algortimos de paso y orden variables. El comando de Matlab que utiliza los métodos de Adams con paso y orden variable es ode113. Empezando con k = 1, puede llegar a utilizar fórmulas de hasta k = 12 pasos.

1.2.

Fórmulas de diferenciación regresiva (BDF)

Las fórmulas de diferenciación regresiva (en inglés “backward differentiation formulae” (BDF)), son los métodos multipaso que se utilizan para problemas “stiff” que explicaremos en la sección siguiente. Datan de 1952, y fueron inventadas por Curtiss y Hirschfelder. Se definen como sigue. Supongamos conocidas las k aproximaciones yn−k+1, . . . , yn . Para obtener la aproximación yn+1 con la fórmula BDF de k pasos, se impone la condición de que el polinomio Qn,k 3

de grado k que interpola a los k + 1 pares (tn−k+1 , yn−k+1), . . . , (tn+1 , yn+1 ) (incluímos el desconocido yn+1 ) satisfaga la ecuación diferencial en tn+1 , esto es Q′n,k (tn+1 ) = f (tn+1 , Qn,k (tn+1 )) = f (tn+1 , yn+1). El polinomio Qn,k , escrito en forma de Newton es  j−1 k Y X (t − tn+1−l ) y[tn+1, . . . , tn+1−j ], Qn,k (t) = j=0

l=0

donde y[tn+1, . . . , tn+1−j ] es la diferencia dividida basada en las aproximaciones yn+1, . . . , yn+1−j , esto es y[tn−j ] = yn−j , j = −1, 0, 1, . . . , k − 1, y[tn+1 , . . . , tn−j+2] − y[tn , . . . , tn−j+1 ] y[tn+1, . . . , tn−j+1] = , tn+1 − tn−j+1

j = 1, 2, . . . .

Derivando Qn,k , evaluando en tn+1 , igualando a f (tn+1 , yn+1) y multiplicando por hn = tn+1 − tn obtenemos la expresión de la fórmula BDF de k pasos hn

j−1 k Y X j=1

l=1

 (tn+1 − tn+1−l ) y[tn+1 , . . . , tn+1−j ] = hn f (tn+1 , yn+1 ).

(5)

Si dejamos a un lado del signo igual la incógnita yn+1 y al otro el resto, tenemos que el método se escribe como αn,k yn+1 − hn f (tn+1 , yn+1 ) = gn,k (yn , . . . , yn−k+1), (6) donde, obviamente, αn,k =

k X j=1

y gn,k (yn , . . . , yn−1) = hn

j−1 k Y X j=1

l=1

hn tn+1 − tn+1−j

 (tn+1 − tn+1−l ) y (0) [tn+1 , . . . , tn+1−j ],

donde, igual que en el apartado anterior, y (0) [tn+1 , . . . , tn+1−j ] denota la diferencia dividida que se obtiene al cambiar yn+1 por 0. Para encontrar yn+1 se resuelve la ecuación en (6) siempre por el método de Newton. Esto es [0] partiendo de una aproximación incial yn+1 (normalmente obtenida por extrapolación de los valores [µ] yn , . . . , yn−k+1), se obtienen aproximaciones sucesivas yn+1 , [µ+1]

[µ]

yn+1 = yn+1 + x[µ] ,

µ = 0, 1, 2, . . . ,

donde x[µ] es la solución del sistema lineal

[µ] [µ] αn,k I − hn fy x[µ] = − αn,k yn+1 − hn f (tn+1 , yn+1) − gn,k (yn , . . . , yn−k+1) , 4

(7)

siendo fy la matriz jacobiana 

[µ]

  fy =  

∂f1 ∂f1 ... ∂y1 ∂ym .. .. . . ∂fm ∂fm ... ∂y1 ∂ym



  , 

evaluada en (tn+1 , yn+1). Por tanto, para poder utilizar el software que emplea este método, no solo hay que proporcionar el segundo miembro f , sino tambien la matriz diferencial fy . Debe saber lo siguiente. La fórmula BDF de k = 1 paso, es el método de Euler implícito yn+1 − yn = hn f (tn+1 , yn+1 ). Siempre que k ≤ 6, la fórmula BDF de k pasos es convergente de orden k. Igual que los métodos de Adams, se utilizan estrategias de orden y paso variables. Las fórmulas BDF son en general menos eficientes que los métodos de Adams (y por supuesto, que los pares encajados modernos de métodos RK) para problemas no stiff. Por eso, sólo se emplean en problemas stiff. El comando de Matlab que utiliza estas fórmulas es ode15s. En realidad, utiliza unas fórmulas similares a las BDF, debidas a Klopfestein y Reiher en 1971 y 1978, denominadas NDF1 (del inglés, “numerical differentiation formulae") y supuestamente más eficientes que las fórmulas BDF. Dado que a día de hoy esa mejor eficiencia no es manifiestamente clara, el comando ode15s permite utilizar las fórmulas BDF siempre que inicialize la opción BDF con el valor ’on’ mediante el comando odeset. Ejemplo 1 Retomemos el sistema del Brusselator visto al principio de la lección, u′ = A + u2 v − (B + 1)u, v ′ = Bu − u2 v con A = 1 y B = 3, en el intervalo [0, 20], y condición inicial     u(0) 1,5 y0 = = . v(0) 3

(8)

(9)

Si abreviamos (8) como y ′ = f (t, y) (con y = [u, v]T ) la función f del segundo miembro es obviamente       u A + u2 v − (B + 1)u, f1 (u, v), , , y= = f (t, y) = f (u, v) = 2 v Bu − u v f2 (u, v) 1

Puede ver una descripción de las fórmulas NDF en [7]

5

y la matriz jacobiana fy es fy =



∂f1 ∂u ∂f2 ∂u

∂f1 ∂v ∂f2 ∂v



=



2uv − (B + 1) B − 2uv

u2 −u2



.

Si queremos emplear el comando ode15s debemos elaborar una función que dados t e y devuelva esta matriz jacobiana. Para los valores de A y B empleados hasta ahora, A = 1 y B = 3, dicha función es function J=brussjac(t,y) % Matriz jacobiana del segundo miembro del Brusselator p=2*y(2)*y(1);q=y(1)*y(1); J=[p-4, q; 3-p, -q]; Para utilizar las fórmulas BDF, en las correspondientes opciones debemos indicar el nombre del fichero que contiene esta función (en nuestro caso brussjac.m) sin la extensión .m. Inicializamos las opciones para utilizar ode15s, opbdf=odeset(’AbsTol’,1e-8,’RelTol’,1e-8,’Stats’,’on’) opbdf=odeset(opbdf,’BDF’,’on’,’Jacobian’,@brussjac) En la primera llamada al comando odeset hemos inicializado opciones ya conocidas, y en la segunda hemos indicado que ode15s debe utilizar las fórmulas BDF inicializando con ’on’ la opción ’BDF’, y hemos indicado que la función que calcula la matriz jacobiana lleva por nombre brussjac.m. Ejecutamos [Tb8,Yb8]=ode15s(@bruss,tiempos,[1.5 3]’,opbdf) y nos responde 901 successful steps 40 failed attempts 1508 function evaluations 1 partial derivatives 123 LU decompositions 1506 solutions of linear systems Podemos observar lo siguiente. En los 901+40 pasos que ha dado y ha aceptado o rechazado, ha resuelto 1506 sistemas lineales. Tal y como señalamos anteriormente, para dar cada paso debe resolver un sistema de ecuaciones no lineales por el método de Newton, resolviendo los sistemas lineales en (7). Vemos que en la mitad de los pasos sólo ha efectuado una iteración. En los 901+40 pasos y rechazos, sólo ha efectuado 123 descomposiciones LU. Esto es debido a que la matriz en el sistema (7), una vez hecha su factorización LU (que recordemos tiene un costo que crece como el cubo de la dimensión de la matriz), aprovecha dicha descomposición no sólo en otras iteraciones del método de Newton, sino en otros pasos. 6

La razón de los ahorros en iteraciones del método de Newton y de factorizaciones LU se debe [0] a que los iterantes iniciales yn+1 para el método de Newton obtenidos por extrapolación son tan correctos, que una sóla iteración del método de Newton y con una matriz que no es la correcta es suficiente para encontrar el valor correcto de yn+1 . El número de pasos, 901, es muy superior al del par encajado DOPRI54 para las mismas tolerancias, que vimos en secciones anteriores. Esto es, el método DOPRI54 puede dar pasos considerablemente mayores que las fórmulas BDF para una misma tolerancia. El gráfico en la Fig. 1, donde se mide la eficiencia relativa de ambos métodos, muestra que, para alcanzar un nivel de precisión determinado, las fórmulas BDF son un promedio de tres veces más caras que el método DOPRI54. Aunque se trata sólamente del ejemplo que venimos tratando, esto −2

10

−4

BDF

10

−6

error

10

−8

10

−10

DOPRI54

10

−12

10

−1

0

10

10

1

10

cputime

Figura 1: Eficiencia de los métodos DOPRI54 y fórmulas BDF en problema (8–9) es, el problema (8–9) asociado al Brusselator, esta situación es típica y se repite en otros problemas, salvo que éstos, como veremos en la sección siguitente, sean stiff.

1.3.

Otros métodos

Dado que los puede encontrar en ciertos contextos determinados así como en diversos libros de texto, comentamos aquí algunos métodos lineales multipaso que no pertenecen a las familias descritas en los apartados anteriores. En lo que sigue y por brevedad, denotaremos fn = f (tn , yn ),

n = 0, 1, . . . N.

7

Para los métodos que siguen a continuación, sólo consideraremos el caso hn = h,

n = 0, 1, . . . , N − 1.

1. Regla de Simpson, Es de la forma yn+1 − yn−1 =

h (fn+1 + 4fn + fn−1 ), 3

n = 1, 2, . . . , N − 1,

y es convergente de orden 4. Note que debido a su alto orden el valor de arranque y1 debe calcularse con un orden similar (por ejemplo con un método Runge-Kutta). Proviene de aplicar la regla de Simpson en la igualdad Z tn +h y(tn + h) − y(tn − h) = f (t, y(t)) dt. tn −h

2. Regla del punto medio explícita o “Leap-Frog”. Es de la forma, yn+1 − yn−1 = 2hfn ,

n = 1, . . . , N − 1.

Es de orden 2. El valor de arranque y1 se suele aproximar con el método de Euler. 3. Método centrado para ecuaciones de segundo orden. Aunque la ecuación (o sistema) de segundo orden y ′′ = f (t, y), se puede escribir como un sistema de primer orden de doble dimensión,   ′   z y = , f (t, y) z′ y tratarlo como tal con cualquiera de los métodos vistos en las secciones anteriores, existen métodos específicos para sistemas de segundo orden. Aunque quedan fuera del alcance de este curso, mencionamos aquí uno que aparece frecuentemente en la aproximación numérica de ecuaciones en derivadas parciales hiperbólicas. yn+1 − 2yn + yn−1 = h2 fn ,

n = 1, . . . , N − 1.

Es de orden 2. El valor de arranque y1 se suele aproximar como y1 = y0 + hy0′ +

2.

h2 f0 . 2

Problemas stiff

Reciben este nombre determinados problemas de valor inicial, donde, para una precisión determinada, un método explícito (como los métodos RK vistos hasta ahora) necesitan para integrarlos tomar pasos mucho más pequeños que lo que la regularidad de la solución demanda. A continuación veremos un ejemplo ilustrativo, pero antes una cuestión de nomenclatura. La palabra “stiff” significa 8

rígido en inglés, pero no se ha traducido al referirnos en español a este tipo de problemas. De hecho, se toma como cursi o ignorante hablar de “problemas rígidos”. Quizá, el ejemplo más sencillo para aclarar la diferencia entre problemas stiff y aquellos que no lo son sea el original de Prothero y Robinson, de 1974. y ′(t) = λ(y(t) − g(t)) + g ′(t), y(0) = g(0).

(10) (11)

Las soluciones de la ecuación (10) son de la forma y(t) = αeλt + g(t),

α ∈ R,

(12)

(o α ∈ C)

y la solución del problema (10–11) es y(t) = g(t). Se supone que g tiene infinitas derivadas y que el tamaño de estas es moderado. Cuando Re(λ) < 0 y |Re(λ)| es muy grande en relación con |Im(λ)| y el tamaño de las derivadas de g, el problema (10–11) es stiff. Ejemplo 2 Para g(t) = cos(t), que es una función cuyas derivadas son del orden de la unidad, consideramos el problema (10–11), con λ = −1 y λ = −100. La Fig. 2 muestra el error frente al BDF (λ=−1)

−4

10

−6

error

10

−8

10

DOPRI54 (λ=−1)

−10

DOPRI54 (λ=−100)

10

−12

BDF (λ=−100)

10

−1

10

0

1

10

10

2

10

cputime

Figura 2: Eficiencia computacional de los métodos DOPRI54 (–) y fórmulas BDF (--) para el problema de Prothero-Robinson (10–11), con g(t) = cos(t), para λ = −1 (*), y para λ = −100 (+). costo para diversos valores T OL = RT OL = 10−4 , 10−6, . . . , 10−12 de los métodos DOPRI54 (en línea continua) y las fórmulas BDF (en línea de trazos) para los casos λ = −1 (marcado con ∗) y λ = −100 (marcado con +). Son oportunas ciertas observaciones que resumimos a continuación. 1. Para λ = −1 (gráficas con *), los resultados son similares a los del problema (8–9) asociado al Brusselator que vimos en la Fig. 1 de la sección anterior, esto es, que para alcanzar un niver de precisión determinado, los fórmulas BDf son tres veces más caras que el DOPRI54. 9

2. Para λ = −100, aunque se trata de la misma solución y(t) = cos(t) del caso λ = −1, la situación es diametralmente opuesta, resultando ser el método DOPRI54 unas 8 veces más costoso que las fórmulas BDF. 3. Los resultados del DOPRI54 muestran que dada una tolerancia T OL, en ambos problemas obtiene los mismos errores, pero para λ = −100 se obtienen con 15 veces más de costo. 4. Los resultados de las fórmulas BDF muestran que dada una tolerancia T OL, en ambos problemas tardan lo mismo en integrar la solución, pero para λ = −100, los errores obtenidos son unas 40 veces más pequeños. A la vista de estas observaciones, el hecho más destacado es el deterioro tan significativo de la eficiencia del método DOPRI54 al pasar de λ = −1 a λ = −100. Este deterioro no puede ser debido a la solución del PVI, pues en ambos casos es la misma, y(t) = cos(t). Debe ser por tanto debido al PVI (ecuación diferencial junto con su condición inicial). La razón de este deterioro es que DOPRI54 se ve forzado a tomar pasos muy pequeños para poder integrar el problema. Por ejemplo, para T OL = 10−8 , la información proporcionada por Matlab es la siguiente. Para λ = −1, 203 successful steps 7 failed attempts 1261 function evaluations 0 partial derivatives 0 LU decompositions 0 solutions of linear systems y para λ = −100, 3457 successful steps 11 failed attempts 20809 function evaluations 0 partial derivatives 0 LU decompositions 0 solutions of linear systems Vemos que efectivamente el cómputo total de pasos aceptados y rechazados se multiplica por 16 al pasar de λ = −1 a λ = −100, y en igual proporción el número de evaluaciones de función. No hay otra descripción de lo que es un problema stiff que la que ya hemos dado y repetimos a continuación. Un problema, para una precisión determinada y en un intervalo de tiempo determinado, es stiff cuando un método explícito debe tomar pasos mucho más pequeños que lo que la regularidad de la solución demanda. Comentamos las palabras en cursiva. Cuando decimos la regularidad de la solución, esto quiere decir que si esa misma solución lo fuese de otro problema, el método tomaría pasos considerablemente más largos manteniendo el nivel de precisión pedido. (Compare si no los pasos que ha necesitado el método DOPRI54 para integrar y(t) = cos(t) según sea λ = −1 o λ = −100). 10

Cuando un problema es stiff, todos los métodos explícitos sufren la misma pérdida de eficiencia. Algunos métodos implícitos también (como por ejemplo los métodos de Adams implícitos) aunque no otros (como las fórmulas BDF). Por eso se utiliza como criterio para determinar si un problema es stiff el comportamiento en cualquier método explícito. No es infrecuente que un problema sea stiff en un intervalo [t0 , tf ] pero deje de serlo a partir de t0 (o caso más frecuente) no lo sea antes de t0 . El caracter stiff de un problema depende de la precisión deseada. Si demandamos mucha precisión (por ejemplo, en un ordenador cuya unidad de redondeo sea del orden de 10−36 ) cualquier método se ve obligado a tomar pasos pequeños para conseguir errores muy bajos, y los problemas dejan de ser stiff. Por último, para hacerse una idea de cuál es la dificultad de los problemas stiff, basta considerar en la ecuación (10) la desviación con respecto a la solución particular z(t) = y(t) − g(t). Note que al ser Re(λ) < 0, y que tal y como vimos en (12), se tiene que (13)

l´ım z(t) = 0.

t→∞

Notemos que z(t) satisface la ecuación z ′ = λz(t). Para el método de Euler explícito (que es el único método Runge-Kutta de una etapa convergente), y sobre una partición uniforme (hn = h, para n = 1, 2, . . . ) se tiene que zn+1 = zn + λhzn = (1 + λh)zn = (1 + λh)2 zn−1 = . . . = (1 + λh)n+1 z0 . Tenemos por tanto que con el metodo de Euler la solución numérica zn reproduce el comportamiento de la solución de verdad en (13), l´ım zn = 0

n→∞

⇔

|1 + λh| < 1.

En el caso de λ ∈ R (y dado que entonces λ < 0) tenemos que l´ım zn = 0

n→∞

⇔

h<

1 . |λ|

Dado que para otros métodos RK la situación es similar, esto explica las diferencias entre λ = −1 y λ = −100 en el Ejemplo 2. Por otro lado para el método de Euler implícito (que es la fórmula BDF de un paso) tenemos que zn+1 = zn + λhzn+1 , esto es, (1 − λh)zn+1 = zn , o, equivalentemente, zn+1 = (1 − λh)−1 zn = (1 − λh)−2 zn−1 = . . . = (1 − λh)−n−1 z0 . 11

Tenemos por tanto que con el metodo de Euler implícito, l´ım zn = 0

n→∞

⇔

1 < 1. |1 + λh|

En el caso de λ ∈ R y negativo tenemos que l´ım zn = 0

⇔

n→∞

h > 0,

esto es, que para reproducir el comportamiento de la solución z(t) no hay ninguna restricción de longitud de paso, salvo la que venga impuesta por la regularidad de la misma. El caso de otras fórmulas BDF de más pasos es similar al del método de Euler implícito, explicando esto que las fórmulas BDF sean insensibles al tamaño de λ en la ecuación (10). Puede encontrar más información sobre la naturaleza de los problemas stiff en [3], § 6.1–6.2, así como amplia documentación sobre este tipo de problemas y su tratamiento numérico en [2]

3.

Otras cuestiones de orden práctico

Los comandos de Matlab para integrar EDOs o sistemas de EDOs, a través de las opciones que les demos mediante el comando odeset, pueden llevar a cabo diversas tareas como dibujar mapas de fase, ayudarnos a determinar el periodo de una solución periódica, etc. Puede que le resulte de utilidad la posibilidad de construir mapas de fase. Mediante la opción OuputFcn podemos indicarle a Matlab que realice determinada acción cada vez que calcula un punto. De hecho, Matlab tiene prevista la función odephas2 que dibuja las trayectorias según las va calculado. El siguiente mapa de fase 1.5

1

0.5

0

−0.5

−1

−1.5 −2

−1.5

del sistema

se obtuvo con la secuencia de comandos

−1

−0.5

0

x′ = 1 − y y ′ = x2 − y 2

12

0.5



1

1.5

2

(14)

% en la opci\’on para la funci\’on de output, damos odephas2 para % que dibuje la curva param\’etrica t -> [x(t),y(t)]’; opciones=odeset(’AbsTol’,1e-6,’RelTol’,1e-6,’OutputFcn’,@odephas2); figure(1);clf % dibujamos los ejes en rojo plot([-2,2],[0,0],’r:’) hold plot([0,0],[-1.5 1.5],’r:’); % fijamos los extremos del dibujo axis([-2 2 -1.5 1.5]) % por cada condici\’on inicial, que avance 5 unidades de tiempo % y que retroceda cinco unidades de tiempo. tiempos=[0 5]; tiempos2=[0 -5]; for j=1:25 % seleccionamos las condiciones inciales con el rat\’on g=ginput(1) % marcamos la condici\’on incial en rojo y con un circulito. plot(g(1),g(2),’ro’) % integramos el sistema (dado que la opci\’on OutputFcn tiene % el valor @odephas2, se dibujar\’an las trayectorias seg\’un se % calculan [T,Y]=ode45(@sistema,tiempos,g,opciones); [T,Y]=ode45(@sistema,tiempos2,g,opciones); end donde la función sistema proporciona el segundo miembro del sistema (14): function f=sistema(t,y) f=[1-y(2);-sum(y)*diff(y)]; Claramente se aprecian dos equilibrios, uno de ellos un punto de silla, y el otro un punto espiral. Puede continuar el estudio de este ejemplo en el Problema ??.

4. 4.1.

El método del disparo para problemas de contorno Cuestiones de repaso

Recordemos que un problema de contorno viene dado por una ecuación diferencial de segundo orden (o de orden par, en general) para la cual se busca la solución u en un intervalo [0, l] o [a, b], en la que en vez de imponer condiciones iniciales en 0 o en a, se imponen tantas condiciones de contorno en 0 o a como en l o b. Un ejemplo puede ser el problema  −kuxx (x) = f (x), x ∈ (0, l), (15) u(0) = α, u(l) = β, 13

Las condiciones de contorno tanto en el extremo x = a como en el extremo x = b pueden ser de uno de los tres tipos siguientes, donde entnderemos que x0 puede ser o bien a o bien b. i) Condición de tipo Dirichlet, u(x0 ) = u0 ii) Condición de tipo Nuemann, βux (x0 ) = v0 , donde β es un valor real. iii) Condición mixta o de tipo Robbin, αu(x0 ) + βux (x0 ) = w0 , donde α y β son valores reales.

Ejemplo 3 Una viga colocada horizontalmente y con una densidad de carga lateral f (x), empotrada en un extremo y apoyada en el otro, f(x) x

sufre desplazamiento transversal w(x) en la dirección vertical con respecto su fibra neutra.

f(x) x w(x)

Dicho desplazamiento w(x) satisface, EIwxxxx = −f (x), w(0) = w′ (0) = 0, w(l) = w′′ (l) = 0,

 x ∈ (0, l), 

(16)



donde E es el módulo de Young e I el momento de inercia de la sección transversal con respecto al plano horizontal (Si la viga no es de sección constante, debe entonces cambiarse la ecuación EIwxxxx = −f (x) por (EI(x)wxx (x))xx = −f (x)). Puede obtener más información [4], § 26.6 y [5], § 7.2.

Los problemas de contorno son lineales si sólo aparecen expresiones lineales de la incógnita y sus derivadas. Tal es el caso de los problemas que hemos escrito hasta ahora. Son lineales y de coeficientes constantes si aparte de ser lineales los coeficientes que multiplican a la incógnita u y sus derivadas no dependen de la variable independiente x. Los problemas en los que nos centraremos serán lineales, aunque deber saber que también los hay no lineales, como por ejemplo el problema de Bratu,  uxx + λeu = 0, x ∈ (0, 1), (17) u(0) = 0, u(1) = 0, donde λ es un parámetro real, o el siguiente problema

u′′ + |u| = 0, u(0) = 0, u(l) = β, 14



,

(18)

donde l > 0. La casuística de los problemas de contorno es bastante más compleja que las de los sistemas de ecuaciones ordinarias. √ Así por ejemplo, el problema (17) tiene dos soluciones distintas para λ ∈ (λ, λ0 ) donde λ0 ≈ 2, una única solución si λ = λ0 y ninguna si λ > λ0 . A su vez, el problema (18) si l > π, tiene dos soluciones distintas si β < 0, una única solución si β = 0, y ninguna si β > 0; En cambio, los problemas (15) y (16) tienen una única solución, sea quien sea f en ambos casos y sean quienes sean α y β en el problema (15). Esta diferencia entre los problemas (17) y (18) por un lado, y (15) y (16) por otro, no se debe a que estos sean lineales de coeficientes constantes y los otros dos sean no lineales. Los problemas lineales y de coeficientes constantes pueden tener también una casuística variada. Así por ejemplo, el problema  y ′′ + y = cos(t), , (19) y(0) = 0, y(π) = 0, tiene infinitas soluciones, pero el problema y ′′ + y = cos( 2t ), y(0) = 0, y(π) = 0,



(20)

,

no tiene ninguna, y ambos son lineales y de coeficientes constantes. Recordemos que para problemas lineales (no necesariamente de coeficientes constantes) la casuística es algo más simple y viene dada por el Teorema de la Alternativa de Fredholm que enunciaremos para un problema general de la forma  uxx + b(x)ux (x) + c(x)u(x) = f (x), x ∈ (0, l) (PC) (21) B0 (u) = a0 , Bl (u) = al , donde    B0 (u)  

Bl (u)

,

puede ser cualquiera de



  u(0), 

u(l),

ux (0),

o

ux (l),

o

 −ux (0) + α0 u(0),  ux (l) + αl u(l),



con α0 6= 0 y αl 6= 0. Recordando que asociado a este problema tenemos el correspondiente problema homogéneo  uxx + b(x)ux (x) + c(x)u(x) = 0, x ∈ (0, l) (PCH). (22) B0 (u) = 0, Bl (u) = 0, Tenemos el siguiente resultado. Teorema 1 Teorema de la Alternativa (de Fredholm). Se verifica una de las dos opciones excluyentes siguientes. i) Para toda f continua el problema (PC) tiene una única solución. ii) El problema (PCH) tiene solución no trivial. Además, en caso de que ii) sea cierto, si el problema (PC) tiene solución (puede no tenerla), el conjunto de soluciones estará constituido la suma de una solución cualquiera de (PC) y todas las de (PCH). 15

Nota 1 El teorema anterior también es válido cuando p y q son periódicas de periodo l y se pide que la solución u sea también periódica mediante las condiciones de contorno B0 (u) = u(l) − u(0) = 0,

B1 (u) = ux (l) − ux (0) = 0.

Nota 2 En el caso ii) del teorema anterior, las soluciones de (PCH) son todas proporcionales a una dada, salvo en el caso de condiciones frontera periódicas que pueden ser las combinaciones lineales de dos soluciones independientes. En el caso de darse la opción ii) del teorema anterior veremos en la próxima sección un criterio simple para averiguar si el problema de contorno tiene solución (en cuyo caso tendrá infinitas) o, por el contrario, no posee ninguna solución.

4.2.

El método del disparo

Recuerde que en el método del disparo para un problema de contorno lineal de la forma  y ′′ + py ′ + qy = r, (PC), y(a) = α, y(b) = β,

(23)

donde p = p(t), q = q(t) y r = r(t) son funciones continuas en todo [a, b], expresamos la solución (en caso de haberla) como combinación líneal donde γ =

y(t) = yp (t) + γy0 (t),

β − yp (b) , y0 (b)

(24)

de las soluciones yp e y0 de dos problemas de valor inicial, yp′′ + pyp′ + qyp = r, yp (a) = α, yp′ (a) = 0,



(25)

.

(26)

y y0′′ + py0′ + qy0 = 0, y0 (a) = 0, y0′ (a) = 1



Recuerde además que si al calcular la constante γ en (24) el denominador satisface que y0 (b) = 0, entonces y0 es solución no trivial (pues y0′ (a) = 1 6= 0) del problema homonéneo asociado a (23), y por tanto éste puede no tener solución o tener infinitas. Este segundo caso se da, obviamente, cuando el numerador de γ, β − yp (b) es nulo (pues yp ya es entonces una solución de (23)). El primer caso se da cuando, siendo y0 (b) = 0, se tiene β − yp (b) 6= 0. Recuerde así mismo que cuando se utiliza este método en el ordenador se obtienen las funciones yp e y0 sobre una red de puntos a = t0 < t1 < . . . < tN = b. En el caso de problemas no lineales, que de forma general se escriben como  y ′′ + f (t, y, y ′) = r(t), . (27) y(a) = α, y(b) = β, el método de disparo consiste en resolver la ecuación g(s) = 0, 16

por algún método para resolver ecuaciones no lineales como puede ser el método de bisección o el método de la secante, donde g(s) = y(b) − β, e y la solución del PVI y ′′ + f (t, y, y ′) = r, y(a) = α, y ′(a) = s.



.

(28)

Para aquél valor de s para el que g(s) = 0, la correspondiente solución de (28) es la solución del problema de contorno (27). Nota 3 Aparte del método del disparo, existen otros métodos para problemas de contorno que trataremos más adelante en el curso. Nota 4 La mayor parte de los problemas de contorno que Vd. se pueda encontrar antes de abandonar esta escuela podrán resolverse sin dificultad con el método del disparo, o con los métodos que explicaremos en lecciones posteriores. No obstante, debe saber que Matlab dispone de un comando muy sofisticado para problemas de contorno, bvp4c. Es más complejo de usar que los comandos para PVI. Si se ve en la necesidad de utilizarlo, encontrará muy útil leer antes [8] (accesible por Internet en la dirección ftp://ftp.mathworks.com/pub/doc/papers/bvp), donde este comando es explicado con detalle por sus propios autores.

5.

Cuestiones y problemas CUESTIONES

Ejercicio 1 Obtener el método de Adams explícito de dos pasos, y particularizar el caso hn = h, para todo n. Ejercicio 2 Idem con método de Adams implícito de dos pasos. Ejercicio 3 Compruebe que efectivamente el método de Euler implícito es la fórmula BDF de un paso. Ejercicio 4 Obtener la fórmula BDF de dos pasos, y particularizar el caso hn = h, para todo n. Ejercicio 5 (Junio 2005). Determinar los valores α, β y γ para que el método de dos pasos
yn+2 − 2yn+1 + yn = h2 αf (yn+2) + βf (yn+1) + γf (yn ) , n = 0, 1, 2, . . . , N − 2,

(tn = t0 + nh, n = 0, 1, 2 . . . , N) sea consistente de orden lo más alto posible para ecuaciones de segundo orden de la forma y ′′ = f (y), esto es que para una solución de y ′′ = f (y) suficientemente regular
τn = y(tn+2 ) − 2y(tn+1 ) + y(tn ) − h2 αy ′′(tn+2 ) + βy ′′ (tn+1 ) + γy ′′ (tn ) sea un infinitésimo en h del orden lo más alto posible. Determinar dicho orden.

17

Ejercicio 6 (Junio 2005). La fórmula BDF de k pasos se obtiene imponiendo que el polinomio interpolador de yn−k+1, . . . , yn , yn+1 satisfaga la ecuación diferencial y ′ = f (t, y) en tn+1 . Para k = 2 (y hn = h para todo n), encuentre el método que se obtiene al imponer que polinomio interpolador de yn−1 , yn e yn+1 satisfaga la ecuación diferencial y ′ = f (t, y) en tn (observe que se pide en tn , y no en tn+1 como en las fórmulas BDF). Ejercicio 7 (Septiembre 2005). Los métodos de Nyström son parecidos a los de Adams pues, aparte de tener el mismo orden, en el método de Nyström explícito de k pasos, dados los valores yn−k+1, . . . , yn en los puntos tn−k+1 , . . . , tn (con tj = t0 + jh, j = 0, 1, 2 . . . ) la aproximación yn+1 se obtiene de la relación Z t n+1

yn+1 − yn−1 =

tn−1

∗ Pn,k (t) dt,

∗ donde, como en los métodos de Adams, Pn,k (t) es el polinomio de grado k − 1 que en las k abcisas tn−k+1 , . . . , tn toma los valores f (tn−k+1, yn−k+1), . . . , f (tn , yn ), respectivamente. Calcule el método de Nyström explícito de 3 pasos.

Ejercicio 8 Resolver −y ′′ + y = 0 sujeto a y(0) = 0 y y(π) = 2. Asimismo, determine los valores de a para los que el siguiente problema de contorno tenga solución.  ′′ y + y = t + asen(t) y(0) = y(π) = 0 Ejercicio 9 Si k(x) es continua y no nula en [0, l], encontrar dos soluciones independientes de (k(x)y ′ (x))′ = 0. Resuelva el problema de contorno  (k(x)y ′ (x))′ = 0 y(0) = 1, y(l) = 0 Ejercicio 10 Encuentre una función u y otra v que satisfagan u(0) = 1,

u′ (1) = 0,

v(0) = 1,

v ′ (1) + 2v(0) = 3.

Ejercicio 11 Basándose en el ejercicio anterior, pruebe que las solución de culquier problema de contorno lineal de segundo orden con condiciones fronteras no homegénes, se pueden encontrar resolviendo un problema con condiciones de contorno homogéneas. Ejercicio 12 Encuentre todas las soluciones de la ecuación −u′′ +u′ = 0 (ya sabe, las combinaciones lineales de las funciones eλx para λ raíz del polinomio p(x) = x2 −x). Determinar entonces los valores de α y β para los cuales el problema  −uxx + ux = 0, ux (0) = α, ux (l) = β, tiene solución.

18

Ejercicio 13 Resuelva el problema x2 uxx (x)+xux (x)−u(x) = 0, x ∈ (0, 1), sujeto a las condiciones de contorno u(1) = 1, y l´ım |u(x)| < ∞. x→0

(Encuentre primero los posibles valores de α para los cuales u = xα es solución de la ecuación diferencial x2 uxx (x) + xux (x) − u(x) = 0, y determine dos soluciones independientes en (0, 1)). Ejercicio 14 Verifique que para cualquier función continua h : [1, 4] → IR , el problema de contorno  2 ′′ 2x y (x) + 3xy ′ (x) − y(x) = h(x) y(1) = 0, y(4) = 1 tiene solución única. Para h(x) = x hallar la correspondiente solución Ejercicio 15 Resolver cuando sea posible el problema de contorno  ′′ y + 9y = t + a y(0) = 1, y(π) = 2 Ejercicio 16 Las fórmulas NDF son como las fórmulas BDF vistas en el curso a las que, para darles mayor estabilidad, se les añade un término proporcional a la diferencia dividida de orden k + 1 (éste término se añade a la fórmula de k pasos). 1. Para niveles de tiempos espaciados tj+1 − tj = h (j = n, n − 1, n − 2), calcule el desarrollo de Taylor en tn+1 de 6y[tn+1, tn , tn−1 , tn−2 ] cuya expresión es 6y[tn+1 , tn , tn−1 , tn−2 ] = y(tn+1) − 3y(tn ) + 3y(tn−1) − y(tn−2 ). 2. Determinar en función de τ el orden de la segunda fórmula NDF, esto es, el orden como infinitésimo en h de la expresión del error de truncación que es 3 1 τn = y(tn+1 ) − 2y(tn ) + y(tn−1 ) − hy ′ (tn+1 ) + 6y[tn+1 , tn , tn−1 , tn−2 ] 2 2

19

PROBLEMAS Problema 1 Consideremos el sistema autónomo no lineal  ′ x = 1−y y ′ = x2 − y 2 1. Determine analíticamente los puntos críticos del sistema y el tipo de su configuración. 2. Esboce el correspondiente plano de fases en el cuadrado [−2, 2] × [0, 2] utilizando la orden ode45 de Matlab. Para ello, diseñe un fichero que permita dibujar en una misma gráfica las órbitas que parten de los puntos (x, 0) con x = −2 : 0,5 : 2 y tales que el instante inicial es cero, el instante final es diez. Posteriormente, añada las órbitas que parten de los puntos (x, 2) con x = −2 : 0,25 : 2, con las mismas especificaciones. 3. Repita el apartado dos utilizando la orden ginput para seleccionar los puntos de arranque de las órbitas. Debe dibujar, al menos, treinta órbitas en el cuadrado [−2, 2]×[0, 2] y las restantes especificaciones son las mismas que en dicho apartado. 4. Repita el apartado dos usando la opción OutputFcn de la orden odeset de Matlab. 5. Para cada uno de los equilibrios del sistema, esboce un plano de fases en un rectángulo adecuado. Describa y justifique las representaciones gráficas obtenidas en este apartado y los tres anteriores. Problema 2 En este problema se estudia la integración numérica de dos ecuaciones diferenciales ordinarias rígidas (stiff) de primer orden. 1. Considere la ecuación diferencial y ′ = −1000(y − sen(t)) + cos(t). Obtenga analíticamente la solución exacta con valor inicial y(0) = 1. ¿Cuál es el comportamiento a largo plazo de dicha solución? 2. Utilizando la orden ode23 de Matlab, resuelva numéricamente el problema anterior en el intervalo [0, 1] y en el intervalo [0, 20]. Compruebe la exactitud de la aproximación y estime el número de evaluaciones funcionales requeridas en ambos casos. Utilizando el estilo de línea ’.’, describa como ha ido cambiando el paso de integración en ambos intervalos. ¿Cambia mucho la situación usando la orden ode45? ¿Y usando la orden ode113? 3. Repita el apartado dos, pero usando ahora la orden ode15s. Justifique los resultados obtenidos. 4. Integre numéricamente el problema 

y ′ = y 2 (1 − y) y(0) = 10−4

en el intervalo [0, 2 · 104 ] usando la orden ode45. Como en el apartado dos, comente las estadísticas de integración y la gráfica en estilo de línea ’.’ de la correspondiente solución. ¿Qué ocurre si repetimos el proceso en el intervalo [0, 104 ]?. 20

5. Repita el apartado cuatro usando la orden ode15s. Justifique los resultados obtenidos. Problema 3 (Junio 2005) En el curso hemos programado el método Runge-Kutta clásico. Parecido es el método Runge-Kutta-Nyström (RKN) que deberá programar ahora. Los métodos RKN son métodos de tipo Runge-Kutta para ecuaciones y sistemas de ecuaciones de segundo orden. Para el problema  y ′′ = f (y), t ∈ [t0 , tf ] (29) y(t0 ) = y0 , y ′(t0 ) = y0′ con un método RKN de s etapas se obtienen aproximaciones yn ≈ y(tn ) e yn′ ≈ y ′ (tn ) (donde tn = t0 + nh), por el siguiente procedimiento, ki = f (yn + ci hyn′ + h2

i−1 X

aij kj ),

i = 1, . . . , s,

j=1

yn+1 = yn + hyn′ + h2

s X

b i ki ,

i=1

′ yn+1 = yn′ + h

s X

bi ki

i=1

(note que hay coeficientes bi para las yn y bi para las yn′ ). Considere entonces el método de Nyström de s = 4 etapas de tablero de Butcher generalizado,

c A T b bT

0 1 5 2 3 dado por

1 bi bi

1 50 −1 27 21 70 14 336 14 336

7 27 −4 70 100 336 125 336

18 70 54 336 162 336

35 336

1. Elabore una función de Matlab que implemente dicho método de encabezamiento function [T,Y,Yp]=elrkn(laf,t0,tf,y0,yp0,h) % devuelve en T los los tiempos T(j)=t0+j*h y en cada fila de Y la % traspuesta de la aproximaci\’on en tiempo T(j) a la soluci’on % del problema y’’=laf(t,y), y(t0)=y0, y’(t0)=yp0 obtenida con el % m\’etodo Runge-Kutta-Nystrom de s=5 etapas con paso h, y en cada fila % de Yp la traspuesta de la aproximaci’on en tiempo T(j) de y’

21

2. Pruebe su función en el problema        ′′  y1 y1 (0) 1 y1 =− , = , y2′′ y2 y2 (0) 0



y1′ (0) y2′ (0)



=



0 1



,

t ∈ [0, 5].

¿Cuál es la solución de este problema? ¿Qué error obtiene con h = 0,1 en t = 5? 3. Averigue el orden p de convergencia del método. Explique cómo lo ha hecho. 4. Escriba con 4 cifras signicativas correctas la solución en tiempo t = 5 de     ′        ′′  −1 0 y1 (0) 1/4 y1 (0) y1 y1 = p , = , . = 2 y2′ (0) 0 y2 (0) y2′′ 5/3 (y1 + y22)3/2 y2

Problema 4 Considere el problema de contorno  ′′ y = 4y + 2t (P 1) ≡ y(0) = α, y(3) = β

1. Calcule la solución exacta del problema (P1) cuando α = 1 y β = 2. 2. Para obtener una aproximación a la solución del problema del apartado anterior sobre una red de espaciado h = 0,01 utilice el método del disparo, con el comando ode45 (con AbsTol y RelTol iguales a 10−8 ) para integrar los correspondientes problemas de valor inicial. Esboce una gráfica de la solución y calcule el máximo error (absoluto) cometido sobre los puntos 0, 0.01, 0.02, . . . , 1.99, 2. 3. Considere ahora el siguiente problema de contorno no lineal  ′′ y − 20yy ′ = 0, (P 2) ≡ y(0) = 1, y(1) = 0. Diseñe una función en Matlab que implemente el método de disparo no lineal, utilizando la orden ode45 para abordar los correspondientes problemas de valor inicial y el método de la secante para determinar la sucesión de parámetros. Los argumentos de entrada deben ser la correspondiente ecuación fun, el paso h para la red sobre la que se desea la solución, la tolerancia (absoluta y relativa) de la orden ode45, el máximo número de iteraciones y los dos valores iniciales s0 y s1 del método de la secante. Se supone que el problema de contorno se ha reformulado de la siguiente manera: determinar el parámetro s tal que  ′′ y = f (t, y, y ′), ; y(1) = 0. y(0) = 1, y ′(0) = s 4. Utilizando la función del apartado anterior, calcule una aproximación a la solución de (P 2) para h = 0,01, tolerancia 10−7 , s0 = 0, s1 = −10−8 y con un máximo de dos iteraciones (además de las asociadas a s0 y s1 ). Esboce una gráfica de la solución y escriba (con seis cifras significativas) el valor de y(1). Repita el proceso pero realizando esta vez cuatro iteraciones. 5. Proporcione el valor de la solución en x = 0,9. 22

Problema 5 1. Diseñe una función en Matlab que implemente el método de disparo no lineal, utilizando la orden ode45 para la resolución de los correspondientes problemas de valor inicial y el método de la secante para determinar la sucesión de parámetros. Los argumentos de entrada deben ser la correspondiente ecuación, los puntos a y b, las condiciones de contorno α y β, la tolerancia de la órden ode45, el máximo número de iteraciones y los dos valores iniciales s0 y s1 del método de la secante. Se supone que el problema de contorno se ha podido formular de la siguiente manera: determinar el parámetro s tal que  ′′ y = f (t, y, y ′), y y(b) = β. y(a) = α, y ′ (a) = s 2. Considere el problema de contorno (

y ′′ = y 3 − yy ′, 1 . 1 y(1) = , y(2) = 2 3

Obtenga una gráfica de la solución de dicho problema usando la función del apartado uno. Compare dicha gráfica con la gráfica de la solución exacta y(t) =

1 . t+1

3. Modifique el fichero del apartado uno para que el método de disparo no lineal incluya, además, como criterio de parada la diferencia en valor absoluto entre el valor obtenido y(sn ) en la iteración n-ésima y el valor deseado β. 4. Considere el siguiente problema de contorno  3 ′′′ t y − 4t2 y ′′ + 8ty ′ − 8y = 4 log(t), y(1) = 0, y ′ (1) = 1, y(2) = 1

.

Obtenga gráficamente una solución del problema anterior en un mallado de paso 0.1, usando la función del apartado tres con tolerancia (para la secante) igual a 10−4 . Calcule el error cometido sabiendo que la solución general de la ecuación es del tipo y(t) = c1 t + c2 t2 + c3 t4 −

1 7 log(t) − . 2 8

Problema 6 (Examan Final 2007-08) Considere el siguiente problema de contorno no lineal  ′′ y + 20y(1 − y 2 ) = 0, (P 1) ≡ y(0) = α, y ′(1) = β. 1. Diseñe una función en Matlab (o modifique la que hicimos en su día en las prácticas del curso) que implemente el método de disparo no lineal para este problema, utilizando la orden ode45 para abordar los correspondientes problemas de valor inicial y el método de la secante para determinar la sucesión de parámetros. Los argumentos de entrada deben ser la correspondiente ecuación fun, el paso h para la red sobre la que se desea la solución, los valores 23

α y β, la tolerancia (absoluta y relativa) de la orden ode45, el máximo número de iteraciones y los dos valores iniciales s0 y s1 del método de la secante. Se supone que el problema de contorno se ha reformulado de la siguiente manera: determinar el parámetro s tal que   y ′ = v, v ′ = f (t, y, v), ; v(1) − β = 0.  y(0) = 1, v(0) = s,

(Obviamente, f (t, y, v) = −20y(1 − y 2).) Escriba el listado de dicho programa en su hoja de respuestas.

2. Utilizando la función del apartado anterior, calcule una aproximación a la solución de (P 1) cuando α = 1 y β = −2; para h = 0,01, tolerancia 10−7 , s0 = 0, s1 = −1/100, y con un máximo de dos iteraciones (además de las asociadas a s0 y s1 ). Esboce una gráfica de la solución y escriba (con seis cifras significativas) el valor de y(1). Repita el proceso pero realizando esta vez cuatro iteraciones. 3. Proporcione el valor de la solución en x = 0,9 con 10 cifras decimales correctas. 4. Considere ahora el mismo problema de contorno pero cambiando la segunda condición de contorno a y ′(1) + y(1) = −2, esto es,  ′′ y + 20y(1 − y 2) = 0, (P 2) ≡ y(0) = 1, y ′(1) + y(1) = −2. Repita el apartado 2 para el problema de contorno con esta nueva condición frontera.

Referencias [1] E. Hairer, S. P. Nørsett & G. Wanner, Solving Ordinary Differential Equations I (2nd. Ed.) (Springer-Verlag, Berlin, 1993). [2] E. Hairer & G. Wanner, Solving Ordinary Differential Equations II (2nd. Ed.) (Springer-Verlag, Berlin, 1997). [3] J. D. Lambert, Numerical Methods for Ordinary Differential Systems. The initial value Problem. John Wiley & Sons, Chichester, 1991. [4] P. Puig Adam, Curso Teórico-Práctico de Cálculo Integral Aplicado a la Física y Técnica (15a Ed.) Ed. Biblioteca Matemática, Madrid, 1976. [5] P. Puig Adam, Curso Teórico-Práctico de Ecuaciones Difrenciales Aplicado a la Física y Técnica (16a Ed.) Roberto Puig Alvarez, Madrid, 1980. [6] J. M. Sanz-Serna, Diez Lecciones de Cálculo Numérico, Universidad de Valladolid, 1998. [7] L. F. Shampine & M. W. Reichelt, The Matlab ODE Suite, SIAM J. Sci. Comput., 18 (1997), pp. 1–22. 24

[8] L. F. Shampine, J. Kierzenka & M. W. Reichelt, Solving Boundary Value Problems for Ordinary Differential Equations in Matlab with bvp4c, manuscrito accesible electrónicamente en la dirección “ftp://ftp.mathworks.com/pub/doc/papers/bvp”, 2000. [9] G. F. Simmons Ecuaciones Diferenciales con Aplicaciones y Notas Históricas (Segund Edición), McGraw-Hill, Madrid, 1991. [10] G. W. Stewart, Afternotes on Numerical Analysis, SIAM, Philadelphia, 1996.

25