Matlab y La Solucion de Ecuaciones Diferenciales

ECUACIONES DIFERENCIALES

A. Cachafeiro L´ op ez – J. Il opez Ill´ l´ an Gon an onz´ z´ al ez alez Departamento de Matem´ Mat em´ ati ca Apl atica Aplic icada ada I

Universidade de Vigo

´ SEGUNDA EDICION A˜ no: no: 2009 Editores: Los autores ISBN: 978-84-611-9916-7

ÍNDICE GENERAL

1 Comandos

simb´ olicos para resolver ecuaciones. olicos Camp os de direcciones. Isoclinas ´ LICOS . . . . . . . . . . . . . . . . . 1.1 1.1 COMA COMAND NDOS OS SIMB SIMBO 1.1.1 1.1.1 Comand Comandos os para para la la deri deriv vaci´ acioń . . . . . . . . . . . . . 1.1.2 Problemas resueltos . . . . . . . . . . . . . . . . . . 1.1.3 El comando comando DSOLVE. DSOLVE. Resoluci´ Resoluci´ on de PVI . . . . . 1.1. 1.1.44 Ejem Ejempl plos os de EE EEDD DD que que no pued pueden en ser ser resu resuel elta tass por por 1.1. 1.1.55 Funci unci´ oń LAMBERTW(Z) . . . . . . . . . . . . . . 1.2 CAMPO DE DIRECCIONES DE UNA ED . . . . . . . . 1.3 TRAYECTORIAS ISOCLINAS . . . . . . . . . . . . . . . 1.3.1 Método odo de las isocl oclinas . . . . . . . . . . . . . . . .

2

. . . . . . . . . . . . . . . . . . . . . . . . DSOL DSOLVE VE . . . . . . . . . . . . . . . . . . . . . . . .

1

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

M´ eto etodo de separ se paraci aci´ ´ on on de variables. ariables. Ecuaciones Ecuaciones homo homog g´ enea e neas, s, exac exacta tas, s, impl impl´ ´ıcit ıcitas as y de Bern Bernou oull llii 2.1 2.2 2.3 2.4 2.4 2.5 2.5 2.6

ECUACION IONES EN VARIABLES SEPARADAS . . . ´ NEAS ECUACION ECUACIONES ES DIFERENCIA DIFERENCIALES LES HOMOG HOMOGE ECUACION IONES DIFERENCIALES EXACTAS . . . . ´ N DE CLAIRAUT . . . . . . . . . . . . . ECUA EC UACI CIO ´ N DE BERNOULLI . . . . . . . . . . . . ECUA EC UACI CIO FACTOR INTEGRANTE . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

11 . . . . . .

11 14 14 18 19 21

3 Resoluci´ on on simb´ olica de sistemas diferenciales ´ 3.1 SISTEM SISTEMA A LINEAL LINEAL HOMOG HOMOGENEO DE COEFICIENTES CONSTANTES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 SISTEM SISTEMA A LINEA LINEAL L COMP COMPLET LETO O DE COEFICIENTES CONSTANTES . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

i

. . . . . .

1 1 2 2 3 4 5 8 8

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

25 27

´ Indice General

ii ´ DE SISTEMAS DIFERENCIALES 3.3 RESOLUCION ´ LINEALES HOMOGENEOS DE COEFICIENTES CONSTANTES POR LA VIA ESPECTRAL . . . . . . . . . . . . . . 3.3.1 Cá lculo de valores y vectores propios via MATLAB . . . . . . 3.3.2 El método espectral . . . . . . . . . . . . . . . . . . . . . . . . 3.4 EEDD LINEALES DE ORDEN SUPERIOR . . . . . . . . . . . . . . 3.4.1 EEDD lineales de orden superior homogéneas . . . . . . . . . 3.4.2 EEDD lineales de orden superior completas . . . . . . . . . . ´ DE EULER Y DE EULER 3.5 ECUACION GENERALIZADA . . . . . . . . . . . . . . . . . . . . . . . . . . . . ´ 3.6 CALCULO DE WRONSKIANOS Y APLICACIONES . . . . . . . . 3.6.1 Wronskiano e independencia lineal . . . . . . . . . . . . . . . 3.7 SOLUCIONES PARTICULARES DE EEDD LINEALES COMPLETAS . . . . . . . . . . . . . . . . . . . . . . . . 3.7.1 Método de coeficientes indeterminados (MCI) . . . . . . . . . 3.7.2 Método de variaci´ on de parámetros (MVP) . . . . . . . . . . . ´ DE ORDEN . . . . . . . . . . . . . . . . . . . . . . . 3.8 REDUCCION 3.9 EJERCICIOS COMPLEMENTARIOS . . . . . . . . . . . . . . . . . ´ 3.10 APENDICE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.1 Notas sobre funciones especiales . . . . . . . . . . . . . . . . . 3.10.2 Sobre el método basado en el cálculo aproximado de las ra´ıces caracter´ıstica . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

30 30 31 32 32 32

. . . . . . . . . 35 . . . . . . . . . 35 . . . . . . . . . 36 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . de la ecuaci´ on . . . . . . . . .

Transformada de Laplace ´ Y CALCULO ´ 4.1 DEFINICION . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Definici´ on de la transformada de Laplace . . . . . . . . . . . . . 4.1.2 Cá lculo mediante comandos MATLAB. Ejemplos . . . . . . . . 4.1.3 Propiedades de la Transformada de Laplace . . . . . . . . . . . 4.2 TRANSFORMADA INVERSA DE LAPLACE . . . . . . . . . . . . . . ´ DE 4.3 APLICACIONES A LA RESOLUCION ´ INICIAL . . . . . . . . . . . . . . . . PROBLEMAS DE CONDICION 4.4 PROBLEMAS DE VALORES INICIALES CON DATOS DISCONTINUOS . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Función escalón unitario (Heaviside) . . . . . . . . . . . . . . . 4.4.2 Representaci´ o n de funciones con discontinuidades de salto finito 4.5 TRANSFORMADA DE LAPLACE DE UNA ´ PERIODICA ´ FUNCION . . . . . . . . . . . . . . . . . . . . . . . . . . ´ DE TRANSFERENCIA . . . . . . . . . . . . . . . . . . 4.6 LA FUNCION

38 38 40 43 44 45 45 45

49 . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. 49 . 49 . 50 . 51 . 52

. . . . . . . .

54

. . . . . . . . 56 . . . . . . . . 56 . . . . . . . . 58 . . . . . . . . . . . . . . . .

62 62

´ Indice general

5

M´ etodos num´ ericos para PVI (I) 5.1 FUNCIONES ESPECIALES . . . . . . . . . . . . . . . . 5.1.1 Funciones de Bessel . . . . . . . . . . . . . . . . . 5.2 ESTUDIO DE LAS SOLUCIONES DE UNA ED . . . . 5.2.1 Métodos numéricos de resolución de PVI . . . . . 5.2.2 Método de Euler o de la tangente . . . . . . . . . 5.2.3 Otra interpretaci´ on del método de Euler . . . . . 5.2.4 Ejemplos en los que se aplica el método de Euler 5.2.5 Precisi´ on de los resultados . . . . . . . . . . . . . ´ 5.3 METODO DE EULER MEJORADO . . . . . . . . . . . 5.3.1 La f´ ormula de los trapecios . . . . . . . . . . . . . 5.3.2 Deducci´ on del método de Euler mejorado . . . . 5.3.3 Euler versus Euler mejorado . . . . . . . . . . . . ´ ´ 5.4 OTROS METODOS NUMERICOS . . . . . . . . . . . . 5.4.1 Ejemplo de método de dos pasos . . . . . . . . . 5.5 CONCLUSIONES . . . . . . . . . . . . . . . . . . . . . . 5.6 EJERCICIOS ADICIONALES . . . . . . . . . . . . . . . 5.6.1 Aplicación: c´ alculo de la elasticidad b(x) . . . . .

6

65 . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

M´ etodos num´ ericos para PVI (II) ´ 6.1 METODOS PARA RESOLVER PROBLEMAS ASOCIADOS A SISTEMAS DE EEDD . . . . . . . . . 6.1.1 Método de Euler mejorado . . . . . . . . . . . . . 6.1.2 Sobre los programas S2 EULER y RK CELL . . . 6.1.3 Ejercicios con indicaciones para ser resueltos . . . 6.2 EJERCICIOS ADICIONALES . . . . . . . . . . . . . . . 6.2.1 Indicaciones para resolver los ejercicios (6.4)-(6.7) 6.3 SERIES DE FOURIER . . . . . . . . . . . . . . . . . . . 6.3.1 Programas . . . . . . . . . . . . . . . . . . . . . . 6.3.2 Ejercicios . . . . . . . . . . . . . . . . . . . . . . ´ 6.4 APENDICE . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.1 Programa SN EULER . . . . . . . . . . . . . . . 6.4.2 Interpolaci´ on . . . . . . . . . . . . . . . . . . . .

7

iii

M´ etodos num´ ericos para PVI (III) 7.1

65 66 67 68 68 69 70 75 78 78 79 79 82 82 83 84 84

87 . . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

87 87 88 89 94 95 99 99 99 101 101 102

105

´ METODOS DE RUNGE-KUTTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.1.1 Métodos de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.1.2 Método de Runge-Kutta de cuarto orden . . . . . . . . . . . . . . . . . . . . . 106

´ Indice General

iv

7.1.3 Descripció n del programa RUNGEKUTTA . . . . . . . . . . . . . . . . . . . . 106 7.1.4 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.1.5 Estabilidad de un método numérico . . . . . . . . . . . . . . . . . . . . . . . . 108 7.1.6 Descripció n del programa EJE7.4 . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.2 PROBLEMAS DE CONTORNO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.2.1 Planteamiento matem´ atico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.2.2 Aplicaci´ on: Viga columna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 7.3 ECUACIONES EN DERIVADAS PARCIALES . . . . . . . . . . . . . . . . . . . . . 111 7.3.1 Ecuació n del calor en una barra de longitud finita . . . . . . . . . . . . . . . . 111 7.3.2 Ecuació n de ondas sobre un intervalo acotado . . . . . . . . . . . . . . . . . . 114 7.3.3 Herramientas simb´ olicas para EDP . . . . . . . . . . . . . . . . . . . . . . . . 116

Bibliograf´ıa de pr´ acticas y sitios WEB

119

LISTA DE FIGURAS 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 1.10

Grá fico de las soluciones de los ejemplos 1.3 y 1.4 . . . . Campo de direcciones y soluci´ on particular. Ejemplo 1.8 Campo de direcciones y soluci´ on particular. Ejemplo 1.9 Representaci´ on (a) . . . . . . . . . . . . . . . . . . . . . Representaci´ on (b) . . . . . . . . . . . . . . . . . . . . . Ejemplo 1.10 . . . . . . . . . . . . . . . . . . . . . . . . Gráfico del ejemplo 1.11 . . . . . . . . . . . . . . . . . . Gráfico del ejemplo 1.12-a . . . . . . . . . . . . . . . . . Gráfico del ejemplo 1.12-b . . . . . . . . . . . . . . . . . Gr´ afico del ejemplo 1.13 . . . . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. 4 . 6 . 7 . 8 . 8 . 9 . 9 . 10 . 10 . 10

2.1 2.2 2.3 2.4 2.5 2.6

Soluci´ on del ejemplo 2.1 Gráfico del ejemplo 2.6 . Gráfico del ejemplo 2.7 . Gráfico del ejemplo 2.8 . Gráfico del ejemplo 2.9 . Gráfico del ejemplo 2.10

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

12 16 17 18 19 20

3.1 Evidencia gráfica de que Prueba(2) debe ser nula en [ 1, 1] . . . . . . . . . . . . . . . iv) iii) ii) i) 3.2 Gráficas de Z = y 1 6y1 + 9y1 24y1 20y1 en [ 2, 2] y [20, 60] . . . . . . . . .

30 47

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

− −

−

− − Representaci´ on gráfica de γ (t) en [−4, 4] . . . . . . . . . . . . . . . . . . . . . . . . .

4.1 52 4.2 Representaci´ on gráfica de la entrada g(t) y la salida I (t) . . . . . . . . . . . . . . . . 59 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8

Gráfica de J 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gráfica de J 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gráficas superpuestas de J 1 y J 0 . . . . . . . . . . . . . . . . . . . . Gráficas de Bessel J 3/2 y J 1/2 . . . . . . . . . . . . . . . . . . . . . Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . Interpretaci´ on del método de Euler . . . . . . . . . . . . . . . . . . Soluciones exactas y aproximadas de los ejemplos 5.1 y 5.2 . . . . . Soluciones exactas y aproximadas del ejemplo 5.2 para n = 100, 200

v

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

66 66 67 67 71 71 73 73

vi

Listado de figuras

5.9 5.10 5.11 5.12 5.13

Soluci´ on aproximada del ejemplo 5.3 . . . . . . . . . . . . . . Soluciones aproximadas del ejemplo 5.4 . . . . . . . . . . . . R egla del trapecio . . . . . . . . . . . . . . . . . . . . . . . . . Soluciones exactas y aproximadas del ejemplo 5.5 (n = 21, 41) Ilustraci´ on de un método de dos pasos . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. 74 . 77 . 78 . 81 . 83

6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11

Grá fico producido por S2 EULER al resolver el ejercicio 6.1 Grá fico producido por S2 EULER al resolver el ejercicio 6.2 Grá fico producido por S2 EULER al resolver el ejercicio 6.3 Gráfico obtenido al resolver el ejemplo 6.1 . . . . . . . . . . Grá ficas de S2 EULER.M e INTERPOLA.M. Ejercicio 6.4 . RK CELL e INTERPOLA.M aplicados al ejercicio 6.5 . . . Soluci´ on gráfica del ejercicio 6.6 . . . . . . . . . . . . . . . . Soluciones DSOLVE y RK CELL. Ejercicio 6.7 . . . . . . . . Ejercicio 6.8. Sumas parciales . . . . . . . . . . . . . . . . . Ejercicio 6.9. Sumas parciales N = 2, 26 . . . . . . . . . . . Ilustraci´ on gráfica del método de interpolación . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

90 91 92 93 96 96 97 98 100 100 103

7.1 Gráficas producidas por EJE7.4.M . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.2 Grá fica producida por CALOR.M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.3 Grá fica producida por ONDA.M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

INDICACIONES GENERALES Antes de comenzar la sesión de trabajo con MATLAB (ML) y Adobe-Reader (PDF) debemos personalizar el entorno de trabajo. 1. (ML) Seleccionar adecuadamente, para no da˜ nar la vista, el tipo, tamaño, color de las fuentes y el color de fondo de las ventanas MATLAB. 2. (ML y PDF) Ajustar el tama˜ no de ambas ventanas, Command Window de MATLAB y el guión de la clase, para trabajar simultáneamente con ambas. 3. (ML) Crear carpetas personales temporales y declarar la correspondiente trayectoria (PATH) mediante el menú principal: FILE\SET PATH. Esto es esencial para la ejecución de programas. Como alternativa basta ejecutar lo siguiente desde la ventana de comandos MATLAB >>!mkdir C:\carpeta_personal >>path(’C:\carpeta_personal’,path)

La primera l´ınea crea nuestra carpeta personal, mientras que la segunda informa a MATLAB dónde se encuentra ésta, y la sitúa al comienzo de la lista de trayectorias. 4. (ML) Para disminuir el espacio vertical inter-l´ınea en la ventana de comandos ejecutar >>format compact

5. (ML) Activar el interruptor MORE mediante >>more on

para impedir el desplazamiento descontrolado de la ventana de comandos. 6. (ML) Desactivar el interruptor WARNING mediante >>warning off

si se desea que el sistema no muestre en pantalla advertencias. El alumno debe entregar por escrito la solución de algunos ejercicios seleccionados, seg´ un se indica en la plantilla correspondiente. Para ello debe basarse en el desarrollo de los ejemplos que aparecen resueltos en el guión de la clase.

vii

´ CLASE PRACTICA 1 Comandos simb´ olicos para resolver ecuaciones. Campos de direcciones. Isoclinas

OBJETIVOS

• Abordar la solución de ecuaciones diferenciales (EEDD) y problemas de valor inicial desde

una perspectiva experimental, mediante la aplicación de los recursos simbólicos del sistema MATLAB.

• Contribuir con el uso de herramientas gráficas al estudio de las EEDD a través de la visualización del campo de direcciones asociado y las isoclinas.

1.1 1.1.1

´ COMANDOS SIMBOLICOS Comandos para la derivaci´ on

La siguiente secuencia halla la función derivada de f respecto a x. >> syms x >> diff(f,x)

Para hallar la derivada enésima hacemos >>diff(f,x,n)

El s´ımbolo Df representa simbólicamente la derivada de f . Para derivadas sucesivas usamos D2f, D3f, etc.

1

Clase Pr´ actica 1

2

1.1.2

Problemas resueltos

on Ejemplo 1.1 Dada la ecuaci´ my = comprobar que

−mg + ky ,

y(t) = A + Be kt/m + mgt/k, es una soluci´ on de dicha E.D

Soluci´ on del ejemplo 1.1 1. Definir simb´ olicamente y(t). 2. Sustituir en la ED >>syms t y k m g A B >>y = A + B*exp(k*t/m) + m*g/k*t; >>a=m*diff(y,t,2)-k*diff(y,t)+m*g a=1/m*B*k^2*exp(k*t/m)-k*(B*k/m*exp(k*t/m)+m*g/k)+m >>simplify(a) ans =0

on Ejemplo 1.2 Comprobar que la funci´ u(x, t) = e −π

2 2

c t/l2

sen (πx/l)

es soluci´ on de la ecuaci´ on del calor ut (x, t) = c 2 uxx(x, t)

Soluci´ on del ejemplo 1.2 >>syms c >>u=’exp(-pi^2*c^2*t/l^2)*sin(pi*x/l)’ >>diff(u,’t’)-c^2*diff(u,’x’,2) ans=0

olicas a las variables u, t, x y l. Nota 1.1 No ha sido necesario declarar simb´

1.1.3

El comando DSOLVE. Resoluci´ on de PVI

El comando DSOLVE es de tipo simbólico y permite obtener la solución exacta de algunas EEDD y de los correspondientes problemas de valor inicial.(1 El argumento de este comando es un STRING de modo que no es necesario declarar simbólicas a las variables involucradas. 1

En lo que sigue usaremos las siglas PVI para problema(s) de valor(es) inicial(es).

Clase Pr´ actica 1

3

Ejemplo 1.3 Resolver la ED y  = 1 + y 2 , obtener la soluci´ on local correspondiente a la C.I. (0, 1), y representarla gr´ aficamente. Resoluci´ on del ejemplo 1.3 Paso 1 >>dsolve(’Dy=1+y^2’) ans =tan(t+C1)

Paso 2 >>dsolve(’Dy=1+y^2’,’y(0)=1’) ans =tan(t+1/4*pi)

Paso 3 >>ezplot(’tan(t+pi/4)’,[-pi/4,pi/4]); >> hold on >> plot(0,1,’*r’)

Ejemplo 1.4 Resolver la ED y  = tarla.

−y − 2t. Obtener la soluci´ on local con la C.I. ( −2, 1) y represen-

Resoluci´ on del ejemplo 1.4 >>dsolve(’Dy=-y-2*t’) ans = -2*t+2+exp(-t)*C1 >>dsolve(’Dy=-y- 2*t’,’y(-2)=1’) ans = -2*t+2-5*exp(-t)/exp(2) >>ezplot(’-2*t+2-5*exp(-t)/exp(2)’,[-4,3]), >>hold on >>plot(-2,1,’*g’)

1.1.4

Ejemplos de EEDD que no pueden ser resueltas por DSOLVE

Ejemplo 1.5 Intentar resolver la ED y  =

−2t − y + cos(y) con DSOLVE.

Desarrollo del ejemplo 1.5 >>dsolve(’Dy=-2*t-y+cos(y)’) Warning: Explicit solution could not be found. >In C:\MATLABR11\toolbox\symbolic\dsolve.m at line 326 ans = [ empty sym ]

Si el comando DSOLVE no puede obtener expl´ıcitamente la solución, entonces el propio DSOLVE ofrece, en algunos casos, la solución en la forma impl´ıcita.

Clase Pr´ actica 1

4

Fig. 1.1: Gráfico de las soluciones de los ejemplos 1.3 y 1.4

Ejemplo 1.6 Intentar resolver y  = y 2 (1

− y) aplicando DSOLVE.

Resoluci´ on del ejemplo 1.6 >>dsolve(’Dy=y^2*(1-y)’) Warning: Explicit solution could not be found; implicit solution returned. > In C:\MATLABR11\toolbox\symbolic\dsolve.m at line 292 ans = t+1/y-log(y)+log(-1+y)+C1=0

Ejemplo 1.7 Resolver la ED y  =

−y + cos(2x), con las condiciones iniciales y(0) = 1, y(0) = 0.

Resoluci´ on del ejemplo 1.7 >> y=dsolve(’D2y=-y+ cos(2*x)’,’y(0)=1’,’Dy(0)=0’,’x’) y = (1/2*sin(x)+1/6*sin(3*x))*sin(x)+(1/6*cos(3*x)-1/2*cos(x))*cos(x)+4/3*cos(x)

Observar en las siguientes l´ıneas el efecto del comando SIMPLIFY. >> y=simplify(y) y= -2/3*cos(x)^2+4/3*cos(x)+1/3

1.1.5

Funci´ on LAMBERTW(Z)

MATLAB se vale de una colecci´ on de funciones básicas elementales más amplia que la que usualmente utilizamos, con la cual nos puede dar una mayor cantidad de respuestas en forma cerrada. Una de

Clase Pr´ actica 1

5

estas funciones es w(z ) = Lambertw(z ) definida impl´ıcitamente por w(z ) exp(w(z )) = z

∗

La definici´ on de cualquier función especial MATLAB puede ser consultada invocando en la l´ınea de comandos >>help+nombre

o bien, para aquellas espec´ıficas del n´ ucleo MAPLE >>mhelp+nombre donde nombre es el nombre de la función.

Otras funciones especiales incorporadas al sistema MATLAB-MAPLE son BERNOULLI, BESSElJ, DILOG, etc. Para ver el catálogo completo MAPLE ejecutamos la orden >>mfunlist

¿C´ omo evitar que Lambertw aparezca en la respuesta final? Un caso en que podemos eliminar a Lambertw de la expresión que hemos obtenido al ejecutar DSOLVE, es cuando podemos despejarle, es decir, cuando somos capaces de llegar a lambertw(E 1 (x)) = E 2 (x), donde E 1 y E 2 son ciertas expresiones en la variable x. Entonces, atendiendo a la formulación impl´ıcita que define a esta función, se tiene que E 2 (x)exp(E 2 (x)) = E 1 (x).

Ejercicio 1.1 Resolver las siguientes EEDD y PVI. a) (x2 + 1)dy + (y2 + 1)dx = 0, y(0) = 1/2 b) y  = y cos(x)/(1 + y 2 ), y(0) = 1

− 3y + 3y − y = x d) y  − 3y + 3y − y = x,  c) y 

y(0) = 0, y (0) = 1, y (0) = 0

1.2

CAMPO DE DIRECCIONES DE UNA ED

Dada la ED y  = f (x, y) , si y = y(x) es solución y (x0 , y0 ) es un punto de la gráfica de y(x), entonces la pendiente de y = y(x) en (x0 , y0 ) es y  (x0 ) que es igual a f (x0 , y0 ).

Clase Pr´ actica 1

6

Fig. 1.2: Campo de direcciones y solución particular. Ejemplo 1.8 Lo anterior significa que, para conocer la pendiente de la solución que pasa por (x0 , y0 ) no hace falta conocer dicha soluci´ on; basta con calcular f (x0 , y0 ). Si esto se hace con todos los puntos del plano, tendremos las pendientes de las soluciones que pasan por cada punto del plano. Naturalmente, en la práctica es imposible hacer esto con todos los puntos del plano, pero nada nos impide hacerlo con tantos puntos como queramos, configurando un gráfico que llamaremos campo de direcciones de la ED.

¿C´ omo construir “el campo de direcciones” asociado a una ED? Consideremos la ED genérica y  (x) = f (x, y(x)). Para construir su campo de direcciones procedemos de la siguiente forma: por cada punto (x, y) de una red de puntos de R2 , convenientemente prefijada, se dibuja un segmento de recta o vector de pendiente f (x, y). El resultado final de este trabajo se interpreta como información gráfica acerca de la solució n de la ED que aún no hemos intentado resolver por otros medios. Para aumentar la eficiencia usaremos un programa MATLAB en el que intervienen los comandos MESHGRID y QUIVER que permiten generar el gráfico del campo de direcciones. Los vectores a considerar para obtener el campo de direcciones serán (1, y  ) = (1, f (x, y)). on particular Ejemplo 1.8 Construir el campo de direcciones de la ED y  = x, y representar la soluci´ que cumple la condici´ on inicial y(0) = 0.

Resoluci´ on del ejemplo 1.8 Las acciones MATLAB son las siguientes: >>f=inline(’x’,’x’,’y’); >>paso=0.5;

Clase Pr´ actica 1

7

Fig. 1.3: Campo de direcciones y solución particular. Ejemplo 1.9 >>iz=-3; >>der=3; >>[x,y]=meshgrid(iz:paso:der,iz:paso:der); >>[n,m]=size(x); >>dx=ones(n,m); >>z=f(x,y); >>dy=z; >>quiver(x,y,dx,dy)

La figura 1.2 (izquierda) muestra el campo de direcciones. La solución particular se obtiene artesanalmente como y = x2 /2, mediante separación de variables, o bien, aplicando DSOLVE. A continuaci´ on aparecen las órdenes MATLAB necesarias (ver figura 1.2-derecha). >>hold on % HOLD ON permite superponer los graficos >>dsolve(’Dy=x’,’y(0)=0’, ’x’) >>ezplot(’1/2*x^2’,[-2.5,2.5]), >>plot(0,0,’*g’)

Ejemplo 1.9 Construir el campo de direcciones asociado a y  = particular que cumple la CI y(1) = 1.

−x/y y representar la soluci´ on

Indicaciones para el ejemplo 1.9 La red creada con MESHGRID para el ejemplo 1.9 debe evitar y = 0 y contener al punto (1, 1) de la condición inicial. La expresión dentro de INLINE que define a la función debe escribirse entre apóstrofes, es decir,  x./y . El punto debe preceder a /, y ^ pues se opera con vectores componente a componente. La solución gráfica está dada por las figuras en el cuadro 1.3.

−

∗

Clase Pr´ actica 1

8

Ejercicio 1.2 Asociar cada una de las representaciones (a) y (b), mostradas en los cuadros 1.4 y 1.5, con el campo de direcciones de alguna de las siguientes EEDD: i) y  = y

− x,

ii) y  = 2x.

Fig. 1.4: Representación (a)

1.3

Fig. 1.5: Representación (b)

TRAYECTORIAS ISOCLINAS

Definici´ on 1.1 Dada la ED y  = f (x, y) , se llaman curvas de nivel o isoclinas a las obtenidas al imponer la condici´ on y  = k.

1.3.1

M´ etodo de las isoclinas

El método de las isoclinas es una variante de las ideas antes descritas. Los puntos del plano por los que pasa una solución con pendiente k, son los puntos de la curva de ecuación f (x, y) = k (isoclina de pendiente k). Dibujando las distintas isoclinas se obtiene una representaci´ on similar a la del campo de direcciones. Puede tener interés identificar la isoclina para la pendiente 0 pues las soluciones tendr´ an generalmente un máximo o un m´ınimo al pasar por esta isoclina.

Ejemplo 1.10 1 Representar las isoclinas de la ED y  = x + y2 . 2 Representar el gr´ afico de contorno (curvas de nivel o isoclinas) de la superficie z = x + y 2 .

Clase Pr´ actica 1

9

Desarrollo de la soluci´ on del ejemplo 1.10 >>[x,y]=meshgrid(0:0.05:3,-2:0.05:2); >> z=x+y.^2; >> isoclinas=contour(x,y,z,20)

El resultado en pantalla es el que ofrece la figura 1.6.

Fig. 1.6: Ejemplo 1.10

Fig. 1.7: Gr´ afico del ejemplo 1.11

Ejemplo 1.11 Representar las isoclinas de la ED y  = x 2 + y2 . Soluci´ on, paso a paso, del ejemplo 1.11 (ver fig. 1.7) >>[x,y]=meshgrid(-4:0.05:4); >>z=x.^2+y.^2; >>isoclinas=contour(x,y,z,20);

Ejemplo 1.12 (a) Representar las isoclinas de la ED y  = 2x

− y.

(b) ¿Qué tipo de curvas son dichas isoclinas? (c) Representar las isoclinas correspondientes a k = 0, y k = 2. (d) ¿Qué particularidad tiene la correspondiente a k = 2?

Clase Pr´ actica 1

10

Fig. 1.8: Gráfico del ejemplo 1.12-a

Fig. 1.9: Gr´ afico del ejemplo 1.12-b

Ejemplo 1.13 Construir el campo de direcciones y las curvas de nivel de la ED y = sen (x) + y. Desarrollo, paso a paso, de la soluci´ o n del ejemplo 1.13 (figura 1.10)

Fig. 1.10: Gráfico del ejemplo 1.13

>>f=inline(’sin(x)+y’,’x’,’y’); >>paso=0.5;iz=-3;der=3; >>[x,y]=meshgrid(iz:paso:der,iz:paso:der); >>[n,m]=size(x);dx=ones(n,m); >>z=f(x,y);dy=z; >>hold on,contour(x,y,z,20); >>quiver(x,y,dx,dy);

´ CLASE PRACTICA 2 M´ etodo de separaci´ on de variables. Ecuaciones homog´ eneas, exactas, impl´ıcitas y de Bernoulli

OBJETIVOS 1. Ejercitar el uso de herramientas simb´ olicas para la resolución de ED. 2. Aplicar métodos de cuadraturas exactas estudiados en clase.

2.1

ECUACIONES EN VARIABLES SEPARADAS

Son de la forma y  (x) = f (x)h(y). Para resolver el PCI y  = f (x)h(y), y(x0 ) = y 0 , escribimos la ED como sigue:

dy = f (x)dx h(y)

y si llamamos g = h1 , entonces queda g(y)dy = f (x)dx. Integrando en ambos miembros de la igualdad anterior obtenemos G(y) = F (x) + C, (Solución general de la ED), y calculamos el parámetro C usando la CI. As´ı C = G(y0 )

− F (x ) 0

y por lo tanto G(y) = F (x) + G(y0 ) F (x0 ), (Solución del PCI). aficamente la soluci´ on. Ejemplo 2.1 Resolver el PVI, y  = 2xy , y(0) = 2, y representar gr´

−

−

11

Clase Pr´ actica 2

12

Fig. 2.1: Solución del ejemplo 2.1

Resoluci´ on del ejemplo 2.1 Resolveremos el PVI con MATLAB siguiendo el método de separaci´ on de variables. >>syms F G x y c >>f=’-2*x’; g=’1/y’; >>x0=0;y0=2; >>F=int(f); G=int(g); >>sol=solve(G-F-c,y)% solucion general >>c=solve(subs(sol,x0,’x’)-y0,c); >>y=subs(sol,c,’c’)% solucion particular

Para obtener el gráfico de la solución del ejemplo 2.1 (figura 2.1), ejecutar en la l´ınea de comandos. >>ezplot(y,[-2,2]), hold on, plot(0,2,’+r’)

Debemos distinguir entre la solución general sol =exp(-x^2+2*c)

y la solución particular del PVI y =exp(-x^2+log(2)). El ejemplo 2.1 puede resolverse directamente usando el comando DSOLVE, en efecto >>y=dsolve(’Dy=-2*x*y’,’y(0)=2’,’x’) y =2*exp(-x^2)

Ejemplo 2.2 Resolver el PVI, y  = sen (x)(y Resoluci´ on del ejemplo 2.2 >>f=’sin(x)’

− 1),

y(0) = 2.

Clase Pr´ actica 2

13

>>g=’1/(y-1)’ >>F=int(f,’x’,0,’X’) F =-cos(X)+1 >>G=int(g,’y’,2,’Y’) G =log(Y-1) >>sol=solve(G-F,’Y’) sol =exp(-cos(X)+1)+1

Ejemplo 2.3 Resolver y cos(x)dx

2

− (1 + y )dy = 0, y(0) = 1.

Resoluci´ on del ejemplo 2.3 Despejando obtenemos y  = y cos(x)/(1 + y2 ) y usando DSOLVE >>y=dsolve(’Dy=y*cos(x)/(1+y^2)’,’x’) y=-i*(-lambertw(exp(2*sin(x)+2*C1)))^(1/2)

La aceptación de la anterior respuesta dependerá del consumidor. Si tenemos en cuenta que Lambertw(x)=w

es la función definida impl´ıcitamente por wew = x, lo que se obtiene es y 2 exp(y2 ) = K exp(2 sen(x)).

∗

∗

∗

Ejemplo 2.4 Resolver la ED anterior, es decir, y cos(x)dx las variables.

2

− (1 + y )dy = 0, y(0) = 1, separando

Nota sobre el ejemplo 2.3 . Supongamos que G(Lambertw(x), y , x , C) = 0, es la solución general (con parámetro C ) , de una ED y  = f (x, y). Entonces podemos hallar una expresi´ on equivalente a la anterior si somos capaces de despejar a Lambertw(x) = w. En efecto, supongamos que Lambertw(x) = F (x). Sabiendo que Lambertw(x) = w se caracteriza por satisfacer la ecuación w(x)ew(x) = x, entonces la nueva expresión de la solución general es F (x)eF (x) = x. Resoluci´ on del ejemplo 2.4 Escribimos de una vez en la l´ınea de comandos >> x=solve(’int(cos(x),x)-int((1+y^2)/y,y)’) x = asin(1/2*y^2+log(y))

para obtener sen(x) = log(y) + y2 /2 + c, que es la solución general de la ED. Para hallar la solución correspondiente a las condiciones iniciales dadas hacemos >>c=simple(sym(’solve(subs(x=0,y=1,sin(x)=log(y)+1/2*y^2+c),c)’)) c=-1/2

Clase Pr´ actica 2

14 es decir, de nuevo obtenemos la solución particular: sen(x) = ln(y) + (y2

− 1)/2.

Ejercicio 2.1 Aplicar herramientas MATLAB para resolver los siguientes ejercicios extra´ıdos del Bolet´ın 2. a) y  + y2 sen (x) = 0 b) 4xy dx + (x2 + 1)dy = 0 c) xy = y + x2 y2 d) (x y)dy = (x + 3y)dx 2 2 e) (x + 1)dy = (y + 1)dx, y(0) = 1/2 f ) (1 + ex )yy  = e y , y(0) = 0

 −

2.2

−

´ ECUACIONES DIFERENCIALES HOMOGENEAS

Ejemplo 2.5 Resolver la ED (x2

2

− y )dx + xydy = 0.

Resoluci´ on del ejemplo 2.5 Es claro que es una ED homogénea pues P (x, y) = x 2 y 2 y Q(x, y) = xy son funciones homogéneas del mismo grado 2. Se sabe entonces que el cambio de variable z = y/x la transforma en una ED en variables separadas. Para resolverla usando objetos simbólicos del MATLAB se escribe la ED en la forma y + (x2 y2 )/(xy) = 0,

−

−

y se siguen los pasos siguientes:

>>syms x y z >>y=x*z; >>y=subs(y,’z(x)’,’z’) y =x*z(x) >>subs(diff(y,x)+(x^2-y^2)/(x*y),y,’y’) ans =z(x)+x*diff(z(x),x)+(x^2-x^2*z(x)^2)/x^2/z(x) La respuesta ans es la ED de variables separables que, para simplificar, resolvemos directamente

usando el comando DSOLVE. >>z=dsolve(’z+x*Dz+(x^2-x^2*z^2)/x^2/z=0’,’x’) z =[ (-2*log(x)+C1)^(1/2)] [ -(-2*log(x)+C1)^(1/2)]

Finalmente se deshace el cambio y se obtiene en forma impl´ıcita la solución (y/x)2 =

2.3

−2log(x) + C 1.

ECUACIONES DIFERENCIALES EXACTAS

Consideremos la una ED de orden uno, que puede escribirse en la forma diferencial, es decir P (x, y)dx + Q(x, y)dy = 0.

Clase Pr´ actica 2

15

Si P y Q tienen derivadas parciales primeras continuas en un abierto simplemente conexo D, entonces en D equivalen: 1. P (x, y)dx + Q(x, y)dy = 0 es exacta. 2. ∂P/∂y = ∂Q/∂x. Si la ED es exacta entonces existe el potencial F (x, y) tal que ∂F/∂x = P y ∂F/∂y = Q y la solución general es F (x, y) = C. Para la obtención de F se procede como sigue. Integramos a F (x, y) respecto a una sola de sus variables, digamos x. F (x, y) =



P (x, y)dx + f (y),

donde f (y) es una constante que depende de la variable y que ha permanecido constante durante la integraci´ on. Notar que para determinar a F (x, y) sólo resta hallar a f (y) Si derivamos con respecto a y



∂ ( P (x, y)dx) + f  (y). Q = ∂y Luego



∂ ( P (x, y)dx) f  (y) = Q . ∂y Integrando respecto a y la anterior expresión obtenemos f (y). La solución general es

−

∈ R.

F (x, y) = C, C

on Ejemplo 2.6 Comprobar que la ED (2x + y)dx + (x 3y)dy = 0, es exacta. Resolver dicha ecuaci´ y representar las soluciones.

−

Resoluci´ on del ejemplo 2.6 Utilizamos las herramientas Matlab de la siguiente manera. >>P=’2*x+y’; >>Q=’x-3*y’; >>test=diff(P,’y’)-diff(Q,’x’) test =0 >>F1=int(P,’x’); F1 =x^2+y*x >>derf=Q-diff(F1,’y’) derf =-3*y >>f=int(derf,’y’) f =-3/2*y^2 >>F=F1+f F =x^2+y*x-3/2*y^2

Clase Pr´ actica 2

16

Fig. 2.2: Gráfico del ejemplo 2.6 Por lo tanto la solución general es x 2 + xy 3/2y2 = C . Para representar gráficamente las soluciones, o curvas de nivel, escribimos (ver figura 2.2)

−

>>[x,y]=meshgrid(0:0.1:3); >>z=x.^2+y.*x-3/2*y.^2; >>contour(x,y,z,15)

Ejemplo 2.7 Comprobar que la ED y  = ciones.

x2 , es exacta. Resolver dicha ED y representar las soluy2

Resoluci´ on del ejemplo 2.7 >>P=’x^2’; >>Q=’-y^2’; >>test=diff(P,’y’)-diff(Q,’x’) test =0 >>F1=int(P,’x’) F1 =1/3*x^3 >>derf=Q-diff(F1,’y’) derf =-y^2 >>f=int(derf,’y’) f =-1/3*y^3 >>F=F1+f F =1/3*x^3-1/3*y^3

Por lo tanto, la solución general de la ED es 1/3x3 1/3y3 = C . Para representar las soluciones escribimos (ver figura 2.3)

−

Clase Pr´ actica 2

17

Fig. 2.3: Gráfico del ejemplo 2.7 >>[x,y]=meshgrid(-3:0.05:3); >>z=z=1/3*x.^3-1/3*y.^3; >>contour(x,y,z,15)

Ejemplo 2.8 Comprobar que la ED cosh(x) sen (y) + senh (x)cos(y)y = 0, es exacta. Resolver dicha ecuaci´ on y representar las soluciones. Resoluci´ on del ejemplo 2.8 >>P=’cosh(x)*sin(y)’; >>Q=’sinh(x)*cos(y)’; >>test=diff(P,’y’)-diff(Q,’x’) test =0 >>F1=int(P,’x’) F1 =sin(y)*sinh(x) >>derf=Q-diff(F1,’y’) derf =0 >>f=int(derf,’y’) f =0 >>F=F1+f o n del ejemplo 2.8 es sen(y) F =sin(y)*sinh(x) Por lo tanto, la soluci´

representar las soluciones escribimos (ver figura 2.4) >>[x,y]=meshgrid(0:0.05:2); >>z=sin(y).*sinh(x); >>contour(x,y,z,15)

∗ senh(x)

= C . Para

Clase Pr´ actica 2

18


2.4

´ DE CLAIRAUT ECUACION

Recordemos que la ecuación de Clairaut es de la forma y = xy  + f (y ).

(2.1)

El método de resolución es hacer y  = p y derivar respecto a x, teniendo en cuenta que p = p(x). Nos queda entonces la expresión siguiente dp (x + f  ( p)) = 0. dx Si

dp = 0 entonces y  = C y por tanto teniendo en cuenta (2.1) se obtiene que dx

∈ R,

y = C x + f (C ), C

(2.2)

siendo (2.2) un haz de rectas, todas ellas solución de la ecuación original (2.1). Si x + f  ( p) = 0, usando (2.1) se obtiene la solución singular en forma paramétrica: x = y =

−f ( p), −f ( p) p + f ( p).

En general no es necesario eliminar p para obtener una ecuación de la forma G(x, y) = 0, y de hecho, podr´ıa resultar muy dif´ıcil o imposible. on de Clairaut y = xy  Ejemplo 2.9 Resolver la ecuaci´ obtenida, as´ı como la soluci´ on singular.

−

1  2 (y ) . 4

Representar la soluci´ on general

Clase Pr´ actica 2

19


Resoluci´ on del ejemplo 2.9 >>sol=dsolve(’y=x*Dy-1/4*Dy^2’,’x’) sol = [ x*C1-1/4*C1^2] [ x^2]

La segunda función es la envolvente de la familia integral de la ecuación. Para dar una idea gráfica aproximada del conjunto solución hacemos (ver figura 2.5) >>x=(-4:0.1:4); >>y=x.^2; z=0; %envolvente y recta para C1=0 >>w=x.*2-1; %recta para C1=2 >>u=x.*8-16; %recta para C1=8 >>v=x.*4- 4; %recta para C1=4 >>t=x.*3-9/4; %recta para C1=3 >>j=x.*(-2)-1; %recta para C1=-2 >>plot(x,y,x,z,x,w,x,u,x,v,x,t,x,j) >>legend(’Envolvente’,’C1=0’,’C1=2’,’C1=8’,’C1=4’,’C1=3’,’C1=-2’)

on de Clairaut y = xy  +ln(y ). Ejercicio 2.2 Resolver utilizando herramientas MATLAB, la ecuaci´

2.5

´ DE BERNOULLI ECUACION

La llamada ecuación de Bernoulli es de la forma y + p(x)y = q (x)yn . on. Ejemplo 2.10 Resolver el PVI y  + y/x = y 2 ln(x), y(1) = 1, y representar la soluci´

Clase Pr´ actica 2

20

Resoluci´ on del ejemplo 2.10 >>sol=dsolve(’Dy+y/x=y^2*log(x)’,’y(1)=1’,’x’) sol =2/x/(-log(x)^2+2)

La gráfica en el intervalo [1, 3] se obtiene ejecutando (ver figura 2.6) >>ezplot(sol,[1,3])

Ejercicio 2.3 1. Resolver la ED 3xy

− 2y = 2. Resolver el PVI 3xy − 2y =

x3 . y2 x3 y2 ,

y( 1) = 1.

−

Fig. 2.6: Gráfico del ejemplo 2.10 La ED de Bernoulli y + p(x)y = q (x)yn se resuelve transformándola en lineal, para lo cual se procede como sigue: Se divide la ED por y n y se hace el cambio z = y 1−1 , obteniendo una ED lineal en z (x). n

>>syms p q n >>y=’z^(1/(1-n))’; >>y=subs(y,’z(x)’,’z’) y =z(x)^(1/(1-n)) >>op=diff(y,’x’)+p*y-q*y^n; >>a=op/’z(x)^(n/(1-n))’ a=(z(x)^(1/(1-n))/(1-n)*diff(z(x),x)/z(x)p*+... z(x)^(1/(1-n))-q*(z(x)^(1/(1-n)))^n)/(z(x)^(n/(1-n)))

on y  + y/x = log(x)y 3 , aplicando el método antes explicado. Ejemplo 2.11 Resolver la ecuaci´

Clase Pr´ actica 2

21

Resoluci´ on del ejemplo 2.11 >>syms x >>y=’z^(1/(1-3))’; >>y=subs(y,’z(x)’,’z’) y =z(x)^(1/(1-3)) >>op=diff(y,x)+y/x-log(x)*y^3 op =-1/2/z(x)^(3/2)*diff(z(x),x)+1/z(x)^(1/2)/x-log(x)/z(x)^(3/2) >>a=op/’z(x)^(3/(1-3))’; >>a=simplify(a) a =-1/2*(diff(z(x),x)*x-2*z(x)+2*log(x)*x)/x

A continuación se resuelve la anterior ED lineal usando DSOLVE >>z=dsolve(’-1/2*Dz+z/x=log(x)’,’x’) z =2*log(x)*x+2*x+x^2*C1

Finalmente se debe deshacer el cambio z = 1/y2 para obtener la solución.

2.6

FACTOR INTEGRANTE

Recordemos que si la ED P (x, y)dx + Q(x, y)dy = 0 no es exacta, siempre existe una función µ = µ(x, y) tal que µP (x, y)dx + µQ(x, y)dy = 0 s´ı es exacta. Los casos aparentemente más sencillos son aquellos en que puede encontrarse un factor µ, llamado “factor integrante” (f.i.), que sólo depende de x o sólo depende de y. Un ejemplo de f.i. que depende de ambas variables x e y es aquél que encontramos para la ED (xy2 yx2 )dx + x3 dy = 0. (2.3)

−

Si llamamos P (x, y) = xy 2

− yx

2

y Q(x, y) = x 3 comprobaremos que 1 1 = 2 2 xP + yQ xy

es un f.i. para la ED (2.3). La comprobación de lo anteriormente dicho podemos hacerla fácilmente con MATLAB. >>syms x y >>P=’x*y^2-y*x^2’;Q=’x^3’; >>Test=diff(P,y)-diff(Q,x); %No es exacta >>M=1/(x*P+y*Q); %f.i. en (x,y) >>P=M*P;Q=M*Q ; >>Test1=simplify(diff(P,y)-diff(Q,x)) Test1=0 La respuesta Test1 =0 demuestra “experimentalmente” que M = x−2 y−2 es ciertamente un f.i.

para (2.3).

Clase Pr´ actica 2

22 El potencial F (x, y) se obtiene mediante las siguientes operaciones. >>F1=int(P,x) >>derf=Q-diff(F1,y); >>f=int(derf,y) >>F=F1+f

para obtener que F = log(x)

− x/y, por lo tanto la solución es log(x) − x/y = C .

on homogénea. El alumno Nota 2.1 El anterior procedimiento puede aplicarse a cualquier ecuaci´ puede comprobar que en el caso P = x2 y 2 , Q = xy, el f.i. que se obtiene depende de una sola variable.

−

Ejemplo 2.12 Usar MATLAB para resolver la siguiente ED del Bolet´ın 2. y  =

√ x + y.

(2.4)

Resoluci´ on del ejemplo 2.12 La solución puede obtenerse mediante el cambio z = x + y, ya que es de la forma: y  = G(ax + by + c) con a = b = 1, c = 0 y G(t) = t. Sin embargo, a´ un no tenemos explicación para el siguiente resultado MATLAB:

√

>>dsolve(’Dy=sqrt(x+y)’,’x’) ans =0

Es simple comprobar que y = 0 no es solución de la ecuación. Sin embargo, si elevamos al cuadrado ambos miembros, aunque introduzcamos soluciones adicionales, entonces s´ı obtenemos respuesta: >>dsolve(’Dy^2=(x+y)’,’x’) ans = [ 1-x] [ (-lambertw(-exp(-1-1/2*x+1/2*C1))-1)^2-x]

Notemos que ahora tenemos todas las soluciones de

±√ x + y. Mediante comprobación directa vemos que y = 1 − x es solamente solución de √ y = − x + y y =

Por otra parte, teniendo en cuenta que lambertw(x)=w(x) cumple w(x)ew(x) = x, obtenemos de la propia respuesta dada por MATLAB (-lambertw(-exp(-1-1/2*x+1/2*C1))-1)^2-x

Clase Pr´ actica 2

23

la siguiente expresión: (1

± √ x + y)exp(−1 ± √ x + y) = exp(−x/2 + C )

Para nuestro problema original, tenemos como respuesta la siguiente ecuación (1 +

√ x + y)exp(−1 − √ x + y) = exp(−x/2 + C )

(2.5)

a la cual podemos llegar con el cambio z = x + y, sin usar herramientas MATLAB. Insistamos en el problema. La respuesta MATLAB se obtuvo al integrar la ecuación y = (y  )2

− x,

que es de la forma, ”resuelta respecto a y”, del tipo ”Lagrange” y = ( 1)x + (y )2 .

−

Hacemos y  = p y derivamos con respecto a x para obtener la ecuación lineal en x = x( p) ( p + 1)x = 2 p, que se integra fácilmente obteniendo x = 2 p y = p 2

− log( p + 1) + C, − 2 p + log( p + 1) − C,

(2.6)

teniendo en cuenta que y = p 2 x. Aunque parezca a primera vista que es dif´ıcil, en realidad es muy fácil eliminar p = x + y de las ecuaciones (2.6), para deducir (2.5). Por lo tanto (2.4) puede ser considerada de dos maneras. A saber, como de la forma y  = G(ax + by + c), y como Lagrange y = xf (y  ) + g(y ). Tambi´ en podemos ver la ED del ejemplo (2.12) como la ecuaci´ on no exacta

−

√

√ x + y dx − dy = 0. Lo más interesante ahora es que resulta infructuoso el intento de hallar un f.i. que sólo dependa de x o que sólo dependa de y. ¡Comprobarlo! En este caso la sugerencia es descubrir que la ED admite un f.i. de la forma µ(x, y) = ϕ(x + y). Para obtenerlo procedemos de la forma siguiente. Tenemos P = x + y y Q = 1. Si µ es de la forma anterior es f.i. entonces debe ser

√

−

∂ (µP ) ∂ (µQ) = , ∂y ∂x es decir,

√

µy x + y + µ

1 √ = −µ , 2 x+y x

Clase Pr´ actica 2

24 o equivalentemente

√

ϕ (x + y) x + y + ϕ(x + y)

1 √ = −ϕ (x + y). 2 x+y

Debido a la forma de µ , y a la expresión anteriormente calculada, si hacemos t = x + y, obtenemos, separando las variables. dϕ dt = (2.7) ϕ 2 t(1 + t)

− √

√

En todo lo anterior, tener en cuenta que hemos usado la regla de la cadena para calcular las derivadas de µ respecto a x e y, esto es, ∂µ ∂µ (x, y) = ϕ (x + y), (x, y) = ϕ (x + y), ∂x ∂y siendo ϕ  la derivada ordinaria de la función ϕ(t). Integrando (2.7) obtenemos que el f.i. de la ED es µ(x, y) =

1 √ . 1+ x+y

√

Notar que en este caso ϕ(t) = 1/(1 + t). Sabiendo que µPdx + µQdy = 0 es exacta procedemos al cálculo del potencial F (x, y). La integral general es F (x, y) = C.

Ejercicio 2.4 Resolver el siguiente PVI (Bolet´ın 2) y  =



y , y(1) = 4. x

´ CLASE PRACTICA 3 Resoluci´ on simb´ olica de sistemas diferenciales

OBJETIVOS Ejercitar el uso combinado de técnicas matemáticas y herramientas simbólicas del MATLAB para el estudio y resoluci´ on de sistemas diferenciales lineales y ecuaciones lineales de orden superior.

3.1

´ SISTEMA LINEAL HOMOGENEO DE COEFICIENTES CONSTANTES

Ejemplo 3.1 Resolver el sistema diferencial lineal siguiente x = 5x + 3y y = 2x 10y

− − −

1. Mediante el comando DSOLVE. 2. Usando herramientas MATLAB pero siguiendo la via matricial.

Resoluci´ on del ejemplo 3.1 I- MATLAB resuelve directamente este tipo de sistemas utilizando el comando DSOLVE con la sintaxis conocida. Si llamamos S a la solución y escribimos >>S=dsolve(’Dx=-5*x+3*y,Dy=-2*x-10*y’)

MATLAB devuelve las componentes simbólicas de la solución: S= x: [1x1 sym] y: [1x1 sym]

25

Clase Pr´ actica 3

26

Para obtener expl´ıcitamente la solución S se pide cada una de sus componentes [S.x,S.y] o pretty([S.x,S.y]) >>pretty([S.x,S.y]) [-2C1 exp(-8t)+3C1 exp(-7t)+3C2 exp(-7t)-3C2 exp(-8t), -2C1 exp(-7t)+2C1 exp(-8t)+3C2 exp(-8t)-2C2 exp(-7t)]

es decir, la respuesta es x = C 1 ( 2e−8t + 3e−7t ) + 3C 2 (e−7t e−8t ) y = 2C 1 ( e−7t + e−8t ) + C 2 (3e−8t 2e−7t )

− −

− −

Alternativamente se puede obtener directamente la respuesta indicando las componentes del vector solución. >>[x,y]=dsolve(’Dx=-5*x+3*y,Dy=-2*x-10*y’) x=-2*C1*exp(-8*t)+3*C1*exp(-7*t)+ 3*C2*exp(- 7*t)-3*C2*exp(-8*t) y=-2*C1*exp(-7*t)+2*C1*exp(-8*t)+ 3*C2*exp(- 8*t)-2*C2*exp(-7*t)

Otra v´ıa para resolver el anterior sistema lineal homogéneo, de coeficientes constantes, es seguir la via matricial directa. Para ello comenzamos escribiendo el sistema en forma matricial:  5 3 x x x = = A 2 10 y y y

   −

− −

   

Conocer la solución general de este sistema significa conocer una matriz fundamental R(t) del mismo. Eso es lo que calcularemos mediante comandos MATLAB. II- Sabemos que R(t) = e tA . Para calcular e tA usaremos el comando EXPM de MATLAB. >>A=[-5 3;-2 -10];%Crear la matriz A >>syms t %declarar simbolica t >>Rt=expm(t*A); %calcular matriz fundamental >>X1=Rt(:,1); %opcional calculo de un SFS, >>X2=Rt(:,2); %es decir, de las columnas de R(t) >>diff(X1)-A*X1 %Pueden hacerse comprobaciones ans = [ 0] [ 0] >>diff(X2)-A*X2 ans = [ 0] [ 0]

Clase Pr´ actica 3

27

Visualizar Rt y comparar con los resultados anteriores >>Rt Rt = [-2*exp(-8*t)+3*exp(-7*t),3*exp(-7*t)-3*exp(-8*t)] [-2*exp(-7*t)+2*exp(-8*t),3*exp(-8*t)-2*exp(-7*t)] >>X1 X1 = [-2*exp(-8*t)+3*exp(-7*t)] [-2*exp(-7*t)+2*exp(-8*t)] >>X2 X2 = [-3*exp(-8*t)+3*exp(-7*t)] [-2*exp(-7*t)+3*exp(-8*t)]

3.2

SISTEMA LINEAL COMPLETO DE COEFICIENTES CONSTANTES

Ejemplo 3.2 Resolver el sistema diferencial lineal siguiente: x y = e t y + 5x + 2y = sen (3 + t),

−

con x(0) = 1, y(0) =

−1.

1. Directamente usando DSOLVE. 2. Utilizando el MVP (1 .

Resoluci´ on del ejemplo 3.2 1. Mediante DSOLVE >>[x,y]=dsolve(’Dx-y=exp(t),Dy+5*x+2*y=sin(3+t)’,’x(0)=1,y(0)=-1’) x= 1/2*exp(-t)*sin(2*t)*(5/8-1/5*sin(3)+1/10*cos(3))+1/2*exp(-t)*... sin(2*t)*(-3/8-1/5*cos(3)-1/10*sin(3))+(5/8-1/5*sin(3)+1/10*... cos(3))*exp(-t)*cos(2*t)+3/8*exp(t)+1/5*sin(3+t)-1/10*cos(3+t) y= 1

Método de Variación de Parámetros

Clase Pr´ actica 3

28 -5/2*exp(-t)*sin(2*t)*(5/8-1/5*sin(3)+1/10*cos(3))-1/2*...

exp(-t)*sin(2*t)*(-3/8-1/5*cos(3)-1/10*sin(3))+1/5*cos(3+t)... +(-3/8-1/5*cos(3)-1/10*sin(3))*exp(-t)*cos(2*t)-5/8*exp(t)... +1/10*sin(3+t)

2. Para la resolución mediante el MVP usaremos la fórmula correspondiente para obtener una soluci´ on particular del sistema diferencial lineal X  = A(t)X + B(t). La soluci´ on general del mismo es t (3.1) X (t) = R(t)V + R(t) R(s)−1 B(s)ds,

 0

siendo R(t) una matriz fundamental y V

∈R . n

La fórmula (3.1) corresponde al caso general A(t) de coeficientes variables. El caso particular de coeficientes constantes aparece a continuación. Recordemos que, en caso de sistemas diferenciales homogéneos con coeficientes constantes, X  = AX , una matriz fundamental está dada por R(t) = e tA . Luego, si se trata de un sistema completo de coeficientes constantes X  = AX + B(t), una soluci´ on particular X p está dada por X p = e

tA



t

e−sA B(s)ds.

0

Para resolver el sistema diferencial lineal del ejemplo 3.2, es decir x y = e t y + 5x + 2y = sen(3 + t),

−

con x(0) = 1, y(0) = 1, usaremos la fórmula general de variación de parámetros para hallar una solución particular X p del sistema no homogéneo.

−

>>syms s t >>A=[0 1;-5 -2]; >>Bt=[exp(t);sin(3+t)]; >>Bs=subs(Bt,s,t); >>Rt=expm(t*A) Rt=[exp(-t)*cos(2*t)+1/2*exp(-t)*sin(2*t),1/2*exp(-t)*sin(2*t)] [-5/2*exp(-t)*sin(2*t), exp(-t)*cos(2*t)-1/2*exp(-t)*sin(2*t)] >>Rs=subs(Rt,s,t) Rs=[exp(-s)*cos(2*s)+1/2*exp(-s)*sin(2*s),1/2*exp(-s)*sin(2*s)]

Clase Pr´ actica 3

29

[-5/2*exp(-s)*sin(2*s),exp(-s)*cos(2*s)-1/2*exp(-s)*sin(2*s)] >>Xp=Rt*int(Rs\Bs,s,0,t); >>Prueba=diff(Xp)-A*Xp-Bt;% Comprobacion

Si >>simplify(Prueba)

no produce [0;0], ejecutar >>ezplot(Prueba(1),[-1,1]), >>ezplot(Prueba(2),[-1,1])

para apreciar, con cierta significaci´ on, que realmente Prueba es cero.(2 La formulación expl´ıcita de la solución particular X p se compone de una cadena de caracteres demasiado larga y por tanto de reducida utilidad, en efecto, si hacemos la conversión de SYM a STRING mediante el comando CHAR, tendremos que >>size(char(Xp)) ans = 1

1669

es decir, ¡Xp está compuesto de 1669 caracteres! La figura 3.1 muestra gráficamente el resultado de ejecutar ezplot(Prueba(2),[-1,1]). El rango de variación en el intervalo [ 1, 1] sugiere la nulidad con una significaci´ on del orden de 10 −15 .

−

Ejemplo 3.3 Resolver el siguiente sistema lineal completo de coeficientes constantes x + x + 2y = cos(t) + sen (t) + e−t y 2x + y = sen (t) cos(t),

−

x(0) = 1, y(0) = 1.

−

Obtenci´ on de la soluci´ on del ejemplo 3.3 mediante DSOLVE >>sis1=’Dx+x+2*y=cos(t)+sin(t)+exp(-t)’; >>sis2=’Dy-2*x+y=sin(t)-cos(t)’; >>[x,y]=dsolve(sis1,sis2,’x(0)=1,y(0)=1’) x =cos(t)-1/2*exp(-t)*sin(2*t) y =1/2*exp(-t)+sin(t)+1/2*exp(-t)*cos(2*t) 2

Los errores menores que eps

≈ 2.2e − 16 se asume que son cero.

Clase Pr´ actica 3

30

Fig. 3.1: Evidencia gráfica de que Prueba(2) debe ser nula en [ 1, 1]

−

3.3

´ DE SISTEMAS DIFERENCIALES RESOLUCION ´ LINEALES HOMOGENEOS DE COEFICIENTES CONSTANTES POR LA VIA ESPECTRAL

Definici´ on 3.1 Valores y vectores propios. Sea A una matriz constante, cuadrada de orden n. Decimos que λ es un valor propio o autovalor de A si existe un vector u no nulo tal que Au = λu. A este u se le llama vector propio o autovalor asociado a λ.

3.3.1

C´ alculo de valores y vectores propios via MATLAB

Si A es una variable que contiene una matriz cuadrada, entonces >>[V,D] = eig(A)

produce la matriz V = [u1 un ] de autovectores, y la matriz diagonal D, tales que AV = V D. Es decir, la diagonal principal de D contiene a los autovalores de A. Si hacemos

···

>>valores_propios=eig(A)

almacenamos en la variable vectorial valores_propios a los valores propios de A.

Clase Pr´ actica 3

3.3.2

31

El m´ etodo espectral

El siguiente resultado nos ofrece la herramienta te´ orica para resolver un sistema diferencial lineal homogéneo, de coeficientes constantes por la via de los valores y vectores propios.

Teorema 3.1 Si la matriz A tiene n valores propios diferentes λ1 ,...,λn , y uk es un vector propio asociado a λk , k = 1,...,n, entonces λ1 t

{e

u1 ,

λn t

··· , e

un ,

}

es un sistema o conjunto fundamental de soluciones (SFS) para el sistema diferencial X  = AX .

Ejemplo 3.4 Resolver por la via espectral, y aplicando herramientas MATLAB, el siguiente problema con valores iniciales X  =

 

1 1 4

2 0 4

−

−1 1 5

 

X (t) X (0) =

−    1 0 0

Indicaciones para resolver el ejemplo 3.4 Usar el comando EIG para calcular los valores y vectores propios de la matriz del sistema y verificar si son diferentes dos a dos para aplicar el Teorema 3.1. Recordar que una vez tengamos una matriz fundamental R(t), entonces c1 X (t) = R(t) c2 c3

     

  −    −   

nos da la solución general del sistema. Luego, la única solución estará dada por la solució n del siguiente sistema 1 c1 0 = X (0) = R(0) c2 0 c3 es decir,

    c1 c2 c3

= R(0)−1

1 0 0

Clase Pr´ actica 3

32

3.4

EEDD LINEALES DE ORDEN SUPERIOR

Las EEDD de orden superior pueden ser tratadas con el comando DSOLVE de una manera muy flexible. Se presentan a continuación diversos ejemplos de EEDD homogéneas y completas de coeficientes constantes para las que se obtiene la solución general as´ı como ciertas soluciones particulares. Se presenta también un ejemplo de coeficientes variables que corresponde a una Euler.

3.4.1

EEDD lineales de orden superior homog´ eneas

Ejemplo 3.5 Resolver la ED y  + 4y + 4y = 0 Resoluci´ on mediante MATLAB del ejemplo 3.5 >>S=dsolve(’D2y+4*Dy+4*y=0’) S =C1*exp(-2*t)+C2*exp(-2*t)*t %(sol. gral.)

Obsérvese que

−2 es una ra´ız doble del polinomio caracter´ıstico

Ejemplo 3.6 Resolver el PVI y  + 4y  + 4y = 0, y(0) = 1, y  (0) = 4. Resoluci´ on mediante MATLAB del ejemplo 3.6 >>S=dsolve(’D2y+4*Dy+4*y=0’,’y(0)=1,Dy(0)=4’) S=exp(-2*t)+6*exp(-2*t)*t %(sol. part.)

Notar que se ha obtenido C 1 = 1 y C 2 = 6. En todos los casos las respuestas son exactas, con fórmulas cerradas.

3.4.2

EEDD lineales de orden superior completas

Ejemplo 3.7 Resolver la ED y  + y + 4y = 1. Resoluci´ on mediante MATLAB del ejemplo 3.7 >>S=dsolve(’D2y+Dy+4*y-1=0’) S=1/4+C1*exp(-1/2*t)*cos(1/2*15^(1/2)*t)+C2*exp(-1/2*t)*... sin(1/2*15^(1/2)*t)

− ± i√ 15)/2.

Notar que el polinomio caracter´ıstico tiene ra´ıces complejas conjugadas ( 1

Ejemplo 3.8 Resolver y  = log(x)/x2 , y(1) = 0, y  (1) = 1, y  (1) = 2.

Clase Pr´ actica 3

33

Resoluci´ on mediante MATLAB del ejemplo 3.8 >>y=dsolve(’D3y=log(x)/x^2’,’y(1)=0’,’Dy(1)=1’,’D2y(1)=2’,’x’) y=-1/2*log(x)^2*x+3/2*x^2-2*x+1/2

´ ¿COMO SOLVENTAR DIFICULTADES RELACIONADAS CON DSOLVE? Si una ED no admite solución exacta en términos del catálogo de funciones elementales MATLAB, entonces DSOLVE no puede ofrecer una solución. Existen otras causas que impiden a DSOLVE llegar a una respuesta exacta pero su consideración queda fuera de nuestras posibilidades. En algunos casos el comando DSOLVE puede producir respuestas poco amigables, es decir, excesivamente largas o expresadas en t´ erminos de funciones especiales poco usuales. En tales casos debemos hallar métodos alternativos que produzcan soluciones, quiz´ as aproximadas pero mucho más potables. Para ilustrar lo dicho, a continuaci´ on se ofrece un ejemplo de EDL homogénea de cuarto orden, cuya soluci´ on se obtiene de forma aproximada, siguiendo un método basado en nuestros conocimientos teóricos. En el apéndice se ampl´ıan los detalles técnicos del método y se ejerce una cr´ıtica necesaria.

Ejemplo 3.9 Considere la EDLH (3 L[y] = 0, donde L[y] = y iv)

− 6y

iii)

+ 9yii)

i)

− 24y − 20y.

Intente resolverla, es decir, halle informaci´ on util ´ sobre la soluci´ on.

Desarrollo del ejemplo 3.9 Tratemos de resolver L[y] = 0 utilizando DSOLVE. >>y=dsolve(’D4y-6*D3y+9*D2y-24*Dy-20*y=0’); >> size(char(y)) ans = 1 2765

¡La respuesta (4 consta de 2765 caracteres! Siendo la EDL de coeficientes constantes, podemos intentar hallar una respuesta, satisfactoria en cierto sentido (ver Apéndice), calculando un SFS a partir del polinomio caracter´ıstico P = m4

3

− 6m

+ 9m2

− 24m − 20.

Calculemos aproximadamente las ra´ıces de P con alguno de los comandos SOLVE y ROOTS. El comando SOLVE resuelve exactamente, en términos de radicales, ecuaciones polinómicas hasta el grado 4. ROOTS calcula aproximadamente las ra´ıces de cualquier polinomio mediante un 3 4

Ecuaci´ on Diferencial Lineal y Homogénea La respuesta puede cambiar según la versión MATLAB o el software que se utilice

Clase Pr´ actica actica 3

34

procedimiento muy eficiente y preciso. No obstante, los resultados num´ ericos ericos obtenidos mediante SOLVE SOLVE deber´ deber´ıan ser superiores, en algunos algunos casos, casos, a los producidos producidos por ROOTS. ROOTS. El nivel nivel de coin− 15 cidencia entre ambos es del orden 10 en el caso que nos ocupa, según un se aprecia en las pruebas numéricas ericas que se muestran a continuación. on. >>form >>format at long long e >>S=numeric(solve(’x^4-6*x^3+9*x^2-24*x-20’)) S = 5.291305619991860e+00 6.653067648628088e-01+2.373811784418384e+00i 6.653067648628088e-01-2.373811784418384e+00i -6.219191497174775e-01 >>R=ro >>R=roots ots([1 ([1 -6 9 -24 -20]) -20]) R = 5.291305619991865e+00 6.653067648628106e-01+2.373811784418383e+00i 6.653067648628106e-01-2.373811784418383e+00i -6.219191497174774e-01 >>max(abs(R-S)) ans= 5.329070518200751 5.329070518200751e-015 e-015

El correspondiente SFS es = e 0.6653067648628088 t sen(2. sen(2.373811784418384 t), y1 = e = e 0.6653067648628088 t cos(2. cos(2.373811784418384 t), y2 = e 5.291305619991860 t = e y3 = e = e −0.6219191497174775 t y4 = e ´ Ver en el AP ENDICE cómo omo intentamos valorar la calidad de esta solución on aproximada.

Ejercicio 3.1 Dadas las siguientes EEDD iv ) a.- y iv) iii) b.- y iii) iv ) c.- y iv)

10y − 10y − 2y − 6y

iii) iii)

ii) + 3y 3yii)

i)

20y = 0, − 2y − 20y

ii) ii)

+ 3y 3yi) = 0, 0,

iii) iii)

ii) + 9y 9yii)

24y − 24y

i)

+ 50y 50y = 1,

obtener una soluci´ on particular para la ED del apartado c) mediante el método etodo de coeficientes coeficientes indeterminados (MCI) sin ayuda del MATLAB, y un SFS para a), b) y c).


35

´ DE EULER Y DE EULER ECUACION GENERALIZADA

3.5

La ecuaci´ ecuación on de Euler de segundo orden es de la forma αx2 y (x) + βxy  (x) + γy( γy (x) = δ (x) Recordemos que ésta se transforma en una ED de coeficiente co eficientess constantes mediante el cambio x cambio x = = e e t , y que la ecuación on de Euler generalizada de segundo orden es de la forma α(ax + b)2 y (x) + β (ax + b)y (x) + γy( γy (x) = δ (x)

Ejemplo 3.10 Resolver (2x (2x + 1)2 y

− (4x (4x + 2)y 2)y − 12y 12y = 8x.

Resoluci´ on del ejemplo 3.10 mediante MATLAB on Aplicando DSOLVE se obtiene que >>y=dsolve(’(2*x+1)^2*D2y-(4*x+2)*Dy-12*y=8*x’,’x’) y=2/3*x^3*(2+x)/(2 y=2/3*x^3*(2+x)/(2*x+1)+ *x+1)+ C1*(2*x+1)^3+C2/( C1*(2*x+1)^3+C2/(2*x+1) 2*x+1)

ales de las siguientes siguientes EEDD, despu´ después es de ser transforma transformadas das en EDL de coeoeEjercicio 3.2 ¿A cu´ ficientes constantes, puede aplicarse el MCI? ¿Cu´ al es la alternativa al MCI? Intentar resolverlas mediante DSOLVE y mediante métodos etodos artesanales. Contrastar los resultados. 1. (2x (2x + 1)2 y (4x (4x + 2)y 2)y 12y 12y = 8x

− − 2. x y + 5xy 5 xy − 2y = −5log(x 5log(x) 3. (10x (10x + 50) y − 3(x 3(x + 5)y 5)y + 7y 7 y = x = x − 4. x y − 6x y  + 2xy 2 xy = = x x + 1 5. 5x y − 6xy + y = e = e 2

2

3

2

2

3.6

2

x

´ CALCULO DE WRONSKIANOS Y APLICACIONES

En lo que sigue se utilizará el wronskiano para el estudio de la independencia lineal de sistemas de funciones. funcion es. También en se aplicar´ aplica r´ a en la obtención on de una EDLH conocidas algunas de sus soluciones. El wronskiano wron skiano juega un u n papel pap el relevante en la l a teor´ıa ıa y pr´ p ráctica actica de las EDL. Si y Si y1 e y e y 2 son funciones definidas sobre el intervalo I intervalo I ,, entonces entonces el wronskiano wronskiano W estas se define como: W ((y1 , y2 ) de éstas )(t) = det W ( W (y1 , y2 )(t



y1 (t) y2 (t) y1 (t) y2 (t)



, t

∈ I.


36

3.6.1 3.6.1

Wronski ronskiano ano e indepe independ nden encia cia linea lineall

Teorema 3.2 Sean x1 (t) y x2 (t) dos soluciones de una ED del tipo + p((t)x (t) + q (t)x(t) = 0, x (t) + p con p(t) y q (t) continuas en I . Entonces se cumple que: x1 y x2 son linealmente independientes en I )(t) = 0 t I W ( W (x1 , x2 )(t )(t) = 0. 0. t I I tal que W ( W (x1 , x2 )(t

∃∈

 ∀ ∈ ⇔ ⇔

⇔ ⇔



Nota 3.1 Sean x1 y x2 de clase 1 en I , no necesariamente soluciones de una EDLH. Entonces 1) W ( )(t) = 0 para alg´ un t I W (x1 , x2 )(t x1 y x2 son linealmente independientes en I . 2) Sin embargo, es posible que x1 y x2 sean linealmente independientes en I y W ( )(t) = 0 W (x1 , x2 )(t 2 2 en alg´ un I . Por ejemplo, t y t son l.i. en [ 1/2, 1/2] y W ( 0. W (t, t )(0) = 0.



∈ ⇒ ⇒

−

senh (t) Ejemplo 3.11 Calcular el Wronskiano de e−2t y senh

− 2 cosh(2 cosh(2tt).

Resoluci´ on on del ejemplo 3.11 >>a=’[exp(-2*t), >>a=’[exp(-2*t), sinh(t)-2*cosh(2*t sinh(t)-2*cosh(2*t)] )] ’; >> b=diff(a b=diff(a,t) ,t) % tambien tambien diff(a) diff(a) b=[-2*exp(-2*t),cosh(t)-4*sinh(2*t)] >> c=[a;b] c=[a;b] c=[exp(-2*t),sinh(t)-2*cosh(2*t)] [-2*exp(-2*t),cosh(t)-4*sinh(2*t)] >>d=det( >>d=det(c) c) %(es el Wronskia Wronskiano no pedido) pedido) d=exp(-2*t)*cosh(t)-4*exp(-2*t)*sinh(2*t)+... 2*exp(-2*t)*sinh(t)-4*exp(-2*t)*cosh(2*t)

Para abreviar puede utilizarse el programa WRONSKY.M, que puede descargarse desde la Web de la asignatura [15].

Ejemplo 3.12 Calcular el Wronskiano de ex , e2x y e3x . Resoluci´ on on del ejemplo 3.12. >>a=’[exp(x),exp(2*x),exp(3*x)]’; >>b=diff(a,’x’); >>c=diff(b,’x’); >>d=[a;b;c]

Clase Pr´ actica 3

37

d=[exp(x),exp(2*x),exp(3*x)] [exp(x),2*exp(2*x),3*exp(3*x)] [exp(x),4*exp(2*x),9*exp(3*x)] >>W=det(d)%%%%(es el Wronskiano) W=2*exp(x)*exp(2*x)*exp(3*x)

Ejemplo 3.13 Determinar si las siguientes funciones son linealmente independientes: f (t) = e t sen (2t) y g(t) = sen (t)cos(t) en el intervalo [1/2, 3/4]. Resoluci´ on del ejemplo 3.13 usando MATLAB >>syms a b c d t >>a=’[exp(t)*sin(2*t),sin(t)*cos(t)]’; >>b=diff(a,t); >>c=[a;b]; >>d=det(c); >>d=simplify(d) d=2*exp(t)*sin(2*t)*cos(t)^2-exp(t)*... sin(2*t)sin(t)*cos(t)*exp(t)*... sin(2*t)-2*sin(t)*cos(t)*exp(t)*cos(2*t) >>W=subs(d,3/4,t) W=-1.05320353709056

Por lo tanto f y g son linealmente independientes.

Ejemplo 3.14 Estudiar la independencia lineal de f (t) = 1, g(t) = ( 1, 2).

−t + 1, h(t) = 3t

−

Resoluci´ on del ejemplo 3.14 usando MATLAB >>syms t >>a=’[1,-t+1,3*t^2+2*t-4]’; >>b=diff(a,t); >>c=diff(b,t); >>d=[a;b;c]; >>e=det(d) e =-6

Por lo tanto las funciones son l.i.

Ejemplo 3.15 Dadas las familias de funciones

2

+ 2t

− 4 en

Clase Pr´ actica 3

38 a) y1 = e x , y2 = e 2x , y3 = e 3x b) y1 = 1, y2 = x, y3 = 1/(1 + x2 ) c) y1 = sen (x), y2 = x d) y1 = e x , y2 = cos(3x) e) y1 = x, y2 = x 2 + 1, y3 = x 2 + 2x + 1

hallar para cada una la EDL homog´ enea con el menor orden posible para la cual las funciones yi forman parte de un conjunto fundamental de soluciones, considerando por separado los casos

{ }

1. coeficientes constantes (cuando sea posible), 2. coeficientes variables.

Construcci´ on de una EDLH de coeficientes variables teniendo a las funciones del ejemplo 3.15, apartado b) como SFS Pueden ejecutarse comandos MATLAB en l´ınea tal como se muestra a continuaci´ on, o utilizar el programa WRONSKY. >>syms x y dy d2y d3y >>A=[1 x 1/(1+x^2)]; >>A1=[y,A]; >>A2=[dy,diff(A)]; >>A3=[d2y,diff(A,2)]; >>A4=[d3y,diff(A,3)]; >>W=det([A1;A2;A3;A4]) W=-2*(12*d2y*x^3-12*d2y*x+2*d3y*x^2+3*d3y*x^4-d3y)/(1+x^2)^4

La ecuación es W = 0, siendo dy = y  , d2y = y  , etc. Si se selecciona un intervalo dominio I 0 , entonces I 0 no debe contener ceros del wronskiano.

3.7 3.7.1

SOLUCIONES PARTICULARES DE EEDD LINEALES COMPLETAS M´ etodo de coeficientes indeterminados (MCI)

El MCI es aplicable cuando la ED tiene coeficientes constantes y la parte completa es de tipo CI. on particular de la ED Ejemplo 3.16 Obtener una soluci´ x + 3x

4t

2t

− 4x = 3 sen (t) + 7e− − e

Clase Pr´ actica 3

39

Resoluci´ on del ejemplo 3.16 con DSOLVE >>y=dsolve(’D2x+3*Dx-4*x=0’) y=C1*exp(t)+C2*exp(-4*t)

que es la solución general de la homogénea. Teniendo en cuenta la forma de la parte completa de la ED, se busca una solución particular x p que sea combinación lineal de x1 = sen(t), x2 = cos(t), x3 = te −4t , y x4 = e 2t . >>syms t A B C D >>xp=A*sin(t)+B*cos(t)+C*t*exp(-4*t)+D*exp(2*t);

que es la solución particular buscada. Para determinar los coeficientes A, B , C y D se debe sustituir x p en la ED. >>l=diff(xp,2,t)+3*diff(xp,t)-4*xp l= -5*A*sin(t)-5*B*cos(t)-5*C*exp(-4*t)+6*D*exp(2*t)+3*A*cos(t)-3*B*sin(t) >>eqn=l-3*sin(t)-7*exp(-4*t)+exp(2*t) eqn= -5*A*sin(t)-5*B*cos(t)-5*C*exp(-4*t)+6*D*exp(2*t)+3*A*cos(t)-3*B*sin(t)-... 3*sin(t)-7*exp(-4*t)+exp(2*t)

Para resolver eqn=0 y as´ı obtener A, B, C , y D, se evalúa eqn y sus derivadas 1ra, 2da y 3ra en un punto (en este caso tomamos el punto 0), obteniéndose un sistema cuya resolución nos da A, B , C , y D. >>eqn1=subs(eqn,0,t) eqn1= -6-5*B-5*C+6*D+3*A >>eqn2=subs(diff(eqn,t),0,t) eqn2= -5*A+27+20*C+12*D-3*B >>eqn3=subs(diff(eqn,t,2),0,t) eqn3= -108+5*B-80*C+24*D-3*A >>eqn4=subs(diff(eqn,t,3),0,t) eqn4= 5*A+459+320*C+48*D+3*B

Clase Pr´ actica 3

40

A continuación se resuelve el sistema lineal formado por las 4 ecuaciones eqn1, eqn2, eqn3 y ognitas A, B, C , y D. eqn4, con las inc´ >>[A B C D]=solve(eqn1,eqn2,eqn3,eqn4) A =-15/34 B =-9/34 C =-7/5 D =-1/6

Una solución particular es xp=-15/34 *sin(t)-9/34 *cos(t)-7/5*t*exp(-4*t)-1/6 *exp(2*t)

El procedimiento anterior usado para el cálculo de los coeficientes no es aconsejable sin la ayuda de instrumentos de cálculo, pero s´ı parece ser más adecuado para simular el proceso en el ordenador. En su defecto puede usarse el programa MCI descargable desde la Web de la asignatura [15].

3.7.2

M´ etodo de variaci´ on de par´ ametros (MVP)

on particular de la ED Ejemplo 3.17 Obtener una soluci´ t2 x + tx

− x =

sen (t).

Se utilizar´ a que x1 = t y x2 = 1/t son soluciones l.i. de la ED homogénea. (Se pueden obtener con DSOLVE)

Resoluci´ on del ejemplo 3.17 La solución general de la homogénea es

∈ R

xh = At + B/t, A, B y se busca una particular de la completa de la forma

x p = A(t)t + B(t)/t. Sabemos que A  (t) y B  (t) son soluciones del sistema



t 1

1/t 1/t2

−

Procedemos a resolver dicho sistema.

   A (t) B  (t)

=

0 sen(t)/t2



Clase Pr´ actica 3

41

>>syms t >>a=[t,1/t] >>b=diff(a); >>c=[a;b] c = [t,1/t ] [1,-1/t^2] >>W=det(c) W =-2/t >>A=int(sin(t)/(2*t^2)) >>B=int(-sin(t)/2) >>xp=A*t+B*1/t xp=(-1/2*sin(t)/t+1/2*cosint(t))*t+1/2*cos(t)/t >>pretty(sym(xp)) / sin(t) \ cos(t) |- 1/2 ------ + 1/2 cosint(t)| t + 1/2 -----\ t / t

on particular de la ED x Ejemplo 3.18 Obtener una soluci´

−

e3t  6x + 9x = 2 . t

Resoluci´ on del ejemplo 3.18 >>y=dsolve(’D2x-6*Dx+9*x=0’) y=C1*exp(3*t)+C2*exp(3*t)*t %(sol. gral. homogenea)

Para obtener xp=A(t)*exp(3*t)+B(t)*exp(3*t)*t,

nos basamos en que A  (t) y B  (t) son soluciones del sistema



e3t te3t (e3t ) (te3t )

  

Por lo tanto, se procede a resolver dicho sistema >>syms t >>a=[exp(3*t),t*exp(3*t)] >>b=diff(a,t); >>c=[a;b] c=[exp(3*t), t*exp(3*t)] [3*exp(3*t),3*exp(3*t)*t+exp(3*t)] >>W=det(c) W =exp(3*t)^2

A (t) B  (t)

=

0 3t 2 e /t



Clase Pr´ actica 3

42 >>A=int(-1/t) A =-log(t) >>B=int(1/t^2) B =-1/t >>xp=A*exp(3*t)+B*t*exp(3*t) xp=-log(t)*exp(3*t)-exp(3*t)

´ Esta es la solución particular obtenida.

Ejemplo 3.19 Resolver mediante el MVP la ED x + 3x + 2x = e 3t cos(5t) + 2t4 .

Resoluci´ on del ejemplo 3.19 >>dsolve(’D2x+3*Dx+2*x=0’) % hallar un SFS ans= C1*exp(-2*t)+C2*exp(-t)

Se buscan A(t) y B(t) para tener una solución particular x p cuya expresión pretendemos que sea x p = A(t)e−2t + B(t)e−t . El procedimiento interactivo (l´ınea de comandos) aparece a continuaci´ on. Puede usarse el programa MVP.M que se descarga desde la Web de la asignatura [15]. >>syms t >>a=[exp(-2*t),exp(-t)];b=diff(a,t);c=[a;b] c =[exp(-2*t),exp(-t)] [-2*exp(-2*t),-exp(-t)] >>H=[0;exp(3*t)*cos(5*t)+2*t^4]; >>S=inv(c)*H;A=int(S(1));B=int(S(2)); >>xp=simplify(A*exp(-2*t) +B*exp(-t)) xp=-8/205*cos(t)^5*exp(3*t)+72/205*cos(t)^4*sin(t)... *exp(3*t)+2/41*cos(t)^3*exp(3*t)-54/205*cos(t)^2... *sin(t)*exp(3*t)-1/82*cos(t)*exp(3*t)+9/410*... sin(t)*exp(3*t)+t^4-6*t^3+21*t^2-45*t+93/2

Como alternativa podemos usar el MCI y el principio de superposición, pero en cualquier caso es mucho más cómodo usar directamente DSOLVE que aplica automáticamente el principio de superposición.

Clase Pr´ actica 3

3.8

43

´ DE ORDEN REDUCCION

on Ejemplo 3.20 Resolver la ED x + (2/t)x + x = 0, usando que f (t) = sen (t)/t es una soluci´ particular.

Resoluci´ on del ejemplo 3.20 Buscamos una soluci´ on de la forma x = sen(t)/t v(t). Derivando y sustituyendo en la ED y haciendo v  = y se llega a una ED de primer orden en y.

∗

y + (2/t + 2f  (t)/f (t))y = 0. Resolviendo esta ED se obtiene su solución general. y = F (t, C 1 ). Como y = v  , integrando se obtiene v = G(t, C 1 , C 2 ). Finalmente haciendo x = f v

∗

se obtiene la solución general de la ED. Los pasos con MATLAB se indican a continuación. >>syms t C2 %ecuacion original x’’+2/tx’+x=0 >>f=sin(t)/t >>diff(f) ans=cos(t)/t-sin(t)/t^2 >>diff(f)/f ans=(cos(t)/t-sin(t)/t^2)/sin(t)*t

Se resuelve la ED y  + (2/t + 2f  /f )y = 0 copiando f’/f =diff(f)/f en DSOLVE >>y=dsolve(’Dy+2*(1/t+((cos(t)/t-sin(t)/... t^2)/sin(t)*t))*y=0’) y =-2*C1/(-1+cos(2*t))

Como y = v  , se integra >>v=int(y,t)+C2

Clase Pr´ actica 3

44 v =-C1/tan(t)+C2

Se usa que x = f v

∗

>> x=sin(t)/t*v x=sin(t)/t*(-C1/tan(t)+C2) >>simplify(x) ans= -(C1*cos(t)-C2*sin(t))/t %(sol. gral.)

3.9

EJERCICIOS COMPLEMENTARIOS

Ejercicio 3.3 1. Investigar si las funciones y 1 (t) = e t , y 2 (t) = tet , e y 3 (t) = t 2 et , son linealmente independientes en [a, b]. 2. Resolver el siguiente problema de CI 2y  + 5y + 5y = 0, y(0) = 0, y  (0) = 1/2. 3. Resolver 9yiv)

− 6y

iii)

+ 46yii)

− 6y

i)

+ 37y = 0.

4. Hallar la soluci´ on general del sistema de ecuaciones y + z  = exp(x) y + z  = x 5. Determinar la soluci´ on del sistema x = y y = z z  = y

−

x(0) = 1, y(0) = 1, z (0) = 0. 6. Calcular una matriz fundamental para el sistema del apartado 5). 7. Intentar hallar la soluci´ on de las siguientes EEDD mediante MATLAB.

− 4ty + y = 1 − 4t y + y = 0 (c) y  − 4ty  + ty = 0

(a) y  (b) y 

2

Clase Pr´ actica 3

45

(d) ty 4y + y = 0 (e) x + 5x 50x = 1/t

−

−

´ Ver en el APENDICE una notas relativas a las funciones especiales que pueden aparecer en la solución de algunos de los problemas anteriores.

3.10

´ APENDICE

3.10.1

Notas sobre funciones especiales

Para complementar lo dicho en la clase No. 1 -pag. 4 de este material- sobre la función Lambertw, mencionaremos otras funciones básicas también incorporadas a MATLAB. Estas son: WhittakerW, WhittakerM , besselj, bessely, Ei, dilog, la función de error erf , Heaviside, Dirac, etc, que tambi´ en pueden ser conocidas mediante MHELP. Se trata de funciones especiales para el matemático y el f´ısico debido a que son la soluci´ on de problemas espec´ıficos, frecuentes y de gran importancia en determinados ámbitos de la matemática, la ingenier´ıa o la f´ısica. La filosof´ıa que sigue el MATLAB respecto a estas funciones es la usual. Expresar las soluciones de una ED en términos de ciertas funciones especiales que hemos adoptado, y que son bien conocidas en el sentido de que tenemos suficiente información acerca de ellas: intervalos de convexidad y de crecimiento, máximos y m´ınimos locales, comportamiento asint´ otico en el infinito, derivabilidad, singularidades, desarrollo en serie de potencias, etc. Las EEDD de coeficientes variables del último listado de ejercicios serán estudiadas en el 2do cuatrimestre utilizando un método relativamente simple, basado en asumir a priori que la solución admite un desarrollo en serie de potencias.(5 Con estos apuntes no se pretende que el alumno estudie detalladamente todas estas funciones sino que aumente su cultura profesional, conociendo aún mejor cuál es la tarea que desempeñan las funciones especiales y que éstas no se reducen a las conocidas exponencial, seno, coseno, etc.

3.10.2

Sobre el m´ etodo basado en el c´ alculo aproximado de las ra´ıces de la ecuaci´ on caracter´ıstica

Para resolver la ecuación lineal homogénea de coeficientes constantes ayiv) + byiii) + cyii) + dyi) + ey = 0 aplicamos el comando DSOLVE. Si la respuesta que se obtiene con DSOLVE no es satisfactoria, una variante poco ortodoxa, que no recomendamos, consiste en tomar el polinomio caracter´ıstico P = as4 + bs3 + cs2 + ds + e 5

Notar que ninguna de las EEDD del problema (7) del ejercicio 3.3 es del tipo Euler.

Clase Pr´ actica 3

46

y calcular aproximadamente sus ra´ıces mediante el comando SOLVE tal como se muestra a continuación >>format long e % se sugiere formato largo >>S=numeric(solve(’a*x^4+b*x^3+c*x^2+d*x+e’))

Finalmente damos como solución el correspondiente SFS. ¿Es exacta la respuesta? Obviamente no lo es. En las siguientes l´ıneas se reconsidera el ejemplo 3.9 de la página 33, para ilustrar el nivel de ajuste que se obtiene al calcular una solución aproximada y˜ por este método. Lo primero a tener en cuenta es que el criterio que aqu´ı se sigue es el de medir el error en la ecuaci´ on: L[˜ on: y˜ y . y ] , en lugar del error en la soluci´ Tomemos por caso a L[˜ y] en los intervalos [ 2, 2] y [20, 60].

|

|

| |

−

 − 

>>syms t; >>y1=exp(t*0.6653067648628088)*sin(2.373811784418384*t);

y calculemos Z = L[y1 ]. >>Z=diff(y1,4)-6*diff(y1,3)+9*diff(y1,2)-24*diff(y1,1)-20*y1;

Hagamos (6 >>subplot(1,2,1);ezplot(Z,[-2,2]),title(’’) >>subplot(1,2,2);ezplot(Z,[20,60]),title(’’)

para obtener los gráficos que se muestran en la figura 3.2.

¿Tienen calidad los resultados? El m´ etodo alternativo a DSOLVE que hemos utilizado es aproximado, y concuerda con nuestros conocimientos teóricos. Sin embargo, en el a´mbito experimental los baremos cambian. La calidad de los resultados obtenidos se mide, en este caso, a partir de la estimación que hagamos del residuo Z = L[˜ y(t)] sobre ciertos intervalos. Hemos comprobado que Z es del orden de 10 −15 10 = 10−14 cuando y˜ está dada por y˜(t) = e α t sen(β t), 2 t 2

| | |

|

| | − ≤ ≤

×

siendo α = 0.6653067648628088, y β = 2.373811784418384. ¿Podemos aceptar que Z = 0 en [ 2, 2]? La respuesta podr´ıa ser afirmativa si tenemos en cuenta que los c´ alculos MATLAB se hacen con (7 una precisión finita. Parece correcto aceptar que los errores experimentales son nulos cuando éstos son < eps 2.22e 016, y éste no es el caso que nos ocupa.(8 ¿Pero acaso 10−14 no es suficientemente pequeño? ¿Podemos aceptar como SFS al obtenido por el método descrito anteriormente?

−

≈

6

−

Agregando title(”) evitamos que en el t´ıtulo de la gráfica aparezcan cadenas interminables de caracteres. Se llama precisión al número de d´ıgitos que forman la mantisa. 8 Por eps se representa al epsilon de máquina, es decir, el número positivo más peque˜ no, representable internamente, tal que 1 + eps > 1. Hacer >> eps para obtener su valor con 16 d´ıgitos. 7

Clase Pr´ actica 3

47

iv)

Fig. 3.2: Gráficas de Z = y 1

− 6y

iii) 1

ii)

+ 9y1

i) 1

− 24y − 20y

1

en [ 2, 2] y [20, 60]

−

En principio la respuesta a la primera pregunta es afirmativa. Podemos aceptar un resultado aproximado siempre que tengamos una medida de la calidad de dicha aproximación y que ésta, a su vez, se corresponda con la tolerancia de nuestros instrumentos. En este caso tendremos presente que no tenemos información acerca del error en la solución

y˜ − y, y que nuestra conclusión dependerá del intervalo considerado. La respuesta a la segunda pregunta es negativa, y es que en realidad hemos avanzado muy poco pues la solución general de la ED es de la forma (9 y˜G = C 1 eα3 t cos(β t) + C 2 eα3 t sen(β t) + C 3 eα1 t + C 4 eα2 t ,

(3.2)

de modo que las constantes C 1 ,...,C 4 pueden modificar sustancialmente la escala y por ello podr´ıa ser a´ un más dif´ıcil aceptar que el residuo relativo a y˜G sea cero. A pesar de lo anterior, el criterio basado en medir directamente las discrepancias en la ecuación es una estrategia universalmente aceptada. El llamado “m´ etodo de los residuos (ponderados)”(10 consiste en aceptar como solución aquella que anula al residuo en un sentido más débil que el aqu´ı utilizado. 9 10

Un serio inconveniente es que la función y˜G depende no linealmente de los parámetros α1 ,..., α3 , y β . Est´ an dise˜ nados para resolver problemas de contorno (Tema 7). Ver 20.3 de [5]

§

48

Clase Pr´ actica 3

Para terminar de momento esta discusión, reiteremos que es posible para EEDD lineales y homogéneas de orden cuatro o menor, calcular exactamente las ra´ıces del polinomio caracter´ıstico utilizando SOLVE sin aplicar NUMERIC, solo que la expresión literal de dichas ra´ıces puede ser un STRING con demasiados caracteres. En cualquier caso, las ventajas que se derivan del trabajo con objetos simbólicos no parecen ser concluyentes, al menos con este enfoque. Para resolver éste y muchos otros problemas nos apuntamos al uso de los métodos num´ ericos tradicionales cuyo estudio será abordado en las clases prácticas 5, 6 y 7.

´ CLASE PRACTICA 4 Transformada de Laplace

OBJETIVOS Ejercitar el uso combinado de técnicas matemáticas y herramientas simbólicas del MATLAB para la resolución, mediante la transformada de Laplace, de problemas de valores iniciales asociados a sistemas diferenciales lineales y ecuaciones lineales de orden superior.

´ Y CALCULO ´ DEFINICION

4.1 4.1.1

Definici´ on de la transformada de Laplace

La transformada de Laplace de una función f (t) definida en [0, + ) es una nueva funció n que ˜ o (f )(s), definida por la integral denotaremos por F (s) o f (s)

∞

L

˜ = f (s)



+∞

e−st f (t)dt.

(4.1)

0

˜ está formado por los valores de s para los cuales la integral (4.1) existe. El dominio de f (s) Recordar que si g es integrable en cada intervalo [0, b], entonces



+∞

0



b

g(t)dt = lim b→+∞

g(t)dt.

(4.2)

0

Decimos que g es integrable Riemann en [0, + ] si el l´ımite de la derecha en (4.2) existe y es finito. El s´ımbolo de la izquierda recibe el nombre de integral impropia de primera especie. Si g es integrable Riemann en [0, + ] decimos que g es absolutamente integrable. En el resto de esta clase sólo trabajaremos con funciones f (t) tales que e −st f (t) es absolutamente integrable. ˜ Atendiendo a los objetivos de la clase, en la siguiente subsección abordaremos el cálculo de f (s) (y de la transformada inversa) utilizando comandos MATLAB-MAPLE.

∞

| |

∞

49

Clase Pr´ actica 4

50

4.1.2

C´ alculo mediante comandos MATLAB. Ejemplos

Para el cálculo automatizado de (4.1) MATLAB ofrece el comando LAPLACE de tipo simbólico, cuya sintaxis es >>F=laplace(f)

donde f es una función escalar de la variable simbólica t, previamente declarada y F es una función cuya variable es por defecto s. También podemos formular >> syms u v >>F=laplace(u^2,v) F = 2/v^3

lo que permite elegir las variables a utilizar. on f (t) = 1. Ejemplo 4.1 Calcular la transformada de Laplace de la funci´

Resoluci´ on del ejemplo 4.1. >>syms t s >>laplace(1,t,s) ans =1/s

Ejemplo 4.2 Calcular la transformada de Laplace de la funci´ on f (t) = e −at y f (t) = e −at sen (bt) . Resoluci´ on del ejemplo 4.2 >>syms a; >>laplace(exp(-a*t),t,s) ans =1/(s+a)

Alternativamente se puede crear la función como una cadena de caracteres (string) >>syms t s >>f=’exp(-a*t)’ >>laplace(f,t,s) ans =1/(s+a)

De esta forma no tenemos que declarar previamente a ningún parámetro como simbólico. Invocando al n´ ucleo MAPLE podemos trabajar exclusivamente con STRINGS. >>maple(’f:=t->exp(-a*t)*sin(b*t)’) >>F=maple(’laplace(f(x),x,s)’) F=b/((s+a)^2+b^2)

Clase Pr´ actica 4

51

Ejercicio 4.1 Calcular, utilizando MATLAB, la transformada de Laplace de las funciones f (t) = t, f (t) = t 2 y f (t) = t 3 . Nota para el ejercicio 4.1 Recordar que para n



∈ N la fórmula general viene dada por la expresión:

tn (s) = y para r >

n! sn+1

,

−1, r ∈ R, se tiene, aún más general, que

 

tr (s) =

Γ(r + 1) , sr+1

donde Γ(r) se define a continuación. La función Gamma se define como Γ(t) =

+∞

ut−1 e−u du, t > 0.

0

MATLAB posee el comando gamma que permite simular a Γ(t).

Ejemplo 4.3 (a) Comprobar experimentalmente la igualdad Γ(n) = (n 1)!, utilizando el comando factorial. (b) Calcular Γ(1/2) y Γ(17/3) por las diferentes v´ıas que se sugieren a continuaci´ on.

−

Resoluci´ on del ejemplo 4.3. >>gamma(1/2) >>maple(’gamma(1/2)’) >>numeric(maple(’gamma(1/2)’)) >>int(’t^(1/2-1)*exp(-t)’,0,inf)

Representación gr´ afica de la funci´ on Γ(t) La figura 4.1 sugiere que las rectas t = n, n = 0, 1, 2, 3 son as´ıntotas verticales de la funci´ on Γ(t). Verdaderamente, esto es cierto para cualquier recta t = n, n = 0, 1, 2,....

− − −

4.1.3

− −

Propiedades de la Transformada de Laplace

El ejemplo 4.2 anterior se corresponde con la Propiedad de Traslación: at f (t)(s) = f (s ˜ e

− a),

que MATLAB aplica automáticamente.

Ejercicio 4.2 Calcular la transformada de Laplace de las funciones f (t) = e −at cos(bt) y f (t) = e at t4 utilizando recursos MATLAB-MAPLE de modo que no intervengan variables SYM.

Clase Pr´ actica 4

52

Fig. 4.1: Representación gr´ afica de γ (t) en [ 4, 4]

−

˜ los siguientes casos y reconocer las propiedades que se han aplicado en Ejercicio 4.3 Calcular f en cada caso. (a) f (t) = t 1/3 (b) f (t) = e −2t t5/6 (c) f (t) = g  (t) (d) f (t) = sen (t + a) (e) f (t) = sen (t 1)u(t 1) (f ) f (t) = t 2 senh (t)

−

−

En (c) hacer >>laplace(diff(sym(’g(t)’))). Notar que esta sintaxis no exige declaración previa SYMS. Aplicar la misma técnica al apartado (e), es decir, laplace(sym(’Heaviside(t-1)*sin(t-1)’)), donde Heaviside(t-a) es el nombre codificado de la función de salto unitario u(t a).

−

4.2

TRANSFORMADA INVERSA DE LAPLACE

El comando ILAPLACE ha sido dise˜ nado para calcular exactamente, cuando ello sea posible, la transformada inversa de Laplace. La sintaxis es como sigue. >>ilaplace(expresion_simbolica,s,t)

o simplemente >>ilaplace(expresion_simbolica)

Notar que ILAPLACE calcula la transformada inversa de Laplace de una expresión F(s), produciendo otra expresión f(t). La versión MAPLE es la que sigue >>maple(’invlaplace(F(s),s,t)’)

Clase Pr´ actica 4

53

Ejemplo 4.4 Obtener la transformada inversa de la funci´ on F (s) =

s+3 s2

Resoluci´ on del ejemplo 4.4. >>syms s >>F=(s+3)/s^2; >>ilaplace(F) ans =3*t+1

Alternativamente con MAPLE >>maple(’invlaplace((s+3)/s^2,s,t)’)

Ejemplo 4.5 Calcular las transformadas inversas de Laplace de las siguientes funciones: F (s) =

1 1 1 = = , F (s) , F (s) . s3 + 1 s4 + 1 s5 + 1

Resoluci´ on del ejemplo 4.5. >>F=1/(s^3+1) >>ilaplace(F) ans = 1/3*exp(-t)- 1/3*exp(1/2*t)*cos(1/2*3^(1/2)*t)+1/3*3^(1/2)*exp(1/2*t)*... sin(1/2*3^(1/2)*t)

También puede utilizarse >>F=1/(s^4+1) >>pretty(ilaplace(F)) ans= 1/2*2^(1/2)*(sin(1/2*2^(1/2)*t)*cosh(1/2*2^(1/2)*t)-cos(1/2*2^(1/2)*t)*... sinh(1/2*2^(1/2)*t))

El caso 1/(s5 + 1) MATLAB lo resuelve en el campo complejo, pues factoriza a s5 + 1 seg´ u n sus 5 ceros simples y complejos, de modo que todo sale por la v´ıa de exponenciales. >>F=1/(s^5+1) >>ilaplace(F) ans=1/5*Sum(exp(1/5*i*pi*(2*k-1))*exp(-exp(1/5*i*pi*(2*k-1))*t),k=1..5)

umero Nota 4.1 exp(1/5*i*pi*(2*k-1)),k=1..5 son las 5 ra´ıces complejas del n´

−1

Clase Pr´ actica 4

54

4.3

´ DE APLICACIONES A LA RESOLUCION ´ INICIAL PROBLEMAS DE CONDICION

Ejemplo 4.6 Resolver el siguiente problema de valor inicial usando transformadas de Laplace. x (t) + 3x(t) = 0, x(0) = 1.

Resoluci´ on del ejemplo 4.6. Se verá que el comando LAPLACE aplica automáticamente la propiedad

 −

x = sx

x(0)

>>syms t s lapx >>diffeqn=’diff(x(t),t)+3*x(t)’ % se introduce la ED diffeqn=diff(x(t),t)+3*x(t) >>a=laplace(diffeqn,t,s) % se transforma la ED a =s*laplace(x(t),t,s)-x(0)+3*laplace(x(t),t,s) >>a=subs(a,’laplace(x(t),t,s)’,lapx) % por comodidad a =s*lapx-x(0)+3*lapx en la expresion transformada se sustituye laplace(x(t),t,s) por lapx >>lapx=solve(a,lapx) % en a=0 se despeja lapx lapx =x(0)/(s+3) >>x=ilaplace(lapx,s,t) %se calcula x mediante x =x(0)*exp(-3*t)

La transformada inversa >>sol=subs(x,’x(0)’,1) en x(t) se sustituye x(0) por 1 sol =exp(-3*t)

Ejemplo 4.7 Resolver el siguiente problema de valor inicial x (t) + 3x (t) + 2x(t) = sen (t), x(0) = x  (0) = 0.

Resoluci´ on del ejemplo 4.7. >>syms s t lapx x >>diffeq=’diff(x(t),t$2)+3*diff(x(t),t)+2*x(t)-sin(t)’; >>a=laplace(diffeq,t,s) a =s*(s*laplace(x(t),t,s)-x(0))-D(x)(0)+... 3*s*laplace(x(t),t,s)-3*x(0)+ 2*laplace(x(t),t,s)-1/(s^2+1)

Clase Pr´ actica 4

55

>>a=subs(a,’laplace(x(t),t,s)’,lapx) a=s*(s*lapx-x(0))-D(x)(0)+3*s*lapx3*x(0)+2*lapx-1/(s^2+1) >>lapx=solve(a,lapx) lapx=(x(0)*s^3+s*x(0)+D(x)(0)*s^2+D(x)(0)+... 3*x(0)*s^2+3*x(0)+1)/(s^4+3*s^2+3*s^3+3*s+2) >>x=ilaplace(lapx,s,t) >>x =-exp(-2*t)*x(0)-1/5*exp(-2*t)exp(-2*t)*D(x)(0)+2*exp(-t)*x(0)+... 1/2*exp(-t)+exp(-t)*D(x)(0)3/10*cos(t)+1/10*sin(t) >>sol=subs(x,’x(0)’,0); >>sol=subs(sol,’D(x)(0)’,0) sol=-1/5*exp(-2*t)+1/2*exp(-t)-3/10*cos(t)+1/10*sin(t)

Ejercicio 4.4 Resolver el siguiente problema de valor inicial x (t) 3x (t) + 2x(t) = 2e−2t , x(0) = 1/6, x (0) = 5/3.

−

Ejemplo 4.8 Resolver el siguiente sistema x = x y y  = y + x x(0) = 1, y(0) = 2

− − − −

Resoluci´ on del ejemplo 4.8 >> syms s t lapx lapy % Se introducen las dos ED. >>diffeq1=’diff(x(t),t)+x(t)+y(t)’; >>diffeq2=’diff(y(t),t)+y(t)-x(t)’;

A continuacion se transforman las dos ED. >> a=laplace(diffeq1,t,s) a =s*laplace(x(t),t,s)-x(0)+laplace(x(t),t,s)+ laplace(y(t),t,s) >> b=laplace(diffeq2,t,s) b=s*laplace(y(t),t,s)-y(0)+laplace(y(t),t,s)- laplace(x(t),t,s) >>a=subs(a,’laplace(x(t),t,s)’,lapx); >>a=subs(a,’laplace(y(t),t,s)’,lapy) a=s*lapx-x(0)+lapx+lapy >>b=subs(b,’laplace(x(t),t,s)’,lapx);

Clase Pr´ actica 4

56 >>b=subs(b,’laplace(y(t),t,s)’,lapy) b=s*lapy-y(0)+lapy-lapx % En las expresiones a=0 y b=0 se despejan % las transformadas de x e y: lapx, lapy >>[lapx,lapy]=solve(a,b,lapx,lapy) lapx = (s*x(0)-y(0)+x(0))/(s^2+2*s+2) lapy =(s*y(0)+x(0)+y(0))/(s^2+2*s+2) % Se obtienen x e y mediante la transformada inversa >>x=ilaplace(lapx,s,t) x = exp(-t)*x(0)*cos(t)-exp(-t)*y(0)*sin(t) >>y=ilaplace(lapy,s,t) y =exp(-t)*x(0)*sin(t)+exp(-t)*y(0)*cos(t) % En las expresiones de x e y se sustituyen las condiciones iniciales >>x=subs(x,’x(0)’,1); >>x=subs(x,’y(0)’,-2) x =exp(-t)*cos(t)+2*exp(-t)*sin(t) >>y=subs(y,’x(0)’,1); y =exp(-t)*sin(t)-2*exp(-t)*cos(t)

Ejercicio 4.5 Resolver el siguiente sistema: x (t) = x  (t) 2y  (t) + 2y(t), y (t) = 3y(t) + x(t), y(0) = 0, x(0) = 2, x (0) = 1.

−

−

4.4 4.4.1

PROBLEMAS DE VALORES INICIALES CON DATOS DISCONTINUOS Funci´ on escal´ on unitario (Heaviside)

La función de salto unitario u(t) =

 

0

t < 0

1

t

≥0

se representa en notación MATLAB como Heaviside(t) y es de gran utilidad para representar otras funciones definidas a tramos. Veamos cómo se calcula su transformada: >>syms t s >>f=’Heaviside(t)’

Clase Pr´ actica 4

57

>>laplace(f,t,s) ans=1/s

on u(t a). Teniendo en cuenta la notaci´ on Ejemplo 4.9 Calcular la transformada de la funci´ τ a (f )(t) = f (t a)u(t a), y la terminolog´ıa al uso, podemos decir que u(t a) = τ a (1), es decir, on f = 1 al punto a. u(t a) es la trasladada paralelamente al eje OX de la funci´

−

−

−

−

−

Resoluci´ on del ejemplo 4.9. >>syms t s >>f=’Heaviside(t-1)’; >>laplace(f,t,s) ans=exp(-s)/s

Lo anterior es un resultado ya visto en clases, y fácil de obtener. Es oportuno recordar una importante propiedad vinculada a la función escalón unitario.

Propiedad de la traslaci´ on en t. ˜ y a > 0. Entonces Sea F (s) = f (s) f (t  a)u(t

−

− a)(s) = e−

as

F (s),

que se transforma en la siguiente igualdad cuando aplicamos en (4.3) la transformada inversa: as F (s) = f (t e−

− a)u(t − a).

on inicial Ejemplo 4.10 Resolver el siguiente problema de condici´ x (t) + x(t) = 1

− u(t − 2),

x(0) = 0. utilizando las herramientas MATLAB que simulan la transformada de Laplace. ¿Es posible resolverlo siguiendo los métodos anteriormente vistos?

Resoluci´ on del ejemplo 4.10. >>syms s t lapx >>diffeqn=’diff(x(t),t)+x(t)-1+Heaviside(t-2)’ diffeqn =diff(x(t),t)+x(t)-1+Heaviside(t-2) >>a=laplace(diffeqn,t,s) a=s*laplace(x(t),t,s)-x(0)+laplace(x(t),t,s) -1/s+exp(-2*s)/s

(4.3)

Clase Pr´ actica 4

58 >>a=subs(a,’laplace(x(t),t,s)’,lapx) a=s*lapx-x(0)+lapx-1/s+exp(-2*s)/s >>lapx=solve(a,lapx) lapx=(x(0)*s+1-exp(-2*s))/s/(1+s) >>x=ilaplace(lapx,s,t) x=x(0)*exp(-t)+1-exp(-t)-Heaviside(t-2)+ Heaviside(t-2)*exp(-t+2) >>sol=subs(x,’x(0)’,0) sol=1-exp(-t)-Heaviside(t-2)+Heaviside(t-2)*exp(-t+2)

4.4.2

Representaci´ on de funciones con discontinuidades de salto finito

a regida por el problema de valor inicial: Ejemplo 4.11 La corriente I (t) en un circuito RLC est´ I  (t) + 4I (t) = g(t), I (0) = I  (0) = 0, siendo

 − 1

g(t) =

1

0

Determinar I (t).

0 < t 1 < t t > 2

≤1 ≤2

on g(t) puede escribirse como g(t) = u(t) Nota 4.2 Obsérvese que la funci´

− 2u(t − 1) + u(t − 2).

Resoluci´ on del ejemplo 4.11. >>syms s t a I lapI >>diffeq=’diff(I(t),t$2)+4*I(t)-Heaviside(t)+2*Heaviside(t-1)-Heaviside(t-2)’; >>a=laplace(diffeq,t,s) a=s*(s*laplace(I(t),t,s)-I(0))-D(I)(0)+... 4*laplace(I(t),t,s)-1/s+2*exp(-s)/s-exp(-2*s)/s >>a=subs(a,’laplace(I(t),t,s)’,lapI) a=s*(s*lapI-I(0))-D(I)(0)+4*lapI-1/s+... 2*exp(-s)/s-exp(-2*s)/s >>lapI=solve(a,lapI) lapI=(s^2*I(0)+D(I)(0)*s+1-2*exp(-s)+exp(-2*s))/s/(s^2+4) >>I=ilaplace(lapI,s,t) I=I(0)*cos(2*t)+1/2*D(I)(0)*sin(2*t)+1/4-1/4*cos(2*t)-1/2*Heaviside(t-1)+... 1/2*Heaviside(t-1)*cos(2*t-2)+1/4*Heaviside(t-2)-1/4*Heaviside(t-2)*cos(2*t-4) >>sol=subs(I,’I(0)’,0) % condicion inicial I(0)=0 >>sol=subs(sol,’D(I)(0)’,0) % condicion inicial I’(0)=0

Clase Pr´ actica 4

59

sol=1/4-1/4*cos(2*t)-1/2*Heaviside(t-1)+1/2*Heaviside(t-1)*cos(2*t-2)+... 1/4*Heaviside(t-2)-1/4*Heaviside(t-2)*cos(2*t-4) % fin de la resolucion del ejemplo 4.11

Fig. 4.2: Representación gráfica de la entrada g(t) y la salida I (t) La representació n gráfica conjunta de la entrada g(t) y la salida I (t), que pueden verse en la figura 4.2, se obtiene como sigue:(1 >>maple(’g:=t->1-2*Heaviside(t-1)+Heaviside(t-2)’); >>maple(’I:=t->1/4-1/4*cos(2*t)-1/2*Heaviside(t-1)+... 1/2*Heaviside(t-1)*cos(2*t-2)+1/4*Heaviside(t-2)-... 1/4*Heaviside(t-2)*cos(2*t-4)’); >>x=linspace(0,6,500); >>for k=1:500,a=num2str(x(k));... h(k)=str2num(maple(strcat(’evalf(g(’,a,’),12)’))); j(k)=str2num(maple(strcat(’evalf(I(’,a,’),12)’)));end >>plot(x,h,’r.’);hold on;plot(x,j,’r.’) >>axis([0 6 -2 2]);title(’funciones g(t) e I(t)’)

La delta de Dirac Sea a

∈ R. La delta de Dirac δ (t − a), se caracteriza por las dos propiedades siguientes:  a 0 t= δ (t − a) = ∞ t = a



1

Lo más sencillo es copiar lo anterior en el Editor-Debugger, b orrar los s´ımbolos >> y ejecutar.

(4.4)

Clase Pr´ actica 4

60 y



∞

f (t)δ (t

− a)dt = f (a),

(4.5)

−∞ para cualquier función f (t) que sea continua en un intervalo abierto que contiene a t = a. Observemos que δ (t a) no es una función en el sentido usual ya que toma el valor punto. De la propiedad (4.5) se obtiene que

−

∞ en un

δ ( t a)(s) = e −as

(4.6)

−

La Delta de Dirac se relaciona con la función escalón unitario de la siguiente forma



t

−∞

δ (u

− a)du =



0 1

t
≥

− a),

(4.7)

lo que, formalmente, significa que δ (t a) = u  (t a). La delta de Dirac sirve para modelizar fenómenos en los que intervienen entidades f´ısicas de ´ındole eléctrica o mecánica, cuyas magnitudes alcanzan en poco tiempo valores relativamente grandes. Ello ocurre cuando se aplica una carga muy pesada concentrada en un punto de una viga, un golpe de martillo a un sistema masa-resorte, picos de voltaje, etc. Debido a razones históricas y pr´ acticas, la delta de Dirac se presenta tal como acabamos de hacer en este guión, pero el alumno puede intentar conciliar las relaciones (4.5-4.7) asumiendo que dδ a (t) = δ (t a)dt es la medida de probabilidad concentrada en t = a. Después de un per´ıodo oscuro pero exitoso, el uso de δ (t a) fue rigurosamente justificado en 1944 por Laurent Schwarz, al crear las distribuciones lo que le valió recibir en 1950 la Medalla Fields. Sobre este tema se recomienda leer la sección 7.8 del libro de Nagle-Saff [8], la 7.8 del Nagle-SaffSnider [9], o la 7.6 de D. G. Zill [11].

−

−

−

−

Ejemplo 4.12 Una masa sujeta a un resorte se suelta a partir del reposo 1 m por debajo de la posici´ on de equilibrio del sistema resorte-masa, y empieza a vibrar. Después de π segundos la masa es golpeada por un martillo que ejerce un impulso sobre la masa. El sistema est´ a regido por el problema de valor inicial siguiente x (t) + 9x(t) = 3δ (t x(0) = 1, x (0) = 0,

− π)

donde x(t) representa el desplazamiento a partir del equilibrio en el instante t. Hallar x(t).

Resoluci´ on del ejemplo 4.12 >>syms s t lapx x % lapx=X >>diffeq=’diff(x(t),t$2)+9*x(t)-3*Dirac(t-pi)’ diffeq =diff(x(t),t$2)+9*x(t)-3*Dirac(t-pi)

Clase Pr´ actica 4

61

>>a=laplace(diffeq,t,s) a=s*(s*laplace(x(t),t,s)-x(0))-D(x)(0)+9*laplace(x(t),t,s)-3*exp(-pi*s) >>a=subs(a,’laplace(x(t),t,s)’,lapx) a=s*(s*lapx-x(0))-D(x)(0)+9*lapx-3*exp(-pi*s) >>lapx=solve(a,lapx) lapx=(s*x(0)+D(x)(0)+3*exp(-pi*s))/(s^2+9) >>lapx=subs(lapx,’x(0)’,1) % cond. inicial x(0)=1 lapx=(s+D(x)(0)+3*exp(-pi*s))/(s^2+9) >>lapx=subs(lapx,’D(x)(0)’,0) %cond. inicial x’(0)=0 lapx=(s+3*exp(-pi*s))/(s^2+9) >>x=ilaplace(lapx,s,t) x=cos(3*t)-Heaviside(t-pi)*sin(3*t) % fin de la resolucion del ejemplo 4.12

Ejemplo 4.13 Resolver la siguiente ED de coeficientes variables, sujeta a ciertas condiciones. tx (t) + (t 1)x (t) x(0) = 5, lim x(t) = 0. t→+∞

−

− x(t) = 0,

on alguna sobre x  (0)! En su lugar Nota 4.3 Esta ED es de 2do orden...¡pero no se establece condici´ se impone una condici´ on en el infinito... Notar que t = 0 es una singularidad (¿por qué?). Quiz´ as esto ultimo ´ esté asociado a esas “extra˜ nas condiciones” que se imponen a la soluci´ on. Una idea: si no, la ED es casi de orden uno...y para t grande es indudablemente t [ ε, ε] con ε > 0 muy peque˜ una ecuaci´ on de orden dos.

∈−

Resoluci´ on del ejemplo 4.13 >>syms s t a x >>diffeq=’t*diff(x(t),t$2)+(t-1)*diff(x(t),t)-x(t)’; >>a=laplace(diffeq,t,s) a=-2*s*laplace(x(t),t,s)+2*x(0)-... s*(laplace(x(t),t,s)+s*diff(laplace(x(t),t,s),s))-... 2*laplace(x(t),t,s)-s*diff(laplace(x(t),t,s),s) >>X=dsolve(’-2*s*X+2*5-s*(X+s*DX)-2*X-s*DX=0’, ’s’) X=5/(s+1)+1/(s+1)/s^2*C1%Ver nota sobre X y x >>x=ilaplace(X,s,t) x=5*exp(-t)+C1*(t-1+exp(-t))

on habitual escribiendo X en lugar de x. Nota 4.4 Hemos sido consecuentes con nuestra notaci´ Aunque esto no es esencial a los efectos de la compilaci´ on y el c´ alculo, s´ı facilita la lectura del c´ odigo. Es un ahorro innecesario de memoria usar la misma variable x para ambos papeles. Notar que X se crea por asignaci´ on y no por la declaraci´ on SYMS.

Clase Pr´ actica 4

62

Finalmente, teniendo en cuenta el comportamiento de x en el punto del infinito se sigue que C 1 debe ser cero. Por lo tanto la solución del ejemplo 4.13 es x(t) = 5e−t .

4.5

TRANSFORMADA DE LAPLACE DE UNA ´ PERIODICA ´ FUNCION

En primer lugar recordamos la expresi´ on obtenida en clase para calcular la transformada de Laplace de una función periódica f (t) con per´ıodo T .



T

˜ = f (s)

e−st f (t)dt

0

1

Ts

− e−

.

A continuación la escribimos de nuevo utilizando ahora la notación MATLAB. 1/(1-exp(-T*s))*(int(exp(-s*t)*f(t),0,T)

on f (t), peri´ odica de per´ıodo 2π, definida por Ejemplo 4.14 Obtener la transformada de la funci´

 −

20

f (t) =

20

0

≤t≤π

π < t < 2π

Resoluci´ on del ejemplo 4.14 >>Lf=1/(1-exp(-2*pi*s))*(int(20*exp(-s*t),0,pi) -int(20*exp(-s*t),pi,2*pi)); >>Lf=simplify(Lf) Lf=-20*(exp(-pi*s)-1)/s/(exp(-pi*s)+1)

4.6

´ DE TRANSFERENCIA LA FUNCION

on de transferencia H (s) de un sistema lineal se define como la razón entre la transformada La funci´ de Laplace de la salida y(t) y la transformada de Laplace de la entrada g(t), asumiendo que las condiciones iniciales son nulas. Supongamos que el sistema está gobernado por la EDL de coeficientes constantes siguiente ay  + by + cy = g(t), t > 0, (4.8) y(0) = y  (0) = 0. Al aplicar transformadas en ambos miembros de la ecuación, y reflejar las condiciones iniciales nulas, se obtiene (as2 + bs + c)Y (s) = G(s).

Clase Pr´ actica 4

63

Luego, obtenemos la siguiente fórmula para H (s) en términos del polinomio caracter´ıstico de la ED H (s) =

1 Y (s) = 2 G(s) as + bs + c

˜ de donde, si h = H , se sigue que Y (s) =

G(s) ˜ = G(s)H (s) = g˜h. 2 as + bs + c

(4.9)

Notar que H (s) no depende de la entrada g(t). on de respuesta al impulso como Definamos ahora la llamada funci´ h(t) = L −1 H (t) Notar que h(t), al igual que H (s), tampoco depende de la entrada g(t). El alumno está en capacidad de deducir de (4.9) que la solución y(t) de la ecuación (4.8) está dada por (aplicar teorema de convolución) (4.10) y(t) = g h (t)

∗

Lo anterior significa que h(t) juega el papel de núcleo resolvente de la ecuación (4.8). Se halla de una vez, y luego, convolucionado con cualquier g(t) de entrada nos produce la respuesta o salida y(t). alculo directo halle la funci´ on de transferencia H (s) del sistema. Con la Ejercicio 4.6 Mediante c´ ayuda del MATLAB encuentre la funci´ on de respuesta al impulso h(t) de los siguientes PVI, con condiciones nulas. (i) y  + 9y = g(t) (iii) y  y 6y = g(t) (v) y  + 2y + 5y = g(t)

− −

(ii) y  9y = g(t) (iv) y  + 2y  15y = g(t) (vi) y  4y + 5y = g(t)

− − −

Los resultados que se piden ya sabemos que no dependen de g(t). Usar la f´ ormula (4.10) para calcular 2 la soluci´ on y(t) de (i) cuando g(t) = 4t y de (vi) cuando g(t) = exp(3t). (Para calcular el producto de convoluci´ on ver las instrucciones que se dan a continuaci´ on.)

Indicaciones sobre el ejercicio 4.6 El siguiente código puede ejecutarse en la l´ınea de comandos, o puede confeccionarse un programa. >>syms s t u >>h=ilaplace(H,s,t); >>f=subs(h,t,t-u); >>g=subs(g,t,u); >>P=f*g; >>y=int(P,u,0,t)

Clase Pr´ actica 4

64

Nota 4.5 Tener en cuenta que eit

it

eit + e−it sen (t) = , cos(t) = 2i 2 et e−t et + e−t senh (t) = , cosh(t) = 2 2 Usar tentativamente el comando SIMPLIFY si fuese necesario

− e− −

´ CLASE PRACTICA 5 M´ etodos num´ ericos para PVI (I)

OBJETIVOS Presentar los métodos numéricos de Euler y Euler mejorado para aproximar la solución de problemas de valores iniciales. Resolver problemas concretos, obteniendo soluciones numéricas y estimados del error, mediante el manejo de programas domésticos MATLAB basados en los algoritmos correspondientes.

5.1

FUNCIONES ESPECIALES

En Matemática Aplicada, F´ısica e Ingenier´ıa aparecen con frecuencia ED de segundo orden. Entre ellas destacamos la Ecuaci´ on de Bessel. Esta ED ha sido ampliamente estudiada y sus soluciones, asi como las de otras EEDD que surgen en las aplicaciones, reciben el nombre de FUNCIONES ESPECIALES. La ED lineal de segundo orden x2 y + xy + (x2

2

− q )y = 0

(5.1)

donde q es un parámetro fijo, se llama ecuaci´ ´ nica singularidad on de Bessel de orden q y tiene una u regular en x = 0. Aplicando el método de Frobenius se obtienen, si 2q / N, las 2 soluciones linealmente independientes siguientes: ∞ ( 1)n x ( )2n+q J q (x) = n!Γ(1 + q + n) 2 n=0 ∞ ( 1)n x ( )2n−q J −q (x) = n!Γ(1 q + n) 2 n=0

∈

 

−

− −

que son las funciones de Bessel de primera especie de órdenes q y q respectivamente. Las Funciones de Bessel verifican diversas relaciones de recurrencia, a saber:

−

J q+1 (x) + J q−1 (x) = 2q/xJ q (x), J q+1 (x)

− J − (x) = 2J  (x). q 1

65

q

Clase Pr´ actica 5

66

Si q Z, entonces J q (x) y J −q (x) son linealmente dependientes. Adem´ as J −q (x) = ( 1)q J q (x). Si q = 1/2, 3/2, 5/2, , las funciones de Bessel J q se pueden expresar como una combinación finita de senos, cosenos y potencias de x; es decir que son funciones elementales. En concreto se cumple que: 1/2 2 sen(x) J 1/2 (x) = πx y 1/2 2 cos(x). J −1/2 (x) = πx Las Funciones de Bessel J q para q N son anal´ıticas en 0. Además son funciones pares si q es par e impares si q es impar. A continuación se representan las funciones de Bessel J 0 , J 1 , J 1/2 y J 3/2 .

∈

±

±

±

−

···

   

∈

5.1.1

Funciones de Bessel

Dentro del grupo de funciones matemáticas especializadas que contiene MATLAB se encuentran las funciones de Bessel de primera clase: besselj(n, x). A continuación se representan la funciónes de Bessel J 1 y J 0 . >>ezplot(’0+x*0’,[-20,20]) >>hold on; ezplot(’besselj(1,x)’,[-20,20]) >>ezplot(’0+x*0’,[-10,10]) >>hold on; ezplot(’besselj(0,x)’,[-10,10])

Fig. 5.1: Gráfica de J 1

Fig. 5.2: Gr´ afica de J 0

Funciones de Bessel J 0 y J 1 Se representan conjuntamente en [ 20, 20] las dos funciones de Bessel J 0 y J 1 . Obsérvese cómo se alternan los ceros de J 0 y J 1 .

−

Clase Pr´ actica 5

67

>>x=linspace(-20, 20,300); >>y=besselj(0,x); >>z=besselj(1,x); >>plot(x, y,’b.’,x,z, ’r-’,x,zeros(1,max(size(x))),’k’)

Funciones de Bessel J 1/2 y J 3/2 Ejecutar en la l´ınea de comandos el código >>subplot(2,1,1),ezplot(’besselj(3/2,x)’,[0,100]) >>subplot(2,1,2),ezplot(’besselj(1/2,x)’,[0,100])

Para lograr los efectos que se muestran en la figura 5.4 deben usarse los recursos que ofrece la propia ventana gráfica.

Fig. 5.3: Gráficas superpuestas de J 1 y J 0

5.2

Fig. 5.4: Gr´ aficas de Bessel J 3/2 y J 1/2

ESTUDIO DE LAS SOLUCIONES DE UNA ED

Los métodos a utilizar se pueden clasificar como: (1) Métodos anal´ıticos de cuadraturas exactas (Temas 2, 3, 4 y 5). (2) Métodos cualitativos. (3) Métodos numéricos (Clases Prácticas 5, 6 y 7). Los métodos exactos hasta ahora estudiados tienen grandes limitaciones, pues la mayor´ıa de las ED que se encuentran en las aplicaciones no se pueden resolver de forma exacta, mediante fórmulas

Clase Pr´ actica 5

68

impl´ıcitas o expl´ıcitas.(1 Para compensar estas limitaciones se introducen los métodos numéricos que son el objeto de esta práctica.

5.2.1

M´ etodos num´ ericos de resoluci´ on de PVI

Consideremos el siguiente PVI y = f (x, y), y(x0 ) = y 0 , siendo x 0

(5.2)

∈ [a, b] y f ∈ C ([a, b] × R).

∂f son continuas en un rectángulo que contiene al punto (x0 , y0 ), de modo que ∂y el problema (5.2) tiene solución u ńica. Se asume que f y

Estrategia de los m´ etodos num´ ericos: Si y(x) es la solución exacta del problema (5.2), se trata de obtener valores aproximados y k de y en ciertos puntos x k de [a, b]. Se trata de crear una tabla de valores de la variable dependiente y correspondientes a valores conocidos de la variable independiente x en un intervalo [x0 , X ] [a, b]. El primer paso para obtener una solución num´ erica consiste en determinar una partici´ o n del subintervalo [x0 , X ], esto es, x0 < x1 < < xN = X.

⊂

···

Se supone que los puntos están distribuidos uniformemente, de modo que xk+1

−x

k

= h, k = 0, 1,...,N

− 1.

Los métodos num´ ericos que estudiaremos permiten determinar los valores yk que constituyen una aproximaci´ on del valor de la solución de (5.2) en los puntos x k .

Definici´ on 5.1 El conjunto

{y }

k k=0,...N

se llama Soluci´ on Numérica de (5.2) relativa al paso h.

5.2.2

M´ etodo de Euler o de la tangente

Reconsideremos el problema (5.2). Sabemos que el punto P = (x0 , y0 ) está en la curva solución del mismo y que la pendiente de la tangente a dicha curva en P es y  (x0 ) = f (x0 , y0 ). El método de Euler aproxima el valor de la soluci´ on exacta en x cercano a x 0 usando la ecuación de la tangente en el punto P = (x0 , y0 ): y(x) = y 0 + f (x0 , y0 )(x 1

−x ) 0

(5.3)

No se pueden resolver de forma exacta porque el catálogo disponible de “funciones elementales” no alcanza para obtener una fórmula cerrada, o porque dicha expresión es “muy dif´ıcil” de encontrar.

Clase Pr´ actica 5

69

Si se denotan por y1 , , yN los valores aproximados correspondientes a los valores exactos y(x1 ), entonces, aplicando (5.3) se sigue que dado y 0 se calcula y1 mediante

···

N ),

··· , y(x

y1 = y 0 + f (x0 , y0 )h, e y 2 mediante y2 = y 1 + f (x1 , y1 )h, y asi sucesivamente hasta y N yN = y N −1 + f (xN −1 , yN −1 )h, En forma compacta se obtiene el siguiente esquema yk+1 = y k + hf (xk , yk ), k = 0,

5.2.3

··· , N − 1.

(5.4)

Otra interpretaci´ on del m´ etodo de Euler

Integrando ambas partes de la ecuación y  (x) = f (x, y(x)) en el intervalo [xk , xk+1 ], obtenemos



xk+1

y(xk+1 )

− y(x ) = k



xk+1

y  (x)dx =

xk

f (x, y(x))dx,

xk

y sustituyendo la integral de la derecha por el área de cierto rectángulo Q se obtiene que y(xk+1 )

≈ y(x ) + f (x , y(x ))(x − x ). k

k

k

k+1

k

Si ponemos el paso h = x k+1 xk y sustituimos y(xk ) por yk obtenemos de nuevo el esquema (5.4) del método de Euler. (5.5) yk+1 = y k + f (xk , yk )h, k = 0, , N 1.

−

··· −

La fórmula (5.5), llamada “ecuaci´ on en diferencias finitas”, permite predecir yk+1 en términos del valor de yk . El dato o condici´ on inicial (CI) es y0 y desde él marchamos en una dirección haciendo cada vez una nueva estimación sobre la base de predicciones anteriores. Esto último introduce desviaciones en nuestros sucesivos cálculos, provocando que éstos sean m´ as imprecisos a medida que nos alejamos del punto inicial x 0 . ´ La formulaci´ on (5.5) anterior corresponde al llamado método de Euler. Este resulta ser el más sencillo entre los métodos que se establecen mediante diferencias finitas para resolver problemas de valores iniciales. La simplicidad del m´ etodo de Euler lo hace especialmente adecuado para introducirnos en el tema num´ erico-experimental, si bien su utilización profesional es nula. Sea y(x; h) la solución aproximada obtenida con paso h, mediante un método numérico, por ejemplo Euler, de un PVI cuya solución exacta es y(x).

Clase Pr´ actica 5

70

Principio 5.1 Si agregamos m´ as puntos en la malla (xk ), para que h sea m´ as peque˜ no, entonces la expectativa es que los errores y(x) y(x; h) sean m´ as peque˜ nos.(2

|

−

|

¿Qué significa decir que el método utilizado converge? (ver el 7.6 de [3])

§

Definici´ on 5.2 Si lim y(x; h) = y(x),

h→0

entonces decimos que el método converge. Lo anterior significa que las soluciones y(x; h) de las ecuaciones en diferencias finitas que caracterizan al método considerado, aproximan a la soluci´ on exacta y(x) de la ED cuando h 0. Luego, teóricamente, deber´ıamos obtener mejores respuestas considerando valores cada vez más peque˜ nos de h. Desafortunadamente, la aritmética de precisi´ on finita del ordenador no lo permite. Para h muy pequeño los valores y k obtenidos podr´ıan estar muy lejos del valor exacto y(xk ).(3 En general, la velocidad con la que y(x; h) se acerca a la solución exacta y(x) se expresa de la siguiente forma: y(x) y(x; h) = O(h p ),

→

−

lo que significa que p

|y(x) − y(x; h)| ≤ M |h| para h pequeño y M > 0 dependiendo de x.

on anterior se le llama orden del método. Definici´ on 5.3 Al exponente p de la relaci´

Nota 5.1 El método de Euler es de orden p = 1.

5.2.4

Ejemplos en los que se aplica el m´ etodo de Euler

on aproximada del PVI Ejemplo 5.1 Obtener una soluci´ y = y 2 + 2x

4

− x , 0 ≤ x ≤ 1, con y(0) = 0

En el programa EULER1.m que indicamos a continuaci´ on, se calcula la soluci´ on aproximada del PVI y = y 2 + 2x

4

− x , a ≤ x ≤ b,

con y(a) = y 0 . Se calcula también el error cometido comparando con la soluci´ on exacta del PVI que es y = x 2 . 2

Un principio es un enunciado metalingu´ıstico, mediante el cual se establecen expectativas, se resumen resultados, se sugieren acciones,... Los principios surgen de la experiencia y conforman el criterio experto. 3 Es decir, la mantisa correspondiente a yk , podr´ıa no tener d´ıgitos significativos coincidentes con la representación de punto flotante de y (xk ).

Clase Pr´ actica 5

71

Fig. 5.5: Método de Euler

Fig. 5.6: Interpretación del método de Euler

Resoluci´ on del ejemplo 5.1 Programa EULER1.m Datos: n es el número de puntos de la malla, a, b: extremos del intervalo dominio. Resultados (salida): arreglo y = [y(1) y(n)] con los valores aproximados de la soluci´ on en los puntos x(k) = a + (k 1) h con k = 1, , n y donde h = (b a)/(n 1) es el paso. Error máximo maxk y(k) y(x(k))

···

−

− ∗ − {

−

··· }

>>clear x y >>n=20 >>a=0, b=1, y0=0 >>tic % Se comienza a calcular el tiempo de ejecucion >>x=linspace(a,b,n), y(1)=y0; h=(b-a)/(n-1) >>fork=1:n-1 >>y(k+1)=y(k)=+h*(y(k)^2+2*x(k)-x(k)^4); >>end >>error=max(abs(x.^2-y)); >>toc % Se termina el computo del tiempo de ejecucion >>y >>disp(’El error maximo es: ’), disp(error) >>xx=linspace(a,b,50); >>plot(x,y,’k*’,xx,xx.^2,’b-’)

Clase Pr´ actica 5

72 >>title(’Solucion exacta y aproximada de y’=y^2+2*x-x^4’) >>gtext(’Solucion exacta : --’) >>gtext(’Solucion aproximada : *’)

Soluci´ on: n = 20, a = 0, b = 1, y0 = 0, h = 5.263157894737e 02 El error máximo es: 8.093460004977e 02 Los valores aproximados se obtienen en la l´ınea de comandos seg´ un se indica a continuación.

−

−

>>[y(1:10)’ y(11:20)’] ans = 0 2.4775e-001 0 3.0234e-001 5.5398e-003 3.6218e-001 1.6615e-002 4.2720e-001 3.3218e-002 4.9729e-001 5.5333e-002 5.7235e-001 8.2942e-002 6.5225e-001 1.1602e-001 7.3682e-001 1.5454e-001 8.2584e-001 1.9847e-001 9.1907e-001

Ejemplo 5.2 Obtener en el ejemplo 5.1 las diferentes aproximaciones correspondientes a n = 10, 100, 200 y 400. Calcular en cada caso el error comparando con la soluci´ on exacta, y relacionar éste con el paso h. Apreciar gr´ aficamente los resultados obtenidos. Resultados del ejemplo 5.2 n 10 100 200 400

paso h Error max 0.111111 0.1579 0.010101 0.0166 0.005025 0.0083 0.002506 0.0042

on aproximada del PVI Ejercicio 5.1 Modificar el programa EULER1 para calcular una soluci´ y = 2xy2 , 0

≤ x ≤ 1,

con y(0) =

−1.

Obtener los errores cometidos, comparando con la soluci´ on exacta, considerendo los pasos h1 = 2−3 y h2 = 2−5 .

Clase Pr´ actica 5

73

Fig. 5.7: Soluciones exactas y aproximadas de los ejemplos 5.1 y 5.2

Fig. 5.8: Soluciones exactas y aproximadas del ejemplo 5.2 para n = 100, 200 A continuación se presenta el siguiente código modificado para el método de Euler que llamaremos EULER2.m. para resolver la ecuaci´ on y  = F (x, y) con la CI y(a) = y0 , en el intervalo [a, b], en n puntos x(k) = a + (k 1)h,

−

con k = 1, , n y donde el paso es h = (b a)/(n 1). La función F se define previamente mediante el comando inline. Se obtiene el gráfico y el vector y de la solución. El programa EULER2.m es un procedimiento de tipo FUNCTION y uno de sus argumentos es una función que se debe crear previamente. Una de las aplicaciones más u ´tiles de los ficheros .M se presenta en la definición de funciones mediante la directiva FUNCTION cuya sintaxis es FUNCION PARAMETROS SALIDA=NOMBRE FUNCION(parametros entrada) El código de EULER2.m es el siguiente

···

−

−

Clase Pr´ actica 5

74

function y=euler2(F,a,b,y0,n) h=(b-a)/(n-1) x=linspace(a,b,n);y(1)=y0; for k=1:n-1 y(k+1)=y(k)+h*F(x(k),y(k)); end disp(’Solucion de dy/dx=F(x,y), donde’), disp(’F es la funcion inline’), disp(F) disp(’El numero de puntos es n=’), disp(n) disp(’Paso h=’), disp(h) plot(x,y,’k’) title(’Solucion de la ED dy/dx=F(x,y)’)

on de la ED Ejemplo 5.3 Calcular aproximadamente, usando EULER2.m, la soluci´ y  = y + sen (x) con y(0) = 1, 0

≤ x ≤ 1.

Fig. 5.9: Solución aproximada del ejemplo 5.3

Soluci´ on del ejemplo 5.3 >>F=inline(’y+sin(x)’) F=Inline function: F(x,y) = y+sin(x) >>y=euler2(F,0,1,1,10)

Clase Pr´ actica 5

75

El numero de puntos es n = 10 Paso h = 0.11111111111111 >>y’= 1.00000000000000 1.11111111111111 1.24688819329123 1.40991996404094 1.60293270413395 1.82880926720753 2.09061200642747 2.39160998526049 2.73531085875018 3.12549783431139

Nota 5.2 Obsérvese que el programa EULER2 no permite juzgar la calidad de los resultados obtenidos. on de la ED Ejercicio 5.2 Aproximar, usando EULER2.m, la soluci´ 1 (y + y2 ), x y(1) = 1, 1 x y =

≤ ≤2

(5.6)

tomando el paso h = 2−2 .

5.2.5

Precisi´ on de los resultados

Todo cálculo aproximado debe tener en cuenta un indicador que permita evaluar la calidad de los resultados. En nuestro caso se trata de calcular la solución aproximada de una ED y de poder decidir sobre la precisión del resultado. El error exacto es imposible de conocer, de modo que se trata de obtener una estimación del mismo. Si dicha estimación se obtiene después de realizar los cálculos y se basa en la utilización de éstos, se dice que es una estimación del error a posteriori. El siguiente programa EULER3 es una versió n más flexible que la anterior, pues permite que el usuario defina el dato F (x, y), el intervalo y la CI. Como salida se obtienen los valores aproximados y las gráficas de dos resultados calculados con diferentes mallas de puntos, lo que permite estimar el error a posteriori si nos basamos en el siguiente on exacta de un PVI, entonces d principio: Si d es la distancia entre dos aproximaciones de la soluci´ es una estimaci´ on de la distancia entre una cualquiera y la soluci´ on exacta.(4 El programa EULER3.m es un procedimiento de tipo FUNCTION y uno de sus argumentos es una función que se debe crear previamente. 4

Se trata de un enunciado obviamente impreciso y aparentemente peligroso. Es frecuente en la práctica cuando no hay otra opción. Se emplea en los comandos Matlab: ODE45 y ODE32.

Clase Pr´ actica 5

76

Resuelve la ecuación y  = F (x, y) con la condici´ on inicial y(a) = y 0 , en el intervalo [a, b], en N puntos x(k), (k = 1, , N ), donde

···

h = (b

− a)/(N − 1) y N = n, 2n − 1.

La función F se define previamente mediante el comando INLINE.

Procedimiento para estimar el error Se estima a posteriori el error tomando la diferencia máxima entre las soluciones obtenidas con los pasos h = (b a)/(n 1) y h1 = h/2 = (b a)/(2n 2).

−

−

−

−

Salida. Se ofrece la estimación del error a posteriori y los gráficos de ambas aproximaciones numéricas. Código de EULER3.M function [error,x1,y1]=euler3(F,a,b,y0,n) N=n; h=(b-a)/(N-1); x=linspace(a,b,N); y(1)=y0; donde y es la solución con N = n puntos. for k=1:N-1 y(k+1)=y(k)+h*F(x(k),y(k)); end N=2*n-1; h=(b-a)/(N-1) x1=linspace(a,b,N) y1(1)=y0 donde y1 es la solución con N = 2n 1 puntos. for k=1:N-1 y1(k+1)=y1(k)+h*F(x1(k),y1(k)); end for k=1:n e(k)=abs(y(k)-y1(2*(k-1)+1); end error=max(e) plot(x,y,’b’,x1,y1,’r’) legend(’Aprox. con n puntos’, ’Aprox. con 2n-1 puntos’) title(’Solucion de la ecuacion dy/dx =F(x,y)’)

−

Clase Pr´ actica 5

77

Ejemplo 5.4 Obtener, usando EULER3.m, una soluci´ on aproximada de la ED y = y 2 + 2x

4

− x , 0 ≤ x ≤ 1,

con la C.I. y(0) = 0. Estimar el error.

Resoluci´ on del ejemplo 5.4 >>F=inline(’y.^2+2.*x-x.^4’) Inline function: F(x,y) = y.^2+2.*x-x.^4 >>[error,x1,y1]=euler3(F,0,1,0,10) h =1.111111111111111e-001 h =5.555555555555555e-002 error = 7.285473724418645e-002

Fig. 5.10: Soluciones aproximadas del ejemplo 5.4

Ejercicio 5.3 Sea el PVI siguiente y =

2

−y , 0 ≤ x ≤ 1/2,

con y(0) =

−1.

• ¿A qué tipos conocidos de ecuaciones corresponde la anterior ED? • ¿Por qué no es lineal? • ¿Qué métodos utilizar´ıas para intentar hal lar una soluci´ on exacta? • Usando EULER3 obtener las soluciones aproximadas con N = 2 estimar el error cometido.

3

y N = 24

(5.7)

− 1 puntos y

• Calcular el error exacto y compararlo con el error estimado teniendo en cuenta que la soluci´ on exacta es y = 1/(x − 1).

Clase Pr´ actica 5

78

5.3

´ METODO DE EULER MEJORADO

Para dise˜ nar otro método numérico más eficiente que el de Euler, utilizaremos una nueva fórmula de integraci´ on num´ erica de mayor orden de exactitud que la utilizada en el método de Euler.

5.3.1

La f´ ormula de los trapecios

La fórmula de los trapecios consiste en sustituir el valor exacto de la integral



b

f (x)dx,

a

por el número (b

− a) (f (a) +2 f (b)) ,

es decir, el área del trapecio cuyos vértices son (a, 0), (a, f (a)), (b, f (b)) y (b, 0). Con ello se trata de efectuar cálculos más precisos en la aproximación de la integral sobre cada uno de los subintervalos [xk , xk+1 ]



xk+1

f (x, y(x))dx.

xk

Notar otra vez que el integrando de la expresi´ on anterior tiene dentro de s´ı a la incógnita y(x) de nuestro problema.

Fig. 5.11: Regla del trapecio En resumen y dicho en términos geométricos, el método de Euler consistió en sustituir el área bajo la curva de g(x) = f (x, y(x)) por el área de cierto rectángulo, mientras que ahora, en lugar del

Clase Pr´ actica 5

79

rectángulo utilizaremos un trapecio para obtener la siguiente expresión(5





xk+1

f (x, y(x))dx

xk



h f (xk , y(xk )) + f (xk+1 , y(xk+1 )) , 2

≈

donde h = x k+1 xk . Surge as´ı un método de cálculo aproximado cuya formulaci´ on es la siguiente.

−





h yk+1 = y k + f (xk , y(xk )) + f (xk+1 , yk+1 ) . 2 Se trata de un método impl´ıcito ya que la incógnita yk+1 , no aparece despejada en la fórmula. Además se prueba que el error global de discretización E (h) cumple E (h) =

2

O(h ),

y por tanto es de orden p = 2.

5.3.2

Deducci´ on del m´ etodo de Euler mejorado

El método de Euler mejorado se obtiene realmente como el método predictor-corrector siguiente (1) Se predice y k+1 mediante Euler

∗ = y + hf (x , y ). yk+1 k k k (2) Se corrige mediante la anterior relación h yk+1 = y k + 2



      

f (xk , yk ) + f (xk + h, yk + hf (xk , yk )) x y∗ k+1

k+1

Resultando definitivamente el algoritmo de Euler mejorado h yk+1 = y k + [f (xk , yk ) + f (xk+1 , yk + hf (xk , yk ))] . 2

5.3.3

Euler versus Euler mejorado

El siguiente programa, EULER M, nos permite comparar entre s´ı los dos métodos estudiados anteriormente. El procedimiento consiste en tratar de resolver una ED, cuya solución es conocida, mediante Euler y Euler mejorado. El experimento culmina con la estimaci´ on del error global con respecto a ambos métodos, y el trazado en pantalla de los respectivos gráficos. 5

Ver sección 3.5 del Nagle Saff y la figura 5.11.

Clase Pr´ actica 5

80

Descripci´ on del programa EULER M El programa EULER M resuelve la ecuación y  = F (x, y) con la C.I. y(a) = y 0 , en el intervalo [a, b], en n puntos x(k), (k = 1, , n), con el paso h = (b a)/(n 1), utilizando los métodos Euler y Euler mejorado. Las funciones F (campo de direcciones) y f (solución exacta conocida) se definen previamente mediante el comando INLINE. El programa estima los respectivos errores absolutos y los ofrece como salida, según el formato que más aba jo se explica. Estos errores se obtienen directamente en términos del máximo de la distancias entre cada estimación yk y el valor exacto f (xk ). Los correspondientes gráficos se despliegan en una misma ventana para facilitar la apreciación visual.

· ··

−

−



C´ odigo EULER M function [es,em]=euler_m(F,f,a,b,y0,n) h=(b-a)/(n-1); x=linspace(a,b,n); phi=f(x); ys(1)=y0; for k=1:n-1 ys(k+1)=ys(k)+h*F(x(k),ys(k)); end ym(1)=y0; for k=1:n-1 ym(k+1)=ym(k)+(h/2)*(F(x(k),ym(k))+F(x(k)+h,ym(k)+h*F(x(k),ym(k)))); end % Errores absolutos em=max(abs(ym-phi)); es=max(abs(ys-phi)); plot(x,phi,’g*’,x,ys,’k’,x,ym,’b’) legend(’Solucion exacta’,’Euler’,’Euler mejorado’)

on aproximada de la ED Ejemplo 5.5 Obtener una soluci´ y + y/x y2 log(x) = 0,

−

en [1, 3] con la C.I. y(1) = 1, usando los métodos de Euler y de Euler mejorado. Para el cálculo del error, usar que la solución exacta es f (x) =

2

−2/x/(log(x) − 2),

la cual se puede obtener usando el comando DSOLVE. >>dsolve(’Dy+y/x-y^2*log(x)=0’,’y(1)=1’,’x’) ans =-2/x/(log(x)^2-2)

Clase Pr´ actica 5

81

Resoluci´ on del ejemplo 5.5 Se resuelve el ejercicio usando EULER M >>F=inline(’-y./x+y.^2.*log(x)’) F =Inline function: F(x,y) = -y./x+y.^2.*log(x) >>f=inline(’-2./x./(log(x).^2-2)’) f =Inline function: f(x) = -2./x./(log(x).^2-2) >>[es,em]= euler_m(F,f,1,3,1,21) es =1.4013e-001 em =4.6401e-003

Fig. 5.12: Soluciones exactas y aproximadas del ejemplo 5.5 (n = 21, 41) Variando el n´ umero de puntos de la partición se obtienen los resultados siguientes: >>[es,em]=euler_m(F,f,1,3,1,41) es = 7.4861e-002 em = 1.1086e-003 >>[es,em]=euler_m(F,f,1,3,1,81) es = 3.8867e-002 em = 2.7031e-004 >>[es,em]=euler_m(F,f,1,3,1,101) es = 3.1342e-002 em = 1.7211e-004

Ejercicio 5.4 Valorar el comportamiento del método Euler mejorado comparando éste con Euler, mediante el c´ odigo EULER M, aplicado a los siguientes datos: (a) F (x, y) = 2y, f (x) = exp(2x), 0 x d, (b) F (x, y) = y, f (x) = exp( x), 0 x d,

−

−

≤ ≤ ≤ ≤

ambos con y(0) = 1, siendo f (x) la soluci´ on exacta en cada caso, considerando (1) n = 24 , d = 1,

Clase Pr´ actica 5

82 (2) n = 24 , d = 100, (3) n = 27 , d = 100.

Valorar los resultados num´ ericos con el error absoluto, para lo cual el alumno deber´ a utilizar la plantilla donde aparece una tabla de resultados con las entradas ya organizadas. Indicaciones para el Ejercicio 5.4. Separar los problemas seg´ un el apartado al que pertenecen, por ejemplo >>Fa=inline(’2*y+0*x’); >>fa=inline(’exp(2*x)’); >>[es,em]=euler_m(Fa,fa,0,1,1,2^4) % no poner ; al final para producir eco en pantalla,

con lo cual obtenemos directamente y sin pausas los errores generados por el método de Euler y Euler mejorado, es y em respectivamente, para el problema (a) con d = 1 y 2 4 nodos. Definir F b y f b y proceder aplicando EULER M. Las plantillas que el profesor ha entregado en clases deben rellenarse con los resultados para ser entregadas al final. La realización de este tipo de problemas debe permitirnos llegar a conclusiones experimentales y teóricas, y además, ganar experiencia en la aplicación de estos métodos.

5.4

´ ´ OTROS METODOS NUMERICOS

Podemos dise˜ nar otros métodos numéricos para resolver y  = f (x, y), y(x0 ) = y 0 , calculando aproximadamente la integral



xk+1

f (x, y(x))dx,

xk

mediante nuevos procedimientos num´ ericos. Seg´ un sea el orden del error producido por la fórmula de integración que utilicemos, as´ı será el orden del método numérico que permitirá calcular aproximadamente la soluci´ on del PVI.

5.4.1

Ejemplo de m´ etodo de dos pasos

Para desarrollar ideas haremos la deducció n rápida de un tercer enfoque para resolver un PVI. El nuevo método que veremos se clasifica como un método de dos pasos, le llamaremos método del punto medio. Al igual que Euler aproximamos el área bajo la curva sustituy´ endola por el a´rea de un rectángulo, esta vez a la altura y(xk ) de la imagen del punto medio x k del intervalo [xk−1 , xk+1 ] (ver figura 5.13). La formulación del método del punto medio es la siguiente yk+1 = y k−1 + 2hf (xk , yk ).

(5.8)

Clase Pr´ actica 5

83

Fig. 5.13: Ilustración de un método de dos pasos Notar que (5.8) es un método expl´ıcito porque yk+1 aparece despejado, y que la estimación yk+1 depende de las estimaciones yk e yk−1 (dos pasos). Esto último indica que para iniciar el proceso, también debemos conocer de antemano el valor de y1 . Recordar que y0 es un dato del propio problema. No continuaremos tratando este método que el alumno interesado puede encontrar en [1].

5.5

CONCLUSIONES

En general podemos postular que 1. El orden de un método es un indicador que est´ a relacionado con el número de operaciones lógico-aritméticas que han de ejecutarse para alcanzar una precisi´ on dada. Es una referencia vaga ya que en la práctica experimental influyen diversos factores en la eficiencia de un método. No obstante, permite hacer una primera comparación entre dos métodos diferentes. 2. No existe un método numérico que sea de aplicabilidad universal. 3. Los estimados del error a posteriori, ya contemplados en el programa EULER M, permiten juzgar acerca de la calidad de los resultados num´ ericos al aplicar un método a un problema concreto cuya solución exacta es desconocida. 4. El método de Euler es un m´ etodo muy simple, que no se utiliza en la pr´ actica, pero resulta muy adecuado para introducirnos en el estudio general de los métodos numéricos que resuelven aproximadamente problemas de valor inicial. El método de Euler pertenece a una clase de métodos llamados de “un paso”, y que en general se caracterizan por la siguiente relaci´ on recurrente: yk+1 = y k + φ(h, xk , yk ).

Clase Pr´ actica 5

84

5.6

EJERCICIOS ADICIONALES

on de Ejercicio 5.5 Utilizar el programa EULER3 para calcular una soluci´ y =

2x 2 y , 0 9

≤ x < 3, con y(0) = 1.

(5.9)

Detectar una anomal´ıa, en este caso la presencia de una “singularidad” en x = 3, calculando la solución para n = 20, 100, 500, 1000, 1500, 3000, y el correspondiente error (que debe crecer con n). Se sugiere poner punto y coma al final de la l´ınea donde se ha invocado EULER3, para evitar el eco incontrolado en pantalla, y dar únicamente salida al contenido de aquellas variables que deseamos conocer, por ejemplo >>[error,x1,y1]=euler3(F,0,3,1,3000); >>error error =2.2165e+002

La solución exacta del problema es y = 9/(9

5.6.1

2

− x ).

Aplicaci´ on: c´ alculo de la elasticidad b(x)

La elasticidad b(x) es un indicador utilizado por los economistas en algunas de sus versiones discretas para analizar la llamada elasticidad de la demanda y = y(x), como función del precio x de un producto. Este an´ alisis puede conducir a decisiones arriesgadas como es la de subir los precios con expectativas de mayores ingresos, a´ un cuando la demanda disminuye. odigo EULER3 para obtener soluciones aproximadas de Ejercicio 5.6 Utilizar el c´ xy = y b(x), y(1) = 1, en los dominios siguientes (I) 1

≤ x ≤ 2, (II) 1 ≤ x ≤ 2 × 10 , (III) 1 ≤ x ≤ 2 × 10 , 2

12

considerando las elasticidades: 1. b(x) =

−0.2,

2. b(x) = sen (x)

− 1.

Estimar el error global de discretizaci´ on para n = 100 puntos.

(5.10)

Clase Pr´ actica 5

85

Ejercicio 5.7 Resuelva la ecuaci´ on (5.10) mediante variables separables y responda con argumentos a las siguientes preguntas. 1. ¿Es (5.10) una ecuaci´ on diferencial lineal? 2. ¿Influyen el intervalo y el signo de b(x) en la calidad de la aproximaci´ on numérica? 3. ¿Se justifica aplicar un método numérico cuando b(x) es constante (apartado 1)?

´ CLASE PRACTICA 6 M´ etodos num´ ericos para PVI (II)

OBJETIVOS 1. Aplicar los métodos numéricos estudiados en la clase pr´ actica No. 5 a problemas de valores o condiciones iniciales asociados a ecuaciones y sistemas diferenciales. 2. Profundizar en el estudio de las series de Fourier mediante ejemplos gráficos. En esta clase el alumno comprobará la importancia que tiene transformar ecuaciones y sistemas diferenciales de orden superior en sistemas de primer orden. Al igual que la clase anterior, sólo utilizaremos herramientas gráficas y numéricas del sistema MATLAB.

6.1

´ METODOS PARA RESOLVER PROBLEMAS ASOCIADOS A SISTEMAS DE EEDD

Para simplificar consideremos un sistema S de dos ecuaciones y dos incógnitas x = x(t) e y = y(t). S :

 

x = f (t,x,y) y  = g(t,x,y) x(t0 ) = x 0 , y(t0 ) = y 0

Notemos que el anterior sistema S puede ser no-lineal. Si los términos de la derecha no dependiesen de t dir´ıamos que S es un “sistema autónomo”.

6.1.1

M´ etodo de Euler mejorado

Aplicaremos métodos como el de Euler mejorado(1 para estudiar num´ ericamente la soluci´ o n del sistema S . No obstante, cualquiera sea el método que apliquemos a un sistema de EDs se tiene que: 1

Euler mejorado pertenece a la familia de métodos de Runge-Kutta que más adelante estudiaremos.

87

Clase Pr´ actica 6

88

1. En cada etapa debe ser aplicado el método a todas las ecuaciones antes de pasar a la siguiente etapa. Además, el paso h debe ser el mismo para todas las ecuaciones. 2. Otra cuesti´ on que no podemos pasar por alto es que la solución numérica de una ecuación o sistema diferencial de orden mayor que uno debe ser primero transformado en un sistema de orden uno. A éste u ´ltimo se le aplica directamente el método numérico. 3. Cualquiera de los m´ etodos de diferencias finitas que estudiemos para resolver el PVI escalar  y = F (x, y), y(x0 ) = y 0 , puede también ser aplicado a un sistema diferencial. A continuación recordamos cómo se transforma una ED de orden dos y = f (t, y , y)

(6.1)

en un sistema de EEDD de orden uno. Si hacemos x = y  , entonces y  = x  . Sustituyendo en (6.1) y agregando la nueva ecuación y  = x, originada por el propio cambio, obtenemos el sistema de primer orden x = f (t,x,y) y = x De su solución (x(t), y(t)) sólo necesitamos la componente y = y(t) que es la solución del problema original. El algoritmo de Euler mejorado aplicado a un sistema diferencial genérico S queda formulado tal como se muestra a continuación

Euler mejorado aplicado a un sistema xk+1 = x k +

 

h ∗ ) g(tk , xk , yk ) + g(t∗k+1 , x∗k+1 , yk+1 2 donde las predicciones para k + 1 están dadas por yk+1 = y k +

 

h ∗ ) f (tk , xk , yk ) + f (t∗k+1 , x∗k+1 , yk+1 2

(6.2)

t∗k+1 = t k + h x∗k+1 = x k + hf (tk , xk , yk )

(6.3)

∗ = y + hg(t , x , y ) yk+1 k k k k

6.1.2

Sobre los programas S2 EULER y RK CELL

El programa S2 EULER, descargable desde la Web de la asignatura [15], se compone de dos módulos. El primero, a modo de interfaz logra que el código sea transparente al usuario. El alumno sólo tiene que saber qué datos debe dar como entrada, e interpretar los resultados.


89

Los datos relativos al campo de direcciones o a la entrada opcional de la solución on exacta se entregan cuando lo solicita el programa, con la sintaxis o notación on usual, sin encerrar entre apóstrofes ostrofes ni poner puntos delante de operadores. Esto ultimo último es necesario necesario para la ejecuci´ ejecuci´ on pero el propio programa se on encarga de hacerlo automáticamente. aticamente. Como sintaxis se admite una cualquiera de las siguientes:(2 s2_e s2_eul uler er

y=s2 y=s2_e _eul uler er

[t,y [t,y]= ]=s2 s2_e _eul uler er

[t,x [t,x,y ,y]= ]=s2 s2_e _eul uler er

El segundo módulo odulo está constituido por dos sub-programas situados al final del c´ odigo. odigo. Uno de ellos, tal como se señal´ naló anteriormente, se encarga de “escribir” el punto delante del operador división on \, la multiplicaci´ on on *, y la exponenciación on ^. La otra sub-función, on, de mayor interés es para p ara el alumno, a lumno, permite simular el método etodo de Euler mejorado aplicado a un sistema diferencial de dos ecuaciones, y se corresponde con la formulación on (6.2-6.3). Los cálculos alculos se realizan atendiendo al paso h paso h dado por el usuario, y para h/ para h/2, 2, con la finalidad de comparar los correspondientes resultados y obtener un estimado del error a posteriori. Pueden entrarse varios pasos h, de modo que al final se dispone de una sola vez de los errores correspondientes. Los gráficos aficos de las soluciones aproximadas aparecen en pantalla superpuestos y sólo olo corresponden al ultimo u ´ltimo paso. El programa programa RK CELL, que tambi´ también en podemos p odemos obtener en la Web de la asignatura asignatura [15], es una interfaz de usuario para facilitar el uso de la herramienta ODE45 de MATLAB aplicada a sistemas de ecuaciones. El dise˜ no de RK CELL se basa en usar objetos CELL para simular funciones vectoriales no de cualquier número umero de dimensione dimensiones. s. El propio programa programa crea autom´ aticamente el fichero FUNCaticamente TION correspondiente al campo vectorial, con las caracter´ caracter´ısticas exigidas por ODE45. Este comando Matlab ejecuta algoritmos correspondientes corresp ondientes a los métodos etodos de Runge-Kutta que estudiaremos en la clase práctica actica 7.

6.1.3 6.1.3

Ejerc Ejercici icios os con indic indicaci acion ones es para para ser resuel resuelto toss

Ejercicio 6.1 Dado el sistema diferencial lineal siguiente = x + y, x = x  = x y, y = x x(0) = 1 + 1, y(0) = 1,

− √ 2,

(6.4)

k

donde 0 t resolverlo.

05/2 , k = 0, 1, 2, 3, 4, para ≤ ≤ 1, emplear el programa S2 EULER con los pasos h = 0.05/

Indicaciones para el ejercicio 6.1 Escribir 2

Para más as información on ejecutar

help s2 euler


90 >>s2_euler

y oprimir INTRO. La solución on exacta del sistema (6.4), que es una entrada opcional, es la siguiente x(t) = (1 +

√ √ 2)e 2)e , t 2

y (t) = e t

√

2

El problema problema (6.4) es lineal lineal y homog´ homogéneo. eneo. El alumno deber´ deber´ıa intent intentar ar resolver resolverlo lo artesanalme artesanalmente nte mediante los métodos etodos estudiados en clases. Ver la figura 6.1.

Fig. 6.1: Gráfico afico producido por S2 EULER al resolver el ejercicio 6.1

ericamente el sistema s istema no lineal: Ejercicio 6.2 Resolver numéricamente = x 2 + x + 5y 5y x = x  2 = y + sen (x) + 3x 3x y = y 0.001, 001, x(0) = 0. 0.001, 001, y(0) = 0.

−y

0

≤ t ≤ 1.

Obtener una estimaci´ on del error para h = 0.05/ 05/2k , k = 0, 1, 2, 3, 4.

Indicaciones para el ejercicio 6.2 Escribimos de nuevo >>s2_euler

oprimimos INTRO e introducimos lo datos del problema para después es obtener los resultados. No cerramos el gráfico afico de salida, cuyo aspecto podemos apreciar anticipadamente en la figura 6.2, para poder comparar con el que produciremos en el ejercicio 6.3.


91

Fig. 6.2: Gráfico afico producido por S2 EULER al resolver el ejercicio 6.2

Ejercicio 6.3 Resolver Resolve r num´ n uméricamente ericamente el sistema sistem a lineal: li neal: = x + 5y 5y x = x  y = 4x y 0.001, 001, x(0) = 0. 0.001, 001, y(0) = 0.

−

0

≤ t ≤ 1.

Obtener una estimaci´ on del error para h = 0.05/ 05/2k , k = 0, 1, 2, 3, 4.

Indicaciones para el ejercicio 6.3 Antes de ejecutar >>s2_euler aplicar HOLD ON para que el gráfico afico anterior, del ejercicio 6.2, quede superpuesto al gráfico afico del ejercicio 6.3.

A continuaci´ on on se muestra c´ omo se ha “perturbado localmente” el problema 6.2 para omo transformarle en el problema 6.3 Los problemas planteados en los ejercicios 6.2 y 6.3 están an relacionados entre s´ s´ı. El estudiante debe apreciar que las ecuaciones diferenciales del ejercicio 6.2 son no lineales , mientras que las del 6.3 son todas lineales. lineales. Mostraremo Mostraremoss cómo omo el campo de direcciones del sistema diferencial del ejercicio 6.2 se aproxima localmente mediante su diferencial, y cómo omo a partir de éste este ultimo u ´ ltimo se producen las ecuaci ecuacione oness lineal lineales es del problema problema 6.3. Se dice dice que el sistem sistemaa del ejercic ejercicio io 6.3 es una aproxima aproximaci´ ci´ on on lineal del sistema sistema del ejercicio ejercicio 6.2. Las expectativ expectativas as son: que la soluci´ on de 6.3 sea una aproximación on on aceptable de la correspondiente a 6.2 en un intervalo [0, [0 , ε], con ε > 0 “suficientemente pequeño”, no”, y que el método etodo numérico erico seleccionado seleccion ado sea más as efectivo resolviendo el problema lineal 6.3, que el (3 problema no lineal 6.2. 3

La frase “suficientemente pequeño” no” encierra otro problema cuya solución no abordaremos en este material.

Clase Pr´ actica 6

92

Fig. 6.3: Gráfico producido por S2 EULER al resolver el ejercicio 6.3 El cálculo de los siguientes diferenciales en (0, 0) D(x2 + x + 5y)(0, 0) = dx + 5dy D(y2 + sen(x) + 3x

− y)(0, 0) = 4dx − dy

y la conocida fórmula f (x, y)

≈ f (0, 0) + Df (0, 0)(x − 0, y − 0)

válida para x e y “peque˜ nos”, nos muestran que x2 + x + 5y y2 +

≈ x + 5y sen(x) + 3x − y ≈ 4x − y

Es oportuno se˜ nalar que los valores iniciales x(0) = 0.001, y(0) = 0.001 se han seleccionado intencionadamente “peque˜ nos” en correspondencia con el proceso anterior.

Ejemplo 6.1 El siguiente PVI, extra´ıdo del Bolet´ın 7, consiste de una ED lineal y homogénea de orden dos con coeficientes variables. y ty + y = 0, y(0) = 2, y  (0) = 1, t

−

∈ [0, 1],

(6.5)

Indicaciones para el ejemplo 6.1 El PVI (6.5) del ejemplo 6.1, puede ser resuelto mediante el método de las series de potencias. Ahora nos proponemos aplicar el programa S2 EULER al sistema de orden uno que se le asocia, y que el propio alumno debe calcular directamente.

Clase Pr´ actica 6

93

Fig. 6.4: Gráfico obtenido al resolver el ejemplo 6.1 Como comprobaci´ on de los resultados apreciaremos que la gr´ afica de la soluci´ on numérica, obtenida seg´ un Euler, es muy similar a la del polinomio y(t) = 2 + t

2

t4 12

− t − −

t6 , 120

(6.6)

obtenido al truncar la serie de potencias que da la soluci´ on exacta de este problema. (4

Resoluci´ on del ejemplo 6.1 A continuaci´ on expondremos el proceso de obtenció n de la solució n exacta y los gráficos de las soluciones aproximadas. Comencemos viendo que la solución exacta del problema (6.5), expresada en términos de una serie de potencias en torno al punto x = 0, es

 ∞

an tn ,

n=0

siendo los a n definidos por la recurrencia a tres términos an+2 =

an (n 1) , n = 0, 1, 2,... (n + 2)(n + 1)

−

donde a 0 = 2 y a 1 = 1, seg´ un se indica en las condiciones iniciales. 4

Observar que siendo t [0, 1], y atendiendo a la precisión deseada, el término con la potencia t6 también podr´ıa despreciarse en la aproximación (6.6).

∈

Clase Pr´ actica 6

94

Sea P=[a_6,...,a_0], el vector formado por los coeficientes an , n = 0,..., 6, de la suma parcial que estamos considerando como una aproximación de la solución del problema (6.5) 6



an tn = a 6 t6 +

n=0

··· + a t + a . 1

0

El comando POLYVAL interpreta a P como un polinomio y permite evaluarlo en variables escalares o vectoriales.(5 Para comparar los resultados obtenidos mediante series de potencias y el método Euler mejorado aplicado al sistema asociado, ejecutar S2 EULER optando por la salida gráfica y hacer >>hold on >>tt=linspace(0,1,50); >>P=[-1/120, 0, -1/12, 0, -1, 1, 2]; >>yy=polyval(P,tt); >>plot(tt,yy,’+g’)

El resultado en pantalla deber´ıa ser similar a la figura 6.4.

RESUMEN Para resolver una ED o un sistema de EEDD de orden mayor que uno, éstos deben ser convertidos previamente en un sistema de orden uno, y aplicarle a este último uno de los métodos estudiados para el caso escalar. La serie de potencias que representa a la solución de un problema con datos anal´ıticos, puede ser simulada numéricamente utilizando una suma parcial seleccionada convenientemente, lo que, en principio, constituye otro método para aproximar la soluci´ on.

6.2

EJERCICIOS ADICIONALES

on del PVI Ejercicio 6.4 Aproximar los valores de la soluci´ y + 2t2 y + 5y = 0, y(0) = 1, y (0) = 1, t [0, 2].

−

∈

en los puntos t = 0.21, 0.95, 1.31, 1.84. ognita y = y(t) Ejercicio 6.5 Aproximar los valores de la inc´ x = tx + 5y + t2 , y = 2x + 6y + 4,

−

x(0) = 0, y(0) =

−1, t ∈ [0, 1].

en los puntos t = 0.11, 0.23, 0.59, 0.82. 5

El comando POLYVAL es compatible con la arquitectura matricial del MatLab, es decir, admite matrices como argumento de entrada.

Clase Pr´ actica 6 Ejercicio 6.6

95

x = tx 10y 4t 2t2 , y = x + 6y, x(0) = 1, y(0) = 1, t [0, 1].

−

− − − ∈

Ejercicio 6.7 y = x + 6y + 3, x = 2x + y, x(0) = 2, y(0) = 4, t

6.2.1

∈ [0, 2].

Indicaciones para resolver los ejercicios (6.4)-(6.7)

Indicaciones para el ejercicio 6.4 Si intentamos resolver el ejercicio 6.4 con DSOLVE obtenemos >> dsolve(’D2y+2*t^2*Dy+5*y=0’,’y(0)=1’,’Dy(0)=-1’) Warning: Explicit solution could not be found. > In C:\MATLABR11\toolbox\symbolic\dsolve.m at line 326 ans = [ empty sym ]

Abordemos la solución por la v´ıa numérica. Para ello escribimos la EDLH de 2do orden del ejercicio 6.4 como un sistema de primer orden x = 2t2 x 5y, y = x,

−

x(0) =

−

−1, y(0) = 1,

para aplicar S2 EULER.M o RK CELL.M. Hagamos >>[t,y]=s2_euler

porque sólo nos interesa almacenar en memoria los valores aproximados de y = y(t), la incógnita del problema original. Entrar h = 0.2 Desplegamos el gráfico para contrastarlo con el de la figura 6.5. Los puntos t donde se nos pide el valor aproximado de la solución no está previsto que coincidan con los de la red uniforme asociada al paso h escogido, por tanto, debemos interpolar. Ejecutamos el programa INTERPOLA con la siguiente sintaxis >>[T1,Y1]=interpola(t,y)

Entramos los puntos donde queremos interpolar seg´ un se nos indica en la plantilla. Entrar el vector T1=[t1,...,tm]=[.21 .95 1.31 1.84]

A continuación aparecen las listas de los puntos interpolados (T 1 , Y 1 ) y el gráfico que muestra su posición (ver figura 6.5). La leyenda que aparece en la esquina superior izquierda puede moverse arrastrándola con el cursor.

Clase Pr´ actica 6

96

Fig. 6.5: Gráficas de S2 EULER.M e INTERPOLA.M. Ejercicio 6.4

Indicaciones para el ejercicio 6.5 Para el ejercicio 6.5 tampoco DSOLVE ofrece respuesta. Es posible que el término tx de la primera ecuaci´ on provoque una complejidad excesiva o no existe, de hecho, una solución representable con el catálogo de funciones MATLAB. Nos queda la alternativa de aplicar S2 EULER o RK CELL al sistema anterior. Los resultados gráficos que se muestran en la figura 6.6 corresponden al método RK CELL. Proceder tal como se hizo en el ejercicio 6.4, empleando el programa INTERPOLA (ver fig. 6.6).

Fig. 6.6: RK CELL e INTERPOLA.M aplicados al ejercicio 6.5

Indicaciones para el ejercicio 6.6 Aplicamos S2 EULER con h = 0.001 al ejercicio 6.6, y ejecutamos en la l´ınea de comandos


97

>>hold >>hold on

Podemos comprobar experimentalmen exp erimentalmente te la coincidencia de dos métodos etodos vistos en clases. Si transformamos el sistema del problema 6.6 en una EDL de 2do orden obtenemos 6)y + (6t (6 t + 10)y 10)y = y (t + 6)y  y(0) = 1, y (0) = 5.

−

−

−

2

−4t − 2t ,

La anterior ecuación on es excesivamente complicada para el n´ ucleo ucleo simbólico olico de MATLAB. La ecuación on homogénea enea asociada aso ciada tiene solución on mediante el comando DSOLVE en términos ermi nos de la l a funci f unci´ón on especial (6 WhittakerM. Notemos que la ED lineal y

2

− (t + 6)y 6)y + (6t (6 t + 10)y 10)y = −4t − 2t ,

puede ser resuelta mediante series de potencias aunque su proceso de resolución es algo más as complicad plicadoo que los proble problemas mas vistos vistos en clases. clases. Sin embarg embargo, o, lo aprend aprendido ido nos sobra sobra para abordar abordar su solución on en torno a t a t = = 0. Después es de asumir que y que y((t) es anal ana l´ıtica ıtic a en x en x 0 = 0, llegamos a que

y(t) =

 ∞

an tn

n=0

donde

Fig. 6.7: Solución on gráfica afica del ejercicio 6.6 6

Usar mhelp en lugar de help para ver más as detalles sobre esta función on especial.


98

Fig. 6.8: Soluciones DSOLVE y RK CELL. Ejercicio 6.7 a0 a1 a2 a3 a4

= 1, = 5, = 3a1 5a0 = = (12a2 9a1 = (18a3 8a2

··· an+2 =

− −

···

−

− −

10, −10, − 6a − 4)/ 4)/6, 2)/12, 12, − 6a − 2)/ + (n ( n − 10)a 10)a − 6a − , 0 1

6(n 6(n + 1)a 1)an+1 n (n + 2)(n 2)(n + 1)

n 1

n = 3, 4, 5,...

Ahora ejecutamos el programa TRUNCADA.M- que descargamos desde [15]- para calcular las sumas parciales de la anterior serie de potencias. Para obtener la suma parcial correspondiente al término N ermino N = = 12, procedemos como sigue. >>t=linspace(0,1,100); >>yy=truncada(t) >>hold >>hold on >>plot(t,yy,’+’)

En la figura 6.7 se comprueba gráficamente aficamente el nivel de coincidencia entre las dos soluciones aproximadas.

Indicaciones para el ejercicio 6.7 El ejercicio 6.7 tiene solución on calculable con DSOLVE pero con una expresión on excesivamente extensa.(7 Los gr´ aficos de las soluciones exactas x(t) e y (t) pueden verse en la figura 6.8 y compararse con aficos las obtenidas usando S2 EULER o RK CELL. Para hacer los gráficos aficos se recomienda usar EZPLOT. 7

Con otra versión on de MATLAB o con otro software el resultado puede variar.


99

>>S=dsolve(’Dy=x+6*y+3’,’Dx=2*x+y’,’x(0)=2’,’y(0)=4’); % En S.x y S.y estan las soluciones x(t) e y(t). >>ezplot >>ezplot(S.x (S.x,[0, ,[0,2])% 2])%camb cambiar iar el color color de la grafica grafica >>hold >>hold on >>ezplot(S.y,[0,2]) >>title( >>title(’ ’ ’)% titulo titulo vacio vacio >>legend(’S.x’,’S.y’)

6.3 6.3 6.3. 6.3.1 1

SERI SERIES ES DE FOURI OURIER ER Prog Progra rama mass

Los ejercicios de esta sección on se basan en la utilización on de los programas SERIE, SERIE1 y SERIE2, todos descargables desde la Web de la asignatura [15]. El programa programa SERIE nos permite dibujar los gráficos aficos superpuestos de las 6 primeras sumas parciales de la serie de Fourier trigonométrica etrica de la función on (6.7) (figura 6.9). Mediante SERIE1 podemos dibujar los gráficos aficos no superpuestos de las 20 primeras sumas parciales de la serie de Fourier trigonométrica etrica de la función on (6.7) (figura 6.9). Con SERIE2 podemos visualizar la convergencia de las series de Fourier viendo los gráficos aficos sobre ( 3π, 3π ), de las primeras 25 sumas parciales de la función f on f ((x) = x para x para π < x < π. π.

−

6.3. 6.3.2 2

−

Ejer Ejerci cici cios os

on Ejercicio 6.8 Dada la funci´ f ( f (x) =

 −

1 si 1 si

0 , −π ≤ x < 0, 0≤x<π

(6.7)

representar gr´ aficamente la funci´ on f y las apr aproximaci oximaciones ones de su serie de Fourier trigonom´ trigonom´ etrica etrica mediante las primeras sumas parciales.

Indicaciones para resolver el ejercicio 6.8 La serie de Fourier de f viene f viene dada por 4 π

 ∞

n=1

sen((2n sen((2n 1)x 1)x) . 2n 1

− −

Basta ejecutar los programas programas SERIE y SERIE1. on Ejercicio 6.9 Representar la funci´ f ( f (x) = x en [ π, π] y su extensi´ on 2π-peri´ odica al intervalo [ 3π, 3π].

−

−

Clase Pr´ actica 6

100

Fig. 6.9: Ejercicio 6.8. Sumas parciales

Fig. 6.10: Ejercicio 6.9. Sumas parciales N = 2, 26

Indicaciones para resolver el ejercicio 6.9 La serie de Fourier de f viene dada por

2

− ∞

n=1

( 1)n+1

sen(nx) . n

Usar el programa SERIE2 para realizar esta tarea (ver figura 6.10).

Clase Pr´ actica 6

6.4 6.4.1

101

´ APENDICE Programa SN EULER

A modo de complemento, el alumno puede utilizar el programa SN EULER para resolver sistemas de ecuaciones diferenciales con un número arbitrario de ecuaciones, utilizando el método Euler mejorado. Este código consta de tres partes: la primera es una interfaz de usuario, y las dos restantes son subfunciones situadas al final del script. Una de las sub-funciones es de especial inter´ es para el alumno pues se corresponde con el procedimiento de Euler mejorado expresado en términos vectoriales. La otra subrutina FUNCTION, llamada CEVAL, permite la evaluación numérica de campos vectoriales almacenados como objetos CELL. Este programa puede ser descargado de la Web de la asignatura [15], y se activa directamente en la l´ınea de comandos haciendo >>sn_euler

En caso de que deseemos almacenar en memoria los valores aproximados de las incógnitas escribimos >>x=sn_euler;

Si además queremos almacenar en la variable t los valores de la red t0 ,...,tm, hacemos >>[x,t]=sn_euler; La variable x es una matriz con tantas columnas como ecuaciones tiene el sistema. El número de filas de x se corresponde con el número de puntos t 0 ,...,tm .

Para visualizar una breve descripción del programa, en particular lo concerniente a la entradasalida, ejecutamos en l´ınea lo siguiente >>help sn_euler

Para obtener la respuesta ................... ................... ENTRADA Numero de ecuaciones Componentes del campo vectorial Intervalo dominio Condicion inicial paso de integracion SALIDA Error a posteriori Graficos de las soluciones aproximadas SALIDAS OPCIONALES -Para almacenar en memoria la solucion [x1(1),...x1(m); x2(1),...,x2(m);... ;xn(1)...xn(m)] escribir >>x=sn_euler;

Clase Pr´ actica 6

102 -Para almacenar en memoria la solucion [x1(1),...x1(m); x2(1),...,x2(m);... ;xn(1)...xn(m)] y ademas: [t(1);t(2);....;t(m)] escribir >>[x,t]=sn_euler;

6.4.2

Interpolaci´ on

Sea el PVI: y  = F (x, y), y(x0 ) = y 0 . Una vez calculados los valores aproximados de la solución y(x) de un PVI, sean estos

≈ y(x ),

yk

k

k = 1,...,n

nos puede interesar aproximar el valor de y(t) para t (xk−1 , xk ). Seguiremos el procedimiento más sencillo consistente en “interpolar linealmente”, es decir uniremos los puntos P k−1 = (xk−1 , yk−1 ), P k = (xk , yk )

∈

mediante el segmento de la recta y = mx + b que dichos puntos determinan. El valor aproximado es mt + b = y ∗ (t)

≈ y(t).

Se trata pues, de construir un aproximante poligonal- lineal a tramos- continuo, de una funci´ on y(x) que admite derivadas hasta el segundo orden (ver figura 6.11). Existen otros métodos de interpolaci´ on que aseguran la suavidad del aproximante, pero no los trataremos en este curso. La figura 6.11 ilustra el proceso de interpolación mediante poligonales. Para realizar eficientemente los cálculos usaremos el programa INTERPOLA.M -descargable desde [15]- basado en el comando INTERP1 de MATLAB. >>Y1=interp1(T,Y,T1);

donde T = [x(1)

··· x(n)], Y = [y(1) · ·· y(n)], T 1 = [t(1) ··· t(n)], Y 1 = [m(1) ∗ t(1) + b(1) ··· m(n) ∗ t(n) + b(n)] Notar que hasta ahora los métodos de diferencias finitas nos han ofrecido la soluci´ on como un vector n de R , y que según el procedimiento basado en interpolación estamos hallando una solución en la forma de una función que, aunque continua, desafortunadamente no es derivable en todos los puntos. A continuación se muestra el código de INTERPOLA.

Clase Pr´ actica 6

103

Fig. 6.11: Ilustración gráfica del método de interpolación function interpola disp(’Se exige que x1
´ CLASE PRACTICA 7 M´ etodos num´ ericos para PVI (III)

OBJETIVOS Adquisici´ on de conocimientos teórico-experimentales sobre los métodos numéricos de resoluci´ o n de problemas de condici´ on inicial. Visualizar las soluciones obtenidas mediante separación de variables para las ecuaciones del calor y de ondas.

7.1

´ METODOS DE RUNGE-KUTTA

7.1.1

M´ etodos de Taylor

Sea la ED y  = f (x, y), sujeta a la CI y(x0 ) = y 0 . Los métodos de Euler y Euler mejorado son de órdenes 1 y 2, respectivamente, y sólo requieren evaluaciones de la función f (x, y). A continuación presentaremos brevemente a los Métodos de Taylor, antes de pasar a los métodos de Runge-Kutta. etodos de Taylor consiguen órdenes altos para el error con el correspondiente Los llamados m´ coste que supone tener que derivar repetidamente la función. Para ello se parte del desarrollo de Taylor de la solución y(x) en x0 y(x0 + h) = y(x0 ) + hy (x0 ) + h2 /2!y (x0 ) +

n

··· + h /n!y

(n)

(x0 ) +

(n+1)

O(h

),

y se sustituyen las derivadas hasta el orden n, y (x0 ) = f (x0 , y0 ) y (x0 ) = f x (x0 , y0 ) + f y (x0 , y0 )f (x0 , y0 ) =

···

· ··

···

El método de Taylor de orden 1 es el ya conocido método de Euler expl´ıcito, estudiado en la anterior clase práctica, y en el mismo no intervienen derivadas. Los métodos de Taylor de orden n 2, que no estudiaremos, s´ı involucran a las derivadas de f (x, y).

≥

105

Clase Pr´ actica 7

106

Sustituyendo las derivadas de y(x) en x0 hasta un orden n 3 arbitrario, podemos obtener un método de mayor orden que el de Euler mejorado. Estos métodos se llaman métodos de Taylor debido al papel jugado por el desarrollo de Taylor de f (x, y) en su deducción. Aunque, tal como se ha mencionado, se pueden teóricamente alcanzar órdenes altos para el error, éstos dependen de que podamos derivar a f (x, y) repetidamente, lo cual en la práctica no siempre es factible, es impracticable con herramientas num´ ericas y muy costoso con simb´ olicas.(1 El camino que se sigue para evitar el inconveniente de calcular derivadas de diferentes órdenes, un proceso en el que ocurren graves problemas num´ ericos, es el diseño de los llamados métodos de Runge-Kutta. Estos combinan ingeniosamente los valores de f (x, y) para evitar la presencia de las derivadas en los cálculos y al mismo tiempo, producir o´rdenes similares a los que se obtienen con los métodos de Taylor.

≥

7.1.2

M´ etodo de Runge-Kutta de cuarto orden

Hay una variedad de métodos de Runge-Kutta que difieren entre s´ı atendiendo a la forma en que combinamos los valores de la función f (x, y) que define el PVI. El que sólo requieran evaluaciones de la función f (x, y) es una cualidad que los Runge-Kutta comparten con Euler y Euler mejorado, pero llegando alcanzar mayores órdenes que estos u ´ ltimos. Los métodos de Runge-Kutta están considerados entre los más adecuados para resolver problemas de valores iniciales. Los de orden n = 2 incluyen entre otros al método de Euler mejorado. El de orden 4 clásico es una buena elección por su precisió n y fácil programación. El esquema del método clásico de Runge-Kutta de 4to orden es el siguiente: yk+1 = y k + (k1 + 2k2 + 2k3 + k4 )/6 donde

7.1.3

k1 = hf (xk , yk ) k2 = hf (xk + h/2, yk + k1 /2) k3 = hf (xk + h/2, yk + k2 /2) k4 = hf (xk + h, yk + k3 )

Descripci´ on del programa RUNGEKUTTA

Se trata de un programa doméstico que podemos descargar desde la Web de la asignatura [15], con el cual intentaremos resolver EEDD de primer orden y  = f (x, y) con una condición inicial y(a) = y 0 , x 1

∈ [a, b].

La derivación numérica es un problema mal acondicionado, lo que significa que pequeños errores en los datos pueden producir grandes errores en los resultados.

Clase Pr´ actica 7

107

Tal como indica su nombre, este código permite simular un método de Runge-Kutta, particularmente el clásico de orden 4. El programa cuenta con una interfaz de usuario desarrollada para facilitar el trabajo, que permite opcionalmente introducir la solución exacta si fuese conocida. En una sola ejecuci´ o n es capaz de calcular los errores a posteriori para un vector de pasos [h1 ,...,hm ]. También muestra los gráficos de las aproximaciones para hm y hm /2. El error a posteriori lo calcula como la mayor de las discrepancias entre las soluciones para h y h/2. Para invocar a RUNGEKUTTA se oprime INTRO después de escribir en la l´ınea de comandos alguna de las siguientes sentencias >>rungekutta

>>Y=rungekutta

>>[T,Y]=rungekutta

dependiendo de la información que deseamos almacenar en memoria. Aqu´ı T representa al vector cuyas componentes ordenadas forman la red de puntos T (k) del dominio seg´ un el paso h, e Y es el vector de las correspondientes aproximaciones Y (k) de la solución. En cualquier caso el programa permite opcionalmente ver en la pantalla los valores T (k) e Y (k) (hasta 50 como máximo), pero no los almacena en memoria.

7.1.4

Ejercicios

on del PVI Ejercicio 7.1 Aproximar la soluci´ y = 1 + y2 , y(0) = 0 para x

∈ [0, 1], con pasos h 0.05 0.05/2 0.05/4 0.05/8 0.05/16

utilizando el programa RUNGEKUTTA. Teniendo en cuenta que la soluci´ on exacta es y = tan(x), calcular el error exacto cometido y compararlo con el estimado a posteriori. on del PVI Ejercicio 7.2 Aproximar la soluci´ 2y = (x y(0) = 1 para x

− y),

∈ [0, 3], con pasos h 0.05 0.05/2 0.05/4 0.05/8 0.05/16

utilizando el programa RUNGEKUTTA. Teniendo en cuenta que la soluci´ on exacta es y = x y compararlo con el estimado a posteriori.

− 2 + 3e−

x/2

calcular el error exacto cometido

Clase Pr´ actica 7

108

Ejercicio 7.3 Aproximar la soluci´ on del PVI (ex + 1)y  + y(ex y(0) = 3 para x

− 1) = 0,

∈ [0, 1], con pasos h 0.05 0.05/2 0.05/4 0.05/8 0.05/16

utilizando el programa RUNGEKUTTA

12ex Teniendo en cuenta que la soluci´ on exacta es y = x calcular el error exacto cometido y (e + 1)2 compararlo con el estimado a posteriori.

7.1.5

Estabilidad de un m´ etodo num´ erico

Un método numérico para resolver una EDO de primer orden, está dado por una ecuación en diferencias finitas (EDF). Decimos que el método es estable si su EDF produce una solución acotada cuando la solución exacta es acotada. Si la soluci´ on obtenida no es acotada, siendo la exacta acotada, decimos que el método es inestable. Sólo analizamos la estabilidad cuando tratamos con EEDD ordinarias lineales y  + py = q (x), p R, en las que el término q (x) no influye. Los siguientes enunciados pueden verse en Hoffman [3], secciones 7.8-7.9.

∈

• Los métodos de Euler y Euler mejorado son estables cuando 0 < ph ≤ 2. • Runge-Kutta es estable cuando 0 < ph ≤ 2.784. on del siguiente problema de valor inicial (Hoffman Ejercicio 7.4 Estudiar numéricamente la soluci´ [3], 7.14) aplicando el programa EJE7.4 (Web de la asignatura [15]). y  = p(y (x + 2)) + 1, y(0) = 1, x [0, + ).

∈

Considere los valores p = indica en la tabla 7.1.5.

− ∞

±1000 considerando los intervalos [0, b], b > 0, y los pasos h seg´ un se

a) ¿Influyen los valores de p, b y h sobre el error y la estabilidad? b) Tratar de explicar el papel que desempe˜ na p en los resultados calculando artesanalmente la soluci´ on de la homogénea asociada. Seg´ un se aprecia en la tabla 7.1.5, el tamaño de h se selecciona más peque˜ no para mostrar: (A-B) c´ omo cambia la precisión ( p = 1000)

Clase Pr´ actica 7

109

b A 0.01 B 0.01 C 0.01 D 1 E 1

Euler RK-2 RK-4 ODE45 h p 1e-03 1000 5e-04 1000 1e-03 -1000 1e-03 -1000 5e-03 -1000

Tabla 7.1: Informe de resultados según los valores de b , h y p

(D-E) estabilidad condicional ( p =

−1000)

(Se detectan en los valores numéricos del error) (ODE45 no se afecta pues selecciona automáticamente el paso h) Los diferentes valores de b revelan: (C-D) la existencia de dos escalas relevantes para la variable independiente ( p =

−1000)

(Se detecta gráficamente con los cuatro programas) (En este ejemplo de PVI, el error producido por ODE45 se modifica algo con el cambio de intervalo) Notar que hB =

7.1.6

hA hE y que h D = . 2 5

Descripci´ on del programa EJE7.4

Tal como sugiere su nombre, este código está especialmente diseñado para resolver el ejercicio 7.4 de esta práctica. Para facilitar el trabajo se cuenta con una interfaz de usuario que permite seleccionar uno de los dos posibles PVI mediante una ventana-menú, atendiendo a que la ecuación que le define depende del valor de un parámetro p. El programa se encarga de aplicar los métodos de Euler, Euler mejorado y Runge-Kutta 4, mostrando al final los respectivos errores y los gráficos (ver figura 7.1). Se aplica adicionalmente la herramienta MATLAB ODE45, basada en Runge-Kutta 4 y 5. Esta u ´ ltima adopta como tolerancia de entrada (error máximo permitido), al menor de los errores producidos por los métodos anteriores. Los principales objetivos de este experimento consisten en la detección experimental de inestabilidad y escalas relevantes en el dominio.

Clase Pr´ actica 7

110

Fig. 7.1: Gráficas producidas por EJE7.4.M

7.2 7.2.1

PROBLEMAS DE CONTORNO Planteamiento matem´ atico

Sea el operador diferencial lineal de segundo orden Λ[y] = a(x)y + b(x)y  + c(x)y, siendo a(x) > 0, para todo x [a, b], y consideremos los operadores lineales de primer orden

∈

B1 [y](x) = α 1 y(x) + α2 y (x), B2 [y](x) = β 1 y(x) + β 2 y (x). Al siguiente problema Λ[y](x) = f (x), x B1 [y](a) = 0, B2 [y](b) = 0,

∈ [a, b],

(7.1)

Clase Pr´ actica 7

111

se le llama problema de contorno con condiciones separadas. El problema de valores propios se establece como el problema de hallar los λ tales que Λ[y] = λy, B1 [y](a) = 0, B2 [y](b) = 0,

−

(7.2)

admite soluciones no triviales y = y(x). Sea w(x) el factor de simetrización asociado a Λ[y], es decir L[y] = w(x)Λ[y] = ( p(x)y ) + q (x)y, entoces (7.2) se convierte en el Problema Regular de Sturm-Liouville L[y](x) + λw(x)y = 0, B1 [y](a) = 0, B2 [y](b) = 0.

7.2.2

(7.3)

Aplicaci´ on: Viga columna

Para calcular las cargas cr´ıticas que actúan longitudinalmente sobre una columna empotrada y las deformaciones correspondientes disponemos del modelo matemático siguiente λ y (x) + y(x) = 0 EI (7.4) y(0) = y(L) = 0, donde EI es el momento flector y L es la longitud de la viga. En la Web de la asignatura, encontramos el programa PANDEO que utiliza las soluciones del sistema de Sturm-Liouville formado por λ y  (x) + y(x) = 0 y diferentes condiciones homogéneas de frontera, para mostrar gráficamente las EI deformaciones (autofunciones) y calcular las cargas que producen deformación (autovalores).

Ejemplo 7.1 Ejecutar en la l´ınea de comandos >>pandeo

e introducir los datos seg´ un se soliciten.

7.3 7.3.1

ECUACIONES EN DERIVADAS PARCIALES Ecuaci´ on del calor en una barra de longitud finita

aficamente la soluci´ on del siguiente problema: Ejemplo 7.2 Representar gr´ ut = 5uxx , 0 x 3, t 0, u(0, t) = u(3, t) = 0, u(x, 0) = 5 sen (4πx) 3 sen (8πx) + 2 sen (10πx).

≤ ≤

≥

−

(7.5)

Clase Pr´ actica 7

112

Fig. 7.2: Gráfica producida por CALOR.M

Resoluci´ on del ejemplo 7.2 La solución exacta u(x, t) de (7.5) es u(x, t) = 5 sen(4πx)exp( 5(122 )π2 t/9)

−

2

2

2

2

− 3 sen(8πx)exp(−5(24 )π t/9) + 2 sen(10πx) exp(−5(30 )π t/9),

y se ha obtenido aplicando el método de separaci´ on de variables. El programa CALOR (se descarga en [15]) utiliza la solución de (7.5) para hacer los gráficos del ejemplo 7.2. on de: Ejemplo 7.3 El flujo de calor a través de una barra cil´ındrica de longitud π, es soluci´ ut = 2uxx ; 0 < x < π, t > 0, u(0, t) = u(π, t) = 0, u(x, 0) =



x π

−x

0 < x π/2, π/2 < x < π.

≤

(7.6)

1. Obtener una aproximaci´ on de la soluci´ on (serie (7.8)) en el instante t = 0, en x = π/2.

(7.7)

Clase Pr´ actica 7

113

2. Comparar la aproximaci´ on con el dato inicial u(x, 0) = f (x) 3. Representar gr´ aficamente la temperatura en el instante t = 1.25 en toda la barra. 4. Ver la evoluci´ on de la temperatura en dos instantes de tiempo. Para ello representar la distribuci´ on de temperaturas para t = 1.25 y t = 2. (Notar el descenso de temperaturas conforme avanza el tiempo.)

Recursos para resolver el ejemplo 7.3 El método de separación de variables permite obtener la solución como la suma de la serie ∞ 2 4 k exp( 2(2k + 1) t) sen((2k + 1)x) ( 1) u(x, t) = (2k + 1) 2 π k=1

−

−

(7.8)

Usar la función TEMP.M para aproximar la temperatura u(x, t) cualquiera sea el instante t y la posición x. >>u=temp(x,t)

El programa F INICIO.M (ver Web de la asignatura [15]) vale para simular el dato inicial f (x). >>y=f_inicio(x)

Los dos anteriores programas FUNCTION son compatibles con la arquitectura matricial de MATLAB, y sus respectivos códigos se muestran a continuación. function u=temp(x,t) mx=max(size(x)); mt=max(size(t)); k=0:400;n=2*k+1; for jx=1:mx for jt=1:mt u(jx,jt)=sum((-1).^k.*exp(-2*n.^2.*... t(jt)).*sin(n.*x(jx))./n.^2); end end u=4*u/pi; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y=f_inicio(x) y=x.*(0
on del calor siguiente Ejercicio 7.5 Considérese el problema de conducci´ 0 < x < π, t > 0, T t (x, t) 3T xx (x, t) = 0, T (0, t) = T (π, t) = 0, T (x, 0) = sen (x) 7 sen (3x) + sen (5x).

−

−

(7.9)

Si el tiempo t se mide en segundos, estime el momento en que la temperatura T se hace m´ axima en los primeros 10 segundos, en el punto π/6 de la barra.

Clase Pr´ actica 7

114

A continuaci´ on se ofrecen indicaciones para seguir un m´ etodo basado en la potencia gr´ afica del Matlab y en nuestra capacidad de apreciación. El alumno puede seguir cualquier otro método, pero, en tal caso, debe ofrecer una breve descripción del mismo. La solución exacta de (7.9) puede obtenerse aplicando el método de separaci´ on de variables y viene dada por T (x, t) = e −3t sen(x) 7e−27t sen(3x) + e−75t sen(5x).

−

Hacer

>>T=’sin(x).*exp(-3*t)-7*sin(3*x).*exp(-27*t)+sin(5*x).*exp(-75*t)’; >>x=pi/6;t=linspace(0,10,50);Z=eval(T);plot(t,Z)

Apreciar el punto donde se sospecha se alcanza el máximo, y hacer de nuevo >>x=pi/6;t=linspace(a,b,50);Z=eval(T);plot(t,Z)

donde esta vez los valores a y b son aproximaciones por defecto y por exceso, respectivamente, del punto de máximo. Repetir el proceso hasta obtener al menos dos d´ıgitos correctos en la aproximaci´ on, es decir, cuando (b a)/a < (1e 02)/2.

−

7.3.2

−

Ecuaci´ on de ondas sobre un intervalo acotado

on del siguiente problema Ejemplo 7.4 Sabiendo que la soluci´ ytt = a 2 yxx , (a = 20) yx (0, t) = y x (π, t) = 0, t 0, y(x, 0) = 1, yt (x, 0) = cos(x), 0

≥

es y(x, t) = 1 +

(7.10)

≤ x ≤ π.

sen (at)cos(x) , a

hacer su gr´ afico en R3 .

Resoluci´ on del ejemplo 7.4 Usaremos el programa ONDA.M. que podemos descargar desde la Web de la asignatura [15] (ver figura 7.3).

Ejercicio 7.6 El movimiento de una cuerda de longitud finita se rige por la EDP ytt (x, t) = 4yxx (x, t) con las condiciones de frontera yx (0, t) = 0, y(π, t) = 0, (t

≥ 0)

Clase Pr´ actica 7

115

Fig. 7.3: Gráfica producida por ONDA.M y las condiciones iniciales y(x, 0) = 0, yt (x, 0) = 2, (0

≤ x ≤ π).

Si la unidad de tiempo es el segundo, se pide: 1. Obtener la trayectoria del movimiento del punto de la cuerda cuya coordenada es x = π/3, durante los diez primeros segundos y estimar los valores m´ aximo y m´ınimo alcanzados por dicho punto en el per´ıodo. 2. Representar gr´ aficamente la posici´ on de la cuerda pasadas dos horas después del momento inicial y calcular el desplazamiento del extremo x = 0 con respecto al nivel y = 0.

Resoluci´ on del ejercicio 7.6 La solución del ejercicio puede obtenerse aplicando el método de separación de variables y viene dada por ∞ 8( 1)n (2n + 1)x cos sen (2n + 1)t . (7.11) y(x, t) = 2 + 1) 2 π(2n n=0



−



 



La serie (7.11) converge para todo x, t R. Un análisis del resto arroja el siguiente resultado: la suma parcial de 400 sumandos aproxima “uni-

∈

Clase Pr´ actica 7

116

formemente” 2 al valor exacto de la serie con dos d´ıgitos significativos correctos. El programa CUERDA se ha diseñado atendiendo al anterior análisis y permite aproximar hasta con dos d´ıgitos correctos a la serie (7.11). A continuación se muestra el código MatLab del programa CUERDA.M. Entrada: vectores x y t. Salida: matriz u. function u=cuerda(x,t) mx=max(size(x)); mt=max(size(t)); k=0:400; n=2*k+1; for jx=1:mx for jt=1:mt u(jx,jt)=sum((-1).^k.*cos(n*x(jx)*.5).*... sin(n*t(jt))./n.^2); end end u=8*u/pi;

Para obtener los gráficos ejecutar: Primera parte del ejercicio 7.6. >>t=linspace(0,10,100); >>u=cuerda(pi/3,t); >>plot(t,u)

Segunda parte del ejercicio 7.6. >>x=linspace(0,pi,100); >>u=cuerda(x,7200); >>plot(x,u)

Obtener las aproximaciones numéricas solicitadas a partir de los gr´ aficos, mediante los comandos max y min o evaluando adecuadamente a CUERDA(x,t).

7.3.3

Herramientas simb´ olicas para EDP

MATLAB implementa varios comandos para EDP, a través del n´ ucleo MAPLE. >>maple(’pdesolve(ecuacion,func(var1,...,varn))’)

resuelve la ecuación en derivadas parciales identificada por ecuacion, para la función f = f unc(x1 ,...,xn ) 2

La precisión que se obtiene no depende de x ni de t.

Clase Pr´ actica 7

117

Ejemplo 7.5 Intentar resolver las siguientes ecuaciones en derivadas parciales aplicando herramientas simb´ olicas MAPLE-MATLAB. a) f xx + 5f xy = 3 b) 3gx + 7gxy = xy c) h tt

xx =

−h

0

d) u xx + uyy = x + y 2

Resoluci´ on del ejemplo 7.5 En lo que sigue _F1 y _F2 son funciones arbitrarias. a) >>maple(’pdesolve(diff(f(x,y),x,x)+5*diff(f(x,y),x,y)=3, f(x,y))’) ans = f(x,y) = 3/2*x^2+_F1(y)+_F2(y-5*x) --------------------------------------

b) >>maple(’pdesolve(3*diff(g(x,y),x)+7*diff(g(x,y),x,y)=x*y,g(x,y))’) ans = g(x,y) = 1/6*x^2*y-7/18*x^2+_F1(y)+exp(-3/7*y)*_F2(x) --------------------------------------

c) >>maple(’pdesolve(diff(h(t,x),t,t)-diff(h(t,x),x,x)=0,h(t,x))’) ans = h(t,x) = _F1(x+t)+_F2(x-t)

on en derivadas parciales del Ejemplo 7.5, apartado d), aplicando Ejercicio 7.7 Resolver la ecuaci´ herramientas simb´ olicas MAPLE-MATLAB. d) u xx + uyy = x + y 2

Emplear la misma sintaxis MAPLE de los ejemplos resueltos 7.5-a), 7.5-b) y 7.5-c).

Matlab y La Solucion de Ecuaciones Diferenciales

Recommend Documents