ECUACIONES DIFERENCIALES

A. Cachafeiro L´ opez – J. Ill´ an Gonz´ alez

Departamento de Matem´ atica Aplicada I

Universidade de Vigo

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

´INDICE GENERAL

1

Comandos simb´ olicos para resolver ecuaciones. Campos de direcciones. Isoclinas 1.1

1.2 1.3

2

. . . . . . . . . . . . . . . . . . . . . . . . DSOLVE . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

M´ etodo de separaci´ on de variables. Ecuaciones homog´ eneas, exactas, impl´ıcitas y de Bernoulli 2.1 2.2 2.3 2.4 2.5 2.6

3

´ COMANDOS SIMBOLICOS . . . . . . . . . . . . . . . . . 1.1.1 Comandos para la derivaci´on . . . . . . . . . . . . . 1.1.2 Problemas resueltos . . . . . . . . . . . . . . . . . . 1.1.3 El comando DSOLVE. Resoluci´on de PVI . . . . . 1.1.4 Ejemplos de EEDD que no pueden ser resueltas por 1.1.5 Funci´on LAMBERTW(Z) . . . . . . . . . . . . . . CAMPO DE DIRECCIONES DE UNA ED . . . . . . . . TRAYECTORIAS ISOCLINAS . . . . . . . . . . . . . . . 1.3.1 M´etodo de las isoclinas . . . . . . . . . . . . . . . .

1

ECUACIONES EN VARIABLES SEPARADAS . . . ´ ECUACIONES DIFERENCIALES HOMOGENEAS ECUACIONES DIFERENCIALES EXACTAS . . . . ´ DE CLAIRAUT . . . . . . . . . . . . . ECUACION ´ DE BERNOULLI . . . . . . . . . . . . ECUACION FACTOR INTEGRANTE . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

11 . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

Resoluci´ on simb´ olica de sistemas diferenciales 3.1 3.2

´ SISTEMA LINEAL HOMOGENEO DE COEFICIENTES CONSTANTES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SISTEMA LINEAL COMPLETO DE COEFICIENTES CONSTANTES . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

i

1 1 2 2 3 4 5 8 8

11 14 14 18 19 21

25 25 27

´Indice General

ii ´ DE SISTEMAS DIFERENCIALES RESOLUCION ´ LINEALES HOMOGENEOS DE COEFICIENTES CONSTANTES POR LA VIA ESPECTRAL . . . . . . . . . . . . . . 3.3.1 C´alculo de valores y vectores propios via MATLAB . . . . . . 3.3.2 El m´etodo espectral . . . . . . . . . . . . . . . . . . . . . . . . 3.4 EEDD LINEALES DE ORDEN SUPERIOR . . . . . . . . . . . . . . 3.4.1 EEDD lineales de orden superior homog´eneas . . . . . . . . . 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´etodo de coeficientes indeterminados (MCI) . . . . . . . . . 3.7.2 M´etodo de variaci´on de par´ametros (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´etodo basado en el c´alculo aproximado de las ra´ıces caracter´ıstica . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3

4

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

30 30 31 32 32 32

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

35 35 36

. . . . . . . . . . . . . . . . . . . . . de la . . .

38 38 40 43 44 45 45

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ecuaci´on . . . . . .

Transformada de Laplace 4.1

4.2 4.3 4.4

4.5 4.6

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

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 5.2

5.3

5.4 5.5 5.6

6

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

M´ etodos num´ ericos para PVI (II) 6.1

6.2 6.3

6.4

7

iii

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

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´etodos de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.1.2 M´etodo de Runge-Kutta de cuarto orden . . . . . . . . . . . . . . . . . . . . . 106

´Indice General

iv

7.2

7.3

7.1.3 Descripci´on del programa RUNGEKUTTA . . . . 7.1.4 Ejercicios . . . . . . . . . . . . . . . . . . . . . . 7.1.5 Estabilidad de un m´etodo num´erico . . . . . . . . 7.1.6 Descripci´on del programa EJE7.4 . . . . . . . . . PROBLEMAS DE CONTORNO . . . . . . . . . . . . . 7.2.1 Planteamiento matem´atico . . . . . . . . . . . . . 7.2.2 Aplicaci´on: Viga columna . . . . . . . . . . . . . ECUACIONES EN DERIVADAS PARCIALES . . . . . 7.3.1 Ecuaci´on del calor en una barra de longitud finita 7.3.2 Ecuaci´on de ondas sobre un intervalo acotado . . 7.3.3 Herramientas simb´olicas para EDP . . . . . . . .

Bibliograf´ıa de pr´ acticas y sitios WEB

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

106 107 108 109 110 110 111 111 111 114 116 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´afico 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´afico del ejemplo 1.11 . . . . . . . . . . . . . . . . . . Gr´afico del ejemplo 1.12-a . . . . . . . . . . . . . . . . . Gr´afico 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´afico del ejemplo 2.6 . Gr´afico del ejemplo 2.7 . Gr´afico del ejemplo 2.8 . Gr´afico del ejemplo 2.9 . Gr´afico del ejemplo 2.10

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

12 16 17 18 19 20

3.1 3.2

Evidencia gr´afica de que Prueba(2) debe ser nula en [−1, 1] . . . . . . . . . . . . . . . iv) iii) ii) i) Gr´aficas de Z = y1 − 6y1 + 9y1 − 24y1 − 20y1 en [−2, 2] y [20, 60] . . . . . . . . .

30 47

4.1 4.2

Representaci´on gr´afica de γ(t) en [−4, 4] . . . . . . . . . . . . . . . . . . . . . . . . . Representaci´on gr´afica de la entrada g(t) y la salida I(t) . . . . . . . . . . . . . . . .

52 59

5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8

Gr´afica de J1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gr´afica de J0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gr´aficas superpuestas de J1 y J0 . . . . . . . . . . . . . . . . . . . . Gr´aficas de Bessel J3/2 y J1/2 . . . . . . . . . . . . . . . . . . . . . M´etodo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . Interpretaci´on del m´etodo 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

66 66 67 67 71 71 73 73

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

v

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

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 . . . . . . . . . . . . Regla del trapecio . . . . . . . . . . . . . . . . . . . . . . . . . Soluciones exactas y aproximadas del ejemplo 5.5 (n = 21, 41) Ilustraci´on de un m´etodo de dos pasos . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11

Gr´afico producido por S2 EULER al resolver el ejercicio 6.1 Gr´afico producido por S2 EULER al resolver el ejercicio 6.2 Gr´afico producido por S2 EULER al resolver el ejercicio 6.3 Gr´afico obtenido al resolver el ejemplo 6.1 . . . . . . . . . . Gr´aficas de S2 EULER.M e INTERPOLA.M. Ejercicio 6.4 . RK CELL e INTERPOLA.M aplicados al ejercicio 6.5 . . . Soluci´on gr´afica 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´afica del m´etodo de interpolaci´on . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. 90 . 91 . 92 . 93 . 96 . 96 . 97 . 98 . 100 . 100 . 103

7.1 7.2 7.3

Gr´aficas producidas por EJE7.4.M . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Gr´afica producida por CALOR.M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Gr´afica producida por ONDA.M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

. . . . . . . . . . .

74 77 78 81 83

INDICACIONES GENERALES Antes de comenzar la sesi´on 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˜ no, 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´on de la clase, para trabajar simult´aneamente con ambas. 3. (ML) Crear carpetas personales temporales y declarar la correspondiente trayectoria (PATH) mediante el men´ u principal: FILE\SET PATH. Esto es esencial para la ejecuci´on 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´onde se encuentra ´esta, y la sit´ ua 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´on 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´on de la clase.

vii

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

OBJETIVOS • Abordar la soluci´on de ecuaciones diferenciales (EEDD) y problemas de valor inicial desde una perspectiva experimental, mediante la aplicaci´on de los recursos simb´olicos del sistema MATLABr .

• Contribuir con el uso de herramientas gr´aficas al estudio de las EEDD a trav´es de la visualizaci´on 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´on derivada de f respecto a x. >> syms x >> diff(f,x) Para hallar la derivada en´esima hacemos >>diff(f,x,n) El s´ımbolo Df representa simb´olicamente la derivada de f . Para derivadas sucesivas usamos D2f, D3f, etc.

1

Clase Pr´actica 1

2

1.1.2

Problemas resueltos

Ejemplo 1.1 Dada la ecuaci´on my 00 = −mg + ky 0 , comprobar que y(t) = A + Bekt/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 Ejemplo 1.2 Comprobar que la funci´on u(x, t) = e−π

2 c2 t/l2

sen(πx/l)

es soluci´on de la ecuaci´on del calor ut (x, t) = c2 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 Nota 1.1 No ha sido necesario declarar simb´olicas a las variables u, t, x y l.

1.1.3

El comando DSOLVE. Resoluci´ on de PVI

El comando DSOLVE es de tipo simb´olico y permite obtener la soluci´on 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´olicas 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 0 = 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 0 = −y − 2t. Obtener la soluci´on local con la C.I. (−2, 1) y representarla. 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 0 = −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´on, entonces el propio DSOLVE ofrece, en algunos casos, la soluci´on en la forma impl´ıcita.

Clase Pr´actica 1

4

Fig. 1.1: Gr´afico de las soluciones de los ejemplos 1.3 y 1.4 Ejemplo 1.6 Intentar resolver y 0 = 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 00 = −y + cos(2x), con las condiciones iniciales y(0) = 1, y 0 (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´asicas elementales m´as 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´on 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´on. Otras funciones especiales incorporadas al sistema MATLAB-MAPLE son BERNOULLI, BESSElJ, DILOG, etc. Para ver el cat´alogo 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´on que hemos obtenido al ejecutar DSOLVE, es cuando podemos despejarle, es decir, cuando somos capaces de llegar a lambertw(E1 (x)) = E2 (x), donde E1 y E2 son ciertas expresiones en la variable x. Entonces, atendiendo a la formulaci´on impl´ıcita que define a esta funci´on, se tiene que E2 (x) exp(E2 (x)) = E1 (x). Ejercicio 1.1 Resolver las siguientes EEDD y PVI. a) (x2 + 1)dy + (y 2 + 1)dx = 0, y(0) = 1/2 b) y 0 = y cos(x)/(1 + y 2 ), y(0) = 1 c) y 000 − 3y 00 + 3y 0 − y = x d) y 000 − 3y 00 + 3y 0 − y = x, y(0) = 0, y 0 (0) = 1, y 00 (0) = 0

1.2

CAMPO DE DIRECCIONES DE UNA ED

Dada la ED y 0 = f (x, y) , si y = y(x) es soluci´on y (x0 , y0 ) es un punto de la gr´afica de y(x), entonces la pendiente de y = y(x) en (x0 , y0 ) es y 0 (x0 ) que es igual a f (x0 , y0 ).

Clase Pr´actica 1

6

Fig. 1.2: Campo de direcciones y soluci´on particular. Ejemplo 1.8 Lo anterior significa que, para conocer la pendiente de la soluci´on 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´actica es imposible hacer esto con todos los puntos del plano, pero nada nos impide hacerlo con tantos puntos como queramos, configurando un gr´afico que llamaremos campo de direcciones de la ED. ¿C´ omo construir “el campo de direcciones” asociado a una ED? Consideremos la ED gen´erica y 0 (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´on gr´afica acerca de la soluci´on de la ED que a´ un 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´afico del campo de direcciones. Los vectores a considerar para obtener el campo de direcciones ser´an (1, y 0 ) = (1, f (x, y)). Ejemplo 1.8 Construir el campo de direcciones de la ED y 0 = x, y representar la soluci´on particular 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´on 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´on particular se obtiene artesanalmente como y = x2 /2, mediante separaci´on de variables, o bien, aplicando DSOLVE. A continuaci´on aparecen las o´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 0 = −x/y y representar la soluci´ on particular que cumple la CI y(1) = 1. 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´on inicial. La expresi´on dentro de INLINE que define a la funci´on debe escribirse entre ap´ostrofes, es decir, 0 −x./y 0 . El punto debe preceder a /, ∗ y ^ pues se opera con vectores componente a componente. La soluci´on gr´afica est´a 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 0 = y − x,

ii) y 0 = 2x.

Fig. 1.4: Representaci´on (a)

1.3

Fig. 1.5: Representaci´on (b)

TRAYECTORIAS ISOCLINAS

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

1.3.1

M´ etodo de las isoclinas

El m´etodo de las isoclinas es una variante de las ideas antes descritas. Los puntos del plano por los que pasa una soluci´on con pendiente k, son los puntos de la curva de ecuaci´on 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´es identificar la isoclina para la pendiente 0 pues las soluciones tendr´an generalmente un m´aximo o un m´ınimo al pasar por esta isoclina. Ejemplo 1.10 1 Representar las isoclinas de la ED y 0 = x + y 2 . 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 0 = x2 + y 2 . 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 0 = 2x − y. (b) ¿Qu´e tipo de curvas son dichas isoclinas? (c) Representar las isoclinas correspondientes a k = 0, y k = 2. (d) ¿Qu´e particularidad tiene la correspondiente a k = 2?

Clase Pr´actica 1

10

Fig. 1.8: Gr´afico 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 0 = sen(x) + y. Desarrollo, paso a paso, de la soluci´ on del ejemplo 1.13 (figura 1.10)

Fig. 1.10: Gr´afico 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´on de ED. 2. Aplicar m´etodos de cuadraturas exactas estudiados en clase.

2.1

ECUACIONES EN VARIABLES SEPARADAS

Son de la forma y 0 (x) = f (x)h(y). Para resolver el PCI y 0 = f (x)h(y), y(x0 ) = y0 , 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´on general de la ED),

y calculamos el par´ametro C usando la CI. As´ı C = G(y0 ) − F (x0 ) y por lo tanto G(y) = F (x) + G(y0 ) − F (x0 ),

(Soluci´on del PCI).

Ejemplo 2.1 Resolver el PVI, y 0 = −2xy, y(0) = 2, y representar gr´aficamente la soluci´ on.

11

Clase Pr´actica 2

12

Fig. 2.1: Soluci´on del ejemplo 2.1 Resoluci´ on del ejemplo 2.1 Resolveremos el PVI con MATLAB siguiendo el m´etodo 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´afico de la soluci´on 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´on general sol =exp(-x^2+2*c) y la soluci´on 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 0 = sen(x)(y − 1), y(0) = 2. Resoluci´ on del ejemplo 2.2 >>f=’sin(x)’

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 − (1 + y 2 )dy = 0, y(0) = 1. Resoluci´ on del ejemplo 2.3 Despejando obtenemos y 0 = y cos(x)/(1 + y 2 ) 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´on de la anterior respuesta depender´a del consumidor. Si tenemos en cuenta que Lambertw(x)=w es la funci´on definida impl´ıcitamente por wew = x, lo que se obtiene es y 2 ∗ exp(y 2 ) = K ∗ exp(2 ∗ sen(x)). Ejemplo 2.4 Resolver la ED anterior, es decir, y cos(x)dx − (1 + y 2 )dy = 0, y(0) = 1, separando las variables. Nota sobre el ejemplo 2.3. Supongamos que G(Lambertw(x), y, x, C) = 0, es la soluci´on general (con par´ametro C), de una ED y 0 = 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´on w(x)ew(x) = x, entonces la nueva expresi´on de la soluci´on 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) + y 2 /2 + c, que es la soluci´on general de la ED. Para hallar la soluci´on 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´on particular: sen(x) = ln(y) + (y 2 − 1)/2. Ejercicio 2.1 Aplicar herramientas MATLAB para resolver los siguientes ejercicios extra´ıdos del Bolet´ın 2. a) y 0 + y 2 sen(x) b) 4xy dx + (x2 + 1)dy = 0 p =0 c) xy 0 = y + x2 − y 2 d) (x − y)dy = (x + 3y)dx 2 2 e) (x + 1)dy = (y + 1)dx, y(0) = 1/2 f ) (1 + ex )yy 0 = ey , y(0) = 0

2.2

´ ECUACIONES DIFERENCIALES HOMOGENEAS

Ejemplo 2.5 Resolver la ED (x2 − y 2 )dx + xydy = 0. Resoluci´ on del ejemplo 2.5 Es claro que es una ED homog´enea pues P (x, y) = x2 − y 2 y Q(x, y) = xy son funciones homog´eneas 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´olicos del MATLAB se escribe la ED en la forma y 0 + (x2 − y 2 )/(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´on (y/x)2 = −2 log(x) + C1.

2.3

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´on general es F (x, y) = C. Para la obtenci´on de F se procede como sigue. Integramos a F (x, y) respecto a una sola de sus variables, digamos x. Z 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´olo resta hallar a f (y) Si derivamos con respecto a y R ∂( P (x, y)dx) + f 0 (y). Q= ∂y Luego R ∂( P (x, y)dx) 0 f (y) = Q − . ∂y Integrando respecto a y la anterior expresi´on obtenemos f (y). La soluci´on general es F (x, y) = C, C ∈ R. Ejemplo 2.6 Comprobar que la ED (2x + y)dx + (x − 3y)dy = 0, es exacta. Resolver dicha ecuaci´ on 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´afico del ejemplo 2.6 Por lo tanto la soluci´on general es x2 + xy − 3/2y 2 = C. Para representar gr´aficamente 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 0 =

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

ciones. 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´on general de la ED es 1/3x3 − 1/3y 3 = C. Para representar las soluciones escribimos (ver figura 2.3)

Clase Pr´actica 2

17

Fig. 2.3: Gr´afico 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 = 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 F =sin(y)*sinh(x) Por lo tanto, la soluci´on del ejemplo 2.8 es sen(y) ∗ senh(x) = C. Para 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)

Clase Pr´actica 2

18

Fig. 2.4: Gr´afico del ejemplo 2.8

2.4

´ DE CLAIRAUT ECUACION

Recordemos que la ecuaci´on de Clairaut es de la forma y = xy 0 + f (y 0 ).

(2.1)

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

dp = 0 entonces y 0 = C y por tanto teniendo en cuenta (2.1) se obtiene que dx y = Cx + f (C), C ∈ R,

(2.2)

siendo (2.2) un haz de rectas, todas ellas soluci´on de la ecuaci´on original (2.1). Si x + f 0 (p) = 0, usando (2.1) se obtiene la soluci´on singular en forma param´etrica: x = −f 0 (p), y = −f 0 (p)p + f (p). En general no es necesario eliminar p para obtener una ecuaci´on de la forma G(x, y) = 0, y de hecho, podr´ıa resultar muy dif´ıcil o imposible. Ejemplo 2.9 Resolver la ecuaci´on de Clairaut y = xy 0 − 41 (y 0 )2 . Representar la soluci´on general obtenida, as´ı como la soluci´on singular.

Clase Pr´actica 2

19

Fig. 2.5: Gr´afico del ejemplo 2.9 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´on es la envolvente de la familia integral de la ecuaci´on. Para dar una idea gr´afica aproximada del conjunto soluci´on 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’) Ejercicio 2.2 Resolver utilizando herramientas MATLAB, la ecuaci´on de Clairaut y = xy 0 + ln(y 0 ).

2.5

´ DE BERNOULLI ECUACION

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

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´afica en el intervalo [1, 3] se obtiene ejecutando (ver figura 2.6) >>ezplot(sol,[1,3]) Ejercicio 2.3 1. Resolver la ED 3xy 0 − 2y = 2. Resolver el PVI 3xy 0 − 2y =

x3 . y2 x3 , y2

y(−1) = 1.

Fig. 2.6: Gr´afico del ejemplo 2.10 La ED de Bernoulli y 0 + p(x)y = q(x)y n se resuelve transform´andola en lineal, para lo cual se procede como sigue: 1 Se divide la ED por y n y se hace el cambio z = yn−1 , obteniendo una ED lineal en z(x). >>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))) Ejemplo 2.11 Resolver la ecuaci´on y 0 + y/x = log(x)y 3 , aplicando el m´etodo antes explicado.

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´on 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/y 2 para obtener la soluci´on.

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´on µ = µ(x, y) tal que µP (x, y)dx + µQ(x, y)dy = 0 s´ı es exacta. Los casos aparentemente m´as sencillos son aquellos en que puede encontrarse un factor µ, llamado “factor integrante” (f.i.), que s´olo depende de x o s´olo depende de y. Un ejemplo de f.i. que depende de ambas variables x e y es aqu´el que encontramos para la ED (xy 2 − yx2 )dx + x3 dy = 0. (2.3) Si llamamos P (x, y) = xy 2 − yx2 y Q(x, y) = x3 comprobaremos que 1 1 = 2 2 xP + yQ xy es un f.i. para la ED (2.3). La comprobaci´on de lo anteriormente dicho podemos hacerla f´acilmente 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´on es log(x) − x/y = C. Nota 2.1 El anterior procedimiento puede aplicarse a cualquier ecuaci´on homog´enea. El alumno 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. y0 =

√

x + y.

(2.4)

Resoluci´ on del ejemplo 2.12 0 La soluci´on puede obtenerse mediante √ el cambio z = x + y, ya que es de la forma: y = G(ax + by + c) un no tenemos explicaci´on para el siguiente con a = b = 1, c = 0 y G(t) = t. Sin embargo, a´ resultado MATLAB: >>dsolve(’Dy=sqrt(x+y)’,’x’) ans =0 Es simple comprobar que y = 0 no es soluci´on de la ecuaci´on. 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

√ y 0 = ± x + y. Mediante comprobaci´on directa vemos que y = 1 − x es solamente soluci´on de √ y0 = − x + 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´on: (1 ±

√

x + y) exp(−1 ±

√

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

Para nuestro problema original, tenemos como respuesta la siguiente ecuaci´on √ √ (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´on y = (y 0 )2 − x, que es de la forma, ”resuelta respecto a y”, del tipo ”Lagrange” y = (−1)x + (y 0 )2 . Hacemos y 0 = p y derivamos con respecto a x para obtener la ecuaci´on lineal en x = x(p) (p + 1)x0 = 2p, que se integra f´acilmente obteniendo x = 2p − log(p + 1) + C, y = p2 − 2p + log(p + 1) − C,

(2.6)

teniendo en cuenta que y = p2 − x. √ Aunque parezca a primera vista que es dif´ıcil, en realidad es muy f´acil 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 0 = G(ax + by + c), y como Lagrange y = xf (y 0 ) + g(y 0 ). Tambi´en podemos ver la ED del ejemplo (2.12) como la ecuaci´on no exacta √ x + y dx − dy = 0. Lo m´as interesante ahora es que resulta infructuoso el intento de hallar un f.i. que s´olo dependa de x o que s´olo 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,

√ 1 = −µx , µy x + y + µ √ 2 x+y

Clase Pr´actica 2

24 o equivalentemente

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

Debido a la forma de µ , y a la expresi´on 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) = ϕ0 (x + y), (x, y) = ϕ0 (x + y), ∂x ∂y siendo ϕ0 la derivada ordinaria de la funci´on ϕ(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 µP dx + µQdy = 0 es exacta procedemos al c´alculo del potencial F (x, y). La integral general es F (x, y) = C. r y 0 Ejercicio 2.4 Resolver el siguiente PVI (Bolet´ın 2) y = , y(1) = 4. x

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

OBJETIVOS Ejercitar el uso combinado de t´ecnicas matem´aticas y herramientas simb´olicas 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 x0 = −5x + 3y y 0 = −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´on y escribimos >>S=dsolve(’Dx=-5*x+3*y,Dy=-2*x-10*y’) MATLAB devuelve las componentes simb´olicas de la soluci´on: S= x: [1x1 sym] y: [1x1 sym]

25

Clase Pr´actica 3

26

Para obtener expl´ıcitamente la soluci´on 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 = C1 (−2e−8t + 3e−7t ) + 3C2 (e−7t − e−8t ) y = 2C1 (−e−7t + e−8t ) + C2 (3e−8t − 2e−7t ) Alternativamente se puede obtener directamente la respuesta indicando las componentes del vector soluci´on. >>[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´eneo, de coeficientes constantes, es seguir la via matricial directa. Para ello comenzamos escribiendo el sistema en forma matricial:      0  −5 3 x x x = =A −2 −10 y y y Conocer la soluci´on 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) = etA . Para calcular etA 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: x0 − y = e t y 0 + 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´etodo de Variaci´ on de Par´ ametros

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´on mediante el MVP usaremos la f´ormula correspondiente para obtener una soluci´on particular del sistema diferencial lineal X 0 = A(t)X + B(t). La soluci´on general del mismo es Z t

R(s)−1 B(s)ds,

X(t) = R(t)V + R(t)

(3.1)

0

siendo R(t) una matriz fundamental y V ∈ Rn . La f´ormula (3.1) corresponde al caso general A(t) de coeficientes variables. El caso particular de coeficientes constantes aparece a continuaci´on. Recordemos que, en caso de sistemas diferenciales homog´eneos con coeficientes constantes, X 0 = AX, una matriz fundamental est´a dada por R(t) = etA . Luego, si se trata de un sistema completo de coeficientes constantes X 0 = AX + B(t), una soluci´on particular Xp est´a dada por Xp = e

tA

Z

t

e−sA B(s)ds.

0

Para resolver el sistema diferencial lineal del ejemplo 3.2, es decir x0 − y = et y 0 + 5x + 2y = sen(3 + t), con x(0) = 1, y(0) = −1, usaremos la f´ormula general de variaci´on de par´ametros para hallar una soluci´on particular Xp del sistema no homog´eneo. >>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´on expl´ıcita de la soluci´on particular Xp se compone de una cadena de caracteres demasiado larga y por tanto de reducida utilidad, en efecto, si hacemos la conversi´on de SYM a STRING mediante el comando CHAR, tendremos que >>size(char(Xp)) ans = 1

1669

es decir, ¡Xp est´a compuesto de 1669 caracteres! La figura 3.1 muestra gr´aficamente el resultado de ejecutar ezplot(Prueba(2),[-1,1]). El rango de variaci´on 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 x0 + x + 2y = cos(t) + sen(t) + e−t y 0 − 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´afica 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´eneo, 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 {eλ1 t u1 , · · · , eλn t un }, es un sistema o conjunto fundamental de soluciones (SFS) para el sistema diferencial X 0 = AX. Ejemplo 3.4 Resolver por la via espectral, y aplicando herramientas MATLAB, el siguiente problema con valores iniciales     1 2 −1 −1 1  X(t) X(0) =  0  X0 =  1 0 4 −4 5 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´on general del sistema. Luego, la u ´nica soluci´on estar´a dada por la soluci´on del siguiente sistema     c1 −1 X(0) = R(0)  c2  =  0  c3 0 es decir, 

   c1 −1  c2  = R(0)−1  0  c3 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´on diversos ejemplos de EEDD homog´eneas y completas de coeficientes constantes para las que se obtiene la soluci´on general as´ı como ciertas soluciones particulares. Se presenta tambi´en 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 00 + 4y 0 + 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´ervese que −2 es una ra´ız doble del polinomio caracter´ıstico Ejemplo 3.6 Resolver el PVI y 00 + 4y 0 + 4y = 0, y(0) = 1, y 0 (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 C1 = 1 y C2 = 6. En todos los casos las respuestas son exactas, con f´ormulas cerradas.

3.4.2

EEDD lineales de orden superior completas

Ejemplo 3.7 Resolver la ED y 00 + y 0 + 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) √ Notar que el polinomio caracter´ıstico tiene ra´ıces complejas conjugadas (−1 ± i 15)/2. Ejemplo 3.8 Resolver y 000 = log(x)/x2 , y(1) = 0, y 0 (1) = 1, y 00 (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´on exacta en t´erminos del cat´alogo de funciones elementales MATLAB, entonces DSOLVE no puede ofrecer una soluci´on. Existen otras causas que impiden a DSOLVE llegar a una respuesta exacta pero su consideraci´on 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´etodos alternativos que produzcan soluciones, quiz´as aproximadas pero mucho m´as potables. Para ilustrar lo dicho, a continuaci´on se ofrece un ejemplo de EDL homog´enea de cuarto orden, cuya soluci´on se obtiene de forma aproximada, siguiendo un m´etodo basado en nuestros conocimientos te´oricos. En el ap´endice se ampl´ıan los detalles t´ecnicos del m´etodo y se ejerce una cr´ıtica necesaria. Ejemplo 3.9 Considere la EDLH (3 L[y] = 0, donde L[y] = y iv) − 6y iii) + 9y ii) − 24y i) − 20y. Intente resolverla, es decir, halle informaci´on u ´til 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´endice), calculando un SFS a partir del polinomio caracter´ıstico P = m4 − 6m3 + 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´erminos de radicales, ecuaciones polin´omicas hasta el grado 4. ROOTS calcula aproximadamente las ra´ıces de cualquier polinomio mediante un 3 4

Ecuaci´ on Diferencial Lineal y Homog´enea La respuesta puede cambiar seg´ un la versi´on MATLAB o el software que se utilice

Clase Pr´actica 3

34

procedimiento muy eficiente y preciso. No obstante, los resultados num´ericos obtenidos mediante SOLVE deber´ıan ser superiores, en algunos casos, a los producidos por ROOTS. El nivel de coincidencia entre ambos es del orden 10−15 en el caso que nos ocupa, seg´ un se aprecia en las pruebas num´ericas que se muestran a continuaci´on. >>format 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=roots([1 -6 9 -24 -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.329070518200751e-015 El correspondiente SFS es y1 y2 y3 y4

= e0.6653067648628088 t sen(2.373811784418384 t), = e 0.6653067648628088 t cos(2.373811784418384 t), = e5.291305619991860 t = e−0.6219191497174775 t

´ Ver en el APENDICE c´omo intentamos valorar la calidad de esta soluci´on aproximada. Ejercicio 3.1 Dadas las siguientes EEDD a.- y iv) − 10y iii) + 3y ii) − 2y i) − 20y = 0, b.- y iii) − 2y ii) + 3y i) = 0, c.- y iv) − 6y iii) + 9y ii) − 24y i) + 50y = 1, obtener una soluci´on particular para la ED del apartado c) mediante el m´etodo de coeficientes indeterminados (MCI) sin ayuda del MATLAB, y un SFS para a), b) y c).

Clase Pr´actica 3

3.5

35

´ DE EULER Y DE EULER ECUACION GENERALIZADA

La ecuaci´on de Euler de segundo orden es de la forma αx2 y 00 (x) + βxy 0 (x) + γy(x) = δ(x) Recordemos que ´esta se transforma en una ED de coeficientes constantes mediante el cambio x = et , y que la ecuaci´on de Euler generalizada de segundo orden es de la forma α(ax + b)2 y 00 (x) + β(ax + b)y 0 (x) + γy(x) = δ(x) Ejemplo 3.10 Resolver (2x + 1)2 y 00 − (4x + 2)y 0 − 12y = 8x. Resoluci´ on del ejemplo 3.10 mediante MATLAB 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*x+1)+ C1*(2*x+1)^3+C2/(2*x+1) Ejercicio 3.2 ¿A cu´ales de las siguientes EEDD, despu´es de ser transformadas en EDL de coeficientes constantes, puede aplicarse el MCI? ¿Cu´al es la alternativa al MCI? Intentar resolverlas mediante DSOLVE y mediante m´etodos artesanales. Contrastar los resultados. 1. (2x + 1)2 y 00 − (4x + 2)y 0 − 12y = 8x 2. x2 y 00 + 5xy 0 − 2y = −5 log(x) 3. (10x + 50)2 y 00 − 3(x + 5)y 0 + 7y = x−2 4. x3 y 00 − 6x2 y 0 + 2xy = x + 1 5. 5x2 y 00 − 6xy 0 + y = ex

3.6

´ CALCULO DE WRONSKIANOS Y APLICACIONES

En lo que sigue se utilizar´a el wronskiano para el estudio de la independencia lineal de sistemas de funciones. Tambi´en se aplicar´a en la obtenci´on de una EDLH conocidas algunas de sus soluciones. El wronskiano juega un papel relevante en la teor´ıa y pr´actica de las EDL. Si y1 e y2 son funciones definidas sobre el intervalo I, entonces el wronskiano W (y1 , y2 ) de ´estas se define como:   y1 (t) y2 (t) W (y1 , y2 )(t) = det , t ∈ I. y10 (t) y20 (t)

Clase Pr´actica 3

36

3.6.1

Wronskiano e independencia lineal

Teorema 3.2 Sean x1 (t) y x2 (t) dos soluciones de una ED del tipo x00 (t) + p(t)x0 (t) + q(t)x(t) = 0, con p(t) y q(t) continuas en I. Entonces se cumple que: x1 y x2 son linealmente independientes en I ⇔ W (x1 , x2 )(t) 6= 0 ∀t ∈ I ⇔ ∃t ∈ I tal que W (x1 , x2 )(t) 6= 0. Nota 3.1 Sean x1 y x2 de clase 1 en I, no necesariamente soluciones de una EDLH. Entonces 1) W (x1 , x2 )(t) 6= 0 para alg´ un t ∈ I ⇒ x1 y x2 son linealmente independientes en I. 2) Sin embargo, es posible que x1 y x2 sean linealmente independientes en I y W (x1 , x2 )(t) = 0 en alg´ un I. Por ejemplo, t y t2 son l.i. en [−1/2, 1/2] y W (t, t2 )(0) = 0. Ejemplo 3.11 Calcular el Wronskiano de e−2t y senh(t) − 2 cosh(2t). Resoluci´ on del ejemplo 3.11 >>a=’[exp(-2*t), sinh(t)-2*cosh(2*t)] ’; >> b=diff(a,t) % tambien diff(a) b=[-2*exp(-2*t),cosh(t)-4*sinh(2*t)] >> 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(c) %(es el Wronskiano 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 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) = et 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) = −t + 1, h(t) = 3t2 + 2t − 4 en (−1, 2). 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

Clase Pr´actica 3

38 a) y1 = ex , y2 = e2x , y3 = e3x b) y1 = 1, y2 = x, y3 = 1/(1 + x2 ) c) y1 = sen(x), y2 = x d) y1 = ex , y2 = cos(3x) e) y1 = x, y2 = x2 + 1, y3 = x2 + 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´on es W = 0, siendo dy = y 0 , d2y = y 00 , etc. Si se selecciona un intervalo dominio I0 , entonces I0 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. Ejemplo 3.16 Obtener una soluci´on particular de la ED x00 + 3x0 − 4x = 3 sen(t) + 7e−4t − e2t

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´on general de la homog´enea. Teniendo en cuenta la forma de la parte completa de la ED, se busca una soluci´on particular xp que sea combinaci´on lineal de x1 = sen(t), x2 = cos(t), x3 = te−4t , y x4 = e2t . >>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´on particular buscada. Para determinar los coeficientes A, B, C y D se debe sustituir xp 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´ ua eqn y sus derivadas 1ra, 2da y 3ra en un punto (en este caso tomamos el punto 0), obteni´endose un sistema cuya resoluci´on 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´on se resuelve el sistema lineal formado por las 4 ecuaciones eqn1, eqn2, eqn3 y eqn4, con las inc´ognitas A, B, C, y D. >>[A B C D]=solve(eqn1,eqn2,eqn3,eqn4) A =-15/34 B =-9/34 C =-7/5 D =-1/6 Una soluci´on 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´alculo de los coeficientes no es aconsejable sin la ayuda de instrumentos de c´alculo, pero s´ı parece ser m´as 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)

Ejemplo 3.17 Obtener una soluci´on particular de la ED t2 x00 + tx0 − x = sen(t). Se utilizar´a que x1 = t y x2 = 1/t son soluciones l.i. de la ED homog´enea. (Se pueden obtener con DSOLVE) Resoluci´ on del ejemplo 3.17 La soluci´on general de la homog´enea es xh = At + B/t, A, B ∈ R y se busca una particular de la completa de la forma xp = A(t)t + B(t)/t. Sabemos que A0 (t) y B 0 (t) son soluciones del sistema   0    t 1/t A (t) 0 = 1 −1/t2 B 0 (t) sen(t)/t2 Procedemos a resolver dicho sistema.

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 Ejemplo 3.18 Obtener una soluci´on particular de la ED x00 − 6x0 + 9x = 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 A0 (t) y B 0 (t) son soluciones del sistema  0     3t e te3t A (t) 0 = (e3t )0 (te3t )0 B 0 (t) e3t /t2 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

e3t . t2

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´on particular obtenida. Ejemplo 3.19 Resolver mediante el MVP la ED x00 + 3x0 + 2x = e3t 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´on particular xp cuya expresi´on pretendemos que sea xp = 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´on, pero en cualquier caso es mucho m´as c´omodo usar directamente DSOLVE que aplica autom´aticamente el principio de superposici´on.

Clase Pr´actica 3

3.8

43

´ DE ORDEN REDUCCION

Ejemplo 3.20 Resolver la ED x00 + (2/t)x0 + x = 0, usando que f (t) = sen(t)/t es una soluci´ on 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 0 = y se llega a una ED de primer orden en y. y 0 + (2/t + 2f 0 (t)/f (t))y = 0. Resolviendo esta ED se obtiene su soluci´on general. y = F (t, C1 ). Como y = v 0 , integrando se obtiene v = G(t, C1 , C2 ). Finalmente haciendo x=f ∗v se obtiene la soluci´on general de la ED. Los pasos con MATLAB se indican a continuaci´on. >>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 0 + (2/t + 2f 0 /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 0 , 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 y1 (t) = et , y2 (t) = tet , e y3 (t) = t2 et , son linealmente independientes en [a, b]. 2. Resolver el siguiente problema de CI 2y 00 + 5y 0 + 5y = 0, y(0) = 0, y 0 (0) = 1/2. 3. Resolver 9y iv) − 6y iii) + 46y ii) − 6y i) + 37y = 0. 4. Hallar la soluci´on general del sistema de ecuaciones y 0 + z 0 = exp(x) y + z 00 = x 5. Determinar la soluci´on del sistema x0 = y y0 = z z 0 = −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. (a) y 00 − 4ty 0 + y = 1 (b) y 00 − 4t2 y 0 + y = 0 (c) y 00 − 4ty 0 + ty = 0

Clase Pr´actica 3

45

(d) ty 00 − 4y 0 + y = 0 (e) x00 + 5x0 − 50x = 1/t ´ Ver en el APENDICE una notas relativas a las funciones especiales que pueden aparecer en la soluci´on 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´on Lambertw, mencionaremos otras funciones b´asicas tambi´en incorporadas a MATLAB. Estas son: WhittakerW, WhittakerM, besselj, bessely, Ei, dilog, la funci´on de error erf, Heaviside, Dirac, etc, que tambi´en pueden ser conocidas mediante MHELP. Se trata de funciones especiales para el matem´atico y el f´ısico debido a que son la soluci´on de problemas espec´ıficos, frecuentes y de gran importancia en determinados ´ambitos de la matem´atica, 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´erminos de ciertas funciones especiales que hemos adoptado, y que son bien conocidas en el sentido de que tenemos suficiente informaci´on acerca de ellas: intervalos de convexidad y de crecimiento, m´aximos y m´ınimos locales, comportamiento asint´otico en el infinito, derivabilidad, singularidades, desarrollo en serie de potencias, etc. Las EEDD de coeficientes variables del u ´ltimo listado de ejercicios ser´an estudiadas en el 2do cuatrimestre utilizando un m´etodo relativamente simple, basado en asumir a priori que la soluci´on 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´ un mejor cu´al es la tarea que desempe˜ nan las funciones especiales y que ´estas 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´on lineal homog´enea de coeficientes constantes ay iv) + by iii) + cy ii) + dy i) + 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´on >>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´on 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´agina 33, para ilustrar el nivel de ajuste que se obtiene al calcular una soluci´on aproximada y˜ por este m´etodo. 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[˜ y ]|, en lugar del error en la soluci´ on: k˜ y − yk. 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´aficos 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´oricos. Sin embargo, en el ´ambito experimental los baremos cambian. La calidad de los resultados obtenidos se mide, en este caso, a partir de la estimaci´on 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´a 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 una precisi´on finita.(7 Parece correcto aceptar que los errores experimentales son nulos cuando ´estos son < eps ≈ 2.22e − 016, y ´este no es el caso que nos ocupa.(8 ¿Pero acaso 10−14 no es suficientemente peque˜ no? ¿Podemos aceptar como SFS al obtenido por el m´etodo descrito anteriormente? 6

Agregando title(”) evitamos que en el t´ıtulo de la gr´afica aparezcan cadenas interminables de caracteres. Se llama precisi´ on al n´ umero de d´ıgitos que forman la mantisa. 8 Por eps se representa al epsilon de m´ aquina, es decir, el n´ umero positivo m´as 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)

iii)

ii)

i)

Fig. 3.2: Gr´aficas de Z = y1 − 6y1 + 9y1 − 24y1 − 20y1 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´on y que ´esta, a su vez, se corresponda con la tolerancia de nuestros instrumentos. En este caso tendremos presente que no tenemos informaci´on acerca del error en la soluci´on k˜ y − yk, y que nuestra conclusi´on depender´a del intervalo considerado. La respuesta a la segunda pregunta es negativa, y es que en realidad hemos avanzado muy poco pues la soluci´on general de la ED es de la forma(9 y˜G = C1 eα3 t cos(β t) + C2 eα3 t sen(β t) + C3 eα1 t + C4 eα2 t ,

(3.2)

de modo que las constantes C1 , ..., C4 pueden modificar sustancialmente la escala y por ello podr´ıa ser a´ un m´as 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´on es una estrategia universalmente aceptada. El llamado “m´etodo de los residuos (ponderados)”(10 consiste en aceptar como soluci´on aquella que anula al residuo en un sentido m´as d´ebil que el aqu´ı utilizado. 9 10

Un serio inconveniente es que la funci´ on y˜G depende no linealmente de los par´ametros α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´on, reiteremos que es posible para EEDD lineales y homog´eneas de orden cuatro o menor, calcular exactamente las ra´ıces del polinomio caracter´ıstico utilizando SOLVE sin aplicar NUMERIC, solo que la expresi´on 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´olicos no parecen ser concluyentes, al menos con este enfoque. Para resolver ´este y muchos otros problemas nos apuntamos al uso de los m´etodos num´ericos tradicionales cuyo estudio ser´a abordado en las clases pr´acticas 5, 6 y 7.

´ CLASE PRACTICA 4 Transformada de Laplace

OBJETIVOS Ejercitar el uso combinado de t´ecnicas matem´aticas y herramientas simb´olicas del MATLAB para la resoluci´on, mediante la transformada de Laplace, de problemas de valores iniciales asociados a sistemas diferenciales lineales y ecuaciones lineales de orden superior.

4.1 4.1.1

´ Y CALCULO ´ DEFINICION Definici´ on de la transformada de Laplace

La transformada de Laplace de una funci´on f (t) definida en [0, +∞) es una nueva funci´on que denotaremos por F (s) o f˜(s) o L(f )(s), definida por la integral f˜(s) =

Z

+∞

e−st f (t)dt.

(4.1)

0

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

+∞

Z g(t)dt = lim

0

b→+∞

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´olo 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´on abordaremos el c´alculo 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´alculo automatizado de (4.1) MATLAB ofrece el comando LAPLACE de tipo simb´olico, cuya sintaxis es >>F=laplace(f) donde f es una funci´on escalar de la variable simb´olica t, previamente declarada y F es una funci´on cuya variable es por defecto s. Tambi´en podemos formular >> syms u v >>F=laplace(u^2,v) F = 2/v^3 lo que permite elegir las variables a utilizar. Ejemplo 4.1 Calcular la transformada de Laplace de la funci´on f (t) = 1. 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´on 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´ un par´ametro como simb´olico. 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) = t2 y f (t) = t3 . Nota para el ejercicio 4.1 Recordar que para n ∈ N la f´ormula general viene dada por la expresi´on: n! ten (s) = n+1 , s y para r > −1, r ∈ R, se tiene, a´ un m´as general, que Γ(r + 1) ter (s) = , sr+1 donde Γ(r) se define a continuaci´on. La funci´on Gamma se define como Z +∞ Γ(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´ on 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´on: at f (t)(s) = f˜(s − a), e^

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

Clase Pr´actica 4

52

Fig. 4.1: Representaci´on gr´afica de γ(t) en [−4, 4] Ejercicio 4.3 Calcular f˜ en los siguientes casos y reconocer las propiedades que se han aplicado en cada caso. (a) f (t) = t1/3 (b) f (t) = e−2t t5/6 (c) f (t) = g 0 (t) (d) f (t) = sen(t + a) (e) f (t) = sen(t − 1)u(t − 1) (f ) f (t) = t2 senh(t) En (c) hacer >>laplace(diff(sym(’g(t)’))). Notar que esta sintaxis no exige declaraci´on previa SYMS. Aplicar la misma t´ecnica al apartado (e), es decir, laplace(sym(’Heaviside(t-1)*sin(t-1)’)), donde Heaviside(t-a) es el nombre codificado de la funci´on 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´on F(s), produciendo otra expresi´on f(t). La versi´on 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) =

s3

1 1 1 , F (s) = 4 , F (s) = 5 . +1 s +1 s +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´en 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´ un 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) Nota 4.1 exp(1/5*i*pi*(2*k-1)),k=1..5 son las 5 ra´ıces complejas del n´ umero −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. x0 (t) + 3x(t) = 0, x(0) = 1. Resoluci´ on del ejemplo 4.6. Se ver´a que el comando LAPLACE aplica autom´aticamente la propiedad xe0 = se x − 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 x00 (t) + 3x0 (t) + 2x(t) = sen(t), x(0) = x0 (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 x00 (t) − 3x0 (t) + 2x(t) = 2e−2t , x(0) = 1/6, x0 (0) = 5/3. Ejemplo 4.8 Resolver el siguiente sistema x0 = −x − y y 0 = −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: x00 (t) = x0 (t) − 2y 0 (t) + 2y(t), y 0 (t) = 3y(t) + x(t), y(0) = 0, x(0) = −2, x0 (0) = 1.

4.4 4.4.1

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

La funci´on de salto unitario u(t) =

  0

t>syms t s >>f=’Heaviside(t)’

Clase Pr´actica 4

57

>>laplace(f,t,s) ans=1/s Ejemplo 4.9 Calcular la transformada de la funci´on u(t − a). Teniendo en cuenta la notaci´ on τ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, u(t − a) es la trasladada paralelamente al eje OX de la funci´on f = 1 al punto a. 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´acil de obtener. Es oportuno recordar una importante propiedad vinculada a la funci´on escal´on unitario. Propiedad de la traslaci´ on en t. Sea F (s) = f˜(s) y a > 0. Entonces 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: ^ e−as F (s) = f (t − a)u(t − a). Ejemplo 4.10 Resolver el siguiente problema de condici´on inicial x0 (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´etodos 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

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

  1 −1 g(t) =  0

0>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´on gr´afica de la entrada g(t) y la salida I(t) La representaci´on gr´afica 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:  0 t 6= a δ(t − a) = ∞ t=a 1

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

(4.4)

Clase Pr´actica 4

60 y Z

∞

f (t)δ(t − a)dt = f (a),

(4.5)

−∞

para cualquier funci´on f (t) que sea continua en un intervalo abierto que contiene a t = a. Observemos que δ(t − a) no es una funci´on en el sentido usual ya que toma el valor ∞ en un punto. De la propiedad (4.5) se obtiene que ^ δ(t − a)(s) = e−as La Delta de Dirac se relaciona con la funci´on escal´on unitario de la siguiente forma  Z t 0 t>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. tx00 (t) + (t − 1)x0 (t) − x(t) = 0, x(0) = 5, lim x(t) = 0. t→+∞

Nota 4.3 Esta ED es de 2do orden...¡pero no se establece condici´on alguna sobre x0 (0)! En su lugar se impone una condici´on en el infinito... Notar que t = 0 es una singularidad (¿por qu´e?). Quiz´ as esto u ´ltimo est´e asociado a esas “extra˜ nas condiciones” que se imponen a la soluci´on. Una idea: si t ∈ [−ε, ε] con ε > 0 muy peque˜ no, la ED es casi de orden uno...y para t grande es indudablemente 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)) Nota 4.4 Hemos sido consecuentes con nuestra notaci´on habitual escribiendo X en lugar de x. 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 C1 debe ser cero. Por lo tanto la soluci´on 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´on peri´odica f (t) con per´ıodo T . Z T e−st f (t)dt . f˜(s) = 0 1 − e−T s A continuaci´on la escribimos de nuevo utilizando ahora la notaci´on MATLAB. 1/(1-exp(-T*s))*(int(exp(-s*t)*f(t),0,T) Ejemplo 4.14 Obtener la transformada de la funci´on f (t), peri´odica de per´ıodo 2π, definida por  0≤t≤π  20 f (t) =  −20 π < 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

La funci´ on de transferencia H(s) de un sistema lineal se define como la raz´on entre la transformada 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´a gobernado por la EDL de coeficientes constantes siguiente ay 00 + by 0 + cy = g(t), t > 0, (4.8) y(0) = y 0 (0) = 0. Al aplicar transformadas en ambos miembros de la ecuaci´on, 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´ormula para H(s) en t´erminos del polinomio caracter´ıstico de la ED H(s) =

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

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

as2

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

(4.9)

Notar que H(s) no depende de la entrada g(t). Definamos ahora la llamada funci´ on de respuesta al impulso como 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´a en capacidad de deducir de (4.9) que la soluci´on y(t) de la ecuaci´on (4.8) est´a dada por (aplicar teorema de convoluci´on)
y(t) = g ∗ h (t) (4.10) Lo anterior significa que h(t) juega el papel de n´ ucleo resolvente de la ecuaci´on (4.8). Se halla de una vez, y luego, convolucionado con cualquier g(t) de entrada nos produce la respuesta o salida y(t). Ejercicio 4.6 Mediante c´alculo directo halle la funci´on de transferencia H(s) del sistema. Con la ayuda del MATLAB encuentre la funci´on de respuesta al impulso h(t) de los siguientes PVI, con condiciones nulas. (i) y 00 + 9y = g(t) (iii) y 00 − y 0 − 6y = g(t) (v) y 00 + 2y 0 + 5y = g(t)

(ii) y 00 − 9y = g(t) (iv) y 00 + 2y 0 − 15y = g(t) (vi) y 00 − 4y 0 + 5y = g(t)

Los resultados que se piden ya sabemos que no dependen de g(t). Usar la f´ormula (4.10) para calcular la soluci´on y(t) de (i) cuando g(t) = 4t2 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´odigo 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 sen(t) =

eit − e−it eit + e−it , cos(t) = 2i 2

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

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

OBJETIVOS Presentar los m´etodos num´ericos de Euler y Euler mejorado para aproximar la soluci´on de problemas de valores iniciales. Resolver problemas concretos, obteniendo soluciones num´ericas y estimados del error, mediante el manejo de programas dom´esticos MATLAB basados en los algoritmos correspondientes.

5.1

FUNCIONES ESPECIALES

En Matem´atica 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 00 + xy 0 + (x2 − q 2 )y = 0

(5.1)

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

∞ X n=0

(−1)n x ( )2n−q n!Γ(1 − q + n) 2

que son las funciones de Bessel de primera especie de ´ordenes q y −q respectivamente. Las Funciones de Bessel verifican diversas relaciones de recurrencia, a saber: Jq+1 (x) + Jq−1 (x) = 2q/xJq (x), Jq+1 (x) − Jq−1 (x) = 2Jq0 (x).

65

Clase Pr´actica 5

66

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

5.1.1

Funciones de Bessel

Dentro del grupo de funciones matem´aticas especializadas que contiene MATLAB se encuentran las funciones de Bessel de primera clase: besselj(n, x). A continuaci´on se representan la funci´ones de Bessel J1 y J0 . >>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´afica de J1

Fig. 5.2: Gr´afica de J0

Funciones de Bessel J0 y J1 Se representan conjuntamente en [−20, 20] las dos funciones de Bessel J0 y J1 . Obs´ervese c´omo se alternan los ceros de J0 y J1 .

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 J1/2 y J3/2 Ejecutar en la l´ınea de comandos el c´odigo >>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´afica.

Fig. 5.3: Gr´aficas superpuestas de J1 y J0

5.2

Fig. 5.4: Gr´aficas de Bessel J3/2 y J1/2

ESTUDIO DE LAS SOLUCIONES DE UNA ED

Los m´etodos a utilizar se pueden clasificar como: (1) M´etodos anal´ıticos de cuadraturas exactas (Temas 2, 3, 4 y 5). (2) M´etodos cualitativos. (3) M´etodos num´ericos (Clases Pr´acticas 5, 6 y 7). Los m´etodos 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´ormulas

Clase Pr´actica 5

68

impl´ıcitas o expl´ıcitas.(1 Para compensar estas limitaciones se introducen los m´etodos num´ericos que son el objeto de esta pr´actica.

5.2.1

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

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

(5.2)

siendo x0 ∈ [a, b] y f ∈ C([a, b] × R). ∂f son continuas en un rect´angulo que contiene al punto (x0 , y0 ), de modo que Se asume que f y ∂y el problema (5.2) tiene soluci´on u ´nica. Estrategia de los m´ etodos num´ ericos: Si y(x) es la soluci´on exacta del problema (5.2), se trata de obtener valores aproximados yk de y en ciertos puntos xk 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´on num´erica consiste en determinar una partici´on del subintervalo [x0 , X], esto es, x0 < x1 < · · · < xN = X. Se supone que los puntos est´an distribuidos uniformemente, de modo que xk+1 − xk = h, k = 0, 1, ..., N − 1. Los m´etodos num´ericos que estudiaremos permiten determinar los valores yk que constituyen una aproximaci´on del valor de la soluci´on de (5.2) en los puntos xk . Definici´ on 5.1 El conjunto {yk }k=0,...N se llama Soluci´on Num´erica 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´a en la curva soluci´on del mismo y que la pendiente de la tangente a dicha curva en P es y 0 (x0 ) = f (x0 , y0 ). El m´etodo de Euler aproxima el valor de la soluci´on exacta en x cercano a x0 usando la ecuaci´on de la tangente en el punto P = (x0 , y0 ): y(x) = y0 + f (x0 , y0 )(x − x0 ) 1

(5.3)

No se pueden resolver de forma exacta porque el cat´alogo disponible de “funciones elementales” no alcanza para obtener una f´ ormula cerrada, o porque dicha expresi´on 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 ), · · · , y(xN ), entonces, aplicando (5.3) se sigue que dado y0 se calcula y1 mediante y1 = y0 + f (x0 , y0 )h, e y2 mediante y2 = y1 + f (x1 , y1 )h, y asi sucesivamente hasta yN yN = yN −1 + f (xN −1 , yN −1 )h, En forma compacta se obtiene el siguiente esquema yk+1 = yk + hf (xk , yk ), k = 0, · · · , N − 1.

5.2.3

(5.4)

Otra interpretaci´ on del m´ etodo de Euler

Integrando ambas partes de la ecuaci´on y 0 (x) = f (x, y(x)) en el intervalo [xk , xk+1 ], obtenemos Z

xk+1

y(xk+1 ) − y(xk ) =

0

Z

xk+1

f (x, y(x))dx,

y (x)dx = xk

xk

y sustituyendo la integral de la derecha por el ´area de cierto rect´angulo Q se obtiene que y(xk+1 ) ≈ y(xk ) + f (xk , y(xk ))(xk+1 − xk ). Si ponemos el paso h = xk+1 − xk y sustituimos y(xk ) por yk obtenemos de nuevo el esquema (5.4) del m´etodo de Euler. yk+1 = yk + f (xk , yk )h, k = 0, · · · , N − 1. (5.5) La f´ormula (5.5), llamada “ecuaci´on en diferencias finitas”, permite predecir yk+1 en t´erminos del valor de yk . El dato o condici´on inicial (CI) es y0 y desde ´el marchamos en una direcci´on haciendo cada vez una nueva estimaci´on sobre la base de predicciones anteriores. Esto u ´ltimo introduce desviaciones en nuestros sucesivos c´alculos, provocando que ´estos sean m´as imprecisos a medida que nos alejamos del punto inicial x0 . ´ La formulaci´on (5.5) anterior corresponde al llamado m´etodo de Euler. Este resulta ser el m´as sencillo entre los m´etodos 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´on profesional es nula. Sea y(x; h) la soluci´on aproximada obtenida con paso h, mediante un m´etodo num´erico, por ejemplo Euler, de un PVI cuya soluci´on 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´e significa decir que el m´etodo 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´etodo converge. Lo anterior significa que las soluciones y(x; h) de las ecuaciones en diferencias finitas que caracterizan al m´etodo considerado, aproximan a la soluci´on exacta y(x) de la ED cuando h → 0. Luego, te´oricamente, deber´ıamos obtener mejores respuestas considerando valores cada vez m´as peque˜ nos de h. Desafortunadamente, la aritm´etica de precisi´on finita del ordenador no lo permite. Para h muy peque˜ no los valores yk 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´on exacta y(x) se expresa de la siguiente forma: y(x) − y(x; h) = O(hp ), lo que significa que |y(x) − y(x; h)| ≤ M |h|p para h peque˜ no y M > 0 dependiendo de x. Definici´ on 5.3 Al exponente p de la relaci´on anterior se le llama orden del m´etodo. Nota 5.1 El m´etodo de Euler es de orden p = 1.

5.2.4

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

Ejemplo 5.1 Obtener una soluci´on aproximada del PVI y 0 = y 2 + 2x − x4 , 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 0 = y 2 + 2x − x4 , a ≤ x ≤ b, con y(a) = y0 . Se calcula tambi´en el error cometido comparando con la soluci´on exacta del PVI que es y = x2 . 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´ on de punto flotante de y(xk ).

Clase Pr´actica 5

Fig. 5.5: M´etodo de Euler

71

Fig. 5.6: Interpretaci´on del m´etodo de Euler

Resoluci´ on del ejemplo 5.1 Programa EULER1.m Datos: n es el n´ umero 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´aximo 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´aximo es: 8.093460004977e − 02 Los valores aproximados se obtienen en la l´ınea de comandos seg´ un se indica a continuaci´on. >>[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 ´este 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

Ejercicio 5.1 Modificar el programa EULER1 para calcular una soluci´on aproximada del PVI y 0 = 2xy 2 , 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´on se presenta el siguiente c´odigo modificado para el m´etodo de Euler que llamaremos EULER2.m. para resolver la ecuaci´on y 0 = 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´on F se define previamente mediante el comando inline. Se obtiene el gr´afico y el vector y de la soluci´on. El programa EULER2.m es un procedimiento de tipo FUNCTION y uno de sus argumentos es una funci´on que se debe crear previamente. Una de las aplicaciones m´as u ´tiles de los ficheros .M se presenta en la definici´on de funciones mediante la directiva FUNCTION cuya sintaxis es FUNCION PARAMETROS SALIDA=NOMBRE FUNCION(parametros entrada) El c´odigo 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)’) Ejemplo 5.3 Calcular aproximadamente, usando EULER2.m, la soluci´on de la ED y 0 = y + sen(x) con y(0) = 1, 0 ≤ x ≤ 1.

Fig. 5.9: Soluci´on 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´ervese que el programa EULER2 no permite juzgar la calidad de los resultados obtenidos. Ejercicio 5.2 Aproximar, usando EULER2.m, la soluci´on de la ED 1 (y + y 2 ), x y(1) = 1, 1 ≤ x ≤ 2 y0 =

(5.6)

tomando el paso h = 2−2 .

5.2.5

Precisi´ on de los resultados

Todo c´alculo aproximado debe tener en cuenta un indicador que permita evaluar la calidad de los resultados. En nuestro caso se trata de calcular la soluci´on aproximada de una ED y de poder decidir sobre la precisi´on del resultado. El error exacto es imposible de conocer, de modo que se trata de obtener una estimaci´on del mismo. Si dicha estimaci´on se obtiene despu´es de realizar los c´alculos y se basa en la utilizaci´on de ´estos, se dice que es una estimaci´on del error a posteriori. El siguiente programa EULER3 es una versi´on m´as 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´aficas de dos resultados calculados con diferentes mallas de puntos, lo que permite estimar el error a posteriori si nos basamos en el siguiente principio: Si d es la distancia entre dos aproximaciones de la soluci´on exacta de un PVI, entonces d 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´on que se debe crear previamente. 4

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

Clase Pr´actica 5

76

Resuelve la ecuaci´on y 0 = F (x, y) con la condici´on inicial y(a) = y0 , 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´on F se define previamente mediante el comando INLINE. Procedimiento para estimar el error Se estima a posteriori el error tomando la diferencia m´axima 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´on del error a posteriori y los gr´aficos de ambas aproximaciones num´ericas. C´odigo 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´on 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´on 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 0 = y 2 + 2x − x4 , 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 0 = −y 2 , 0 ≤ x ≤ 1/2, con y(0) = −1.

(5.7)

• ¿A qu´e tipos conocidos de ecuaciones corresponde la anterior ED? • ¿Por qu´e no es lineal? • ¿Qu´e m´etodos utilizar´ıas para intentar hallar una soluci´on exacta? • Usando EULER3 obtener las soluciones aproximadas con N = 23 y N = 24 − 1 puntos y estimar el error cometido. • 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´etodo num´erico m´as eficiente que el de Euler, utilizaremos una nueva f´ormula de integraci´on num´erica de mayor orden de exactitud que la utilizada en el m´etodo de Euler.

5.3.1

La f´ ormula de los trapecios

La f´ormula de los trapecios consiste en sustituir el valor exacto de la integral b

Z

f (x)dx, a

por el n´ umero (b − a)

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

es decir, el a´rea del trapecio cuyos v´ertices son (a, 0), (a, f (a)), (b, f (b)) y (b, 0). Con ello se trata de efectuar c´alculos m´as precisos en la aproximaci´on de la integral sobre cada uno de los subintervalos [xk , xk+1 ] Z 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´ognita y(x) de nuestro problema.

Fig. 5.11: Regla del trapecio En resumen y dicho en t´erminos geom´etricos, el m´etodo de Euler consisti´o en sustituir el a´rea bajo la curva de g(x) = f (x, y(x)) por el ´area de cierto rect´angulo, mientras que ahora, en lugar del

Clase Pr´actica 5

79

rect´angulo utilizaremos un trapecio para obtener la siguiente expresi´on(5 Z xk+1  h f (x, y(x))dx ≈ f (xk , y(xk )) + f (xk+1 , y(xk+1 )) , 2 xk donde h = xk+1 − xk . Surge as´ı un m´etodo de c´alculo aproximado cuya formulaci´on es la siguiente. yk+1 = yk +

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

Se trata de un m´etodo impl´ıcito ya que la inc´ognita yk+1 , no aparece despejada en la f´ormula. Adem´as se prueba que el error global de discretizaci´on E(h) cumple E(h) = O(h2 ), y por tanto es de orden p = 2.

5.3.2

Deducci´ on del m´ etodo de Euler mejorado

El m´etodo de Euler mejorado se obtiene realmente como el m´etodo predictor-corrector siguiente (1) Se predice yk+1 mediante Euler ∗ yk+1 = yk + hf (xk , yk ).

(2) Se corrige mediante la anterior relaci´on   h yk+1 = yk + f (xk , yk ) + f (xk + h, yk + hf (xk , yk )) | {z } | {z } 2 ∗ xk+1 yk+1 Resultando definitivamente el algoritmo de Euler mejorado yk+1 = yk +

5.3.3

h [f (xk , yk ) + f (xk+1 , yk + hf (xk , yk ))] . 2

Euler versus Euler mejorado

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

Ver secci´ on 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´on y 0 = F (x, y) con la C.I. y(a) = y0 , en el intervalo [a, b], en n puntos x(k), (k = 1, · · · , n), con el paso h = (b − a)/(n − 1), utilizando los m´etodos Euler y Euler mejorado. Las funciones F (campo de direcciones) y f (soluci´on exacta conocida) se definen previamente mediante el comando INLINE. El programa estima los respectivos errores absolutos y los ofrece como salida, seg´ un el formato que m´as abajo se explica. Estos errores se obtienen directamente en t´erminos del m´aximo de la distancias entre cada estimaci´on yek y el valor exacto f (xk ). Los correspondientes gr´aficos se despliegan en una misma ventana para facilitar la apreciaci´on 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’) Ejemplo 5.5 Obtener una soluci´on aproximada de la ED y 0 + y/x − y 2 log(x) = 0, en [1, 3] con la C.I. y(1) = 1, usando los m´etodos de Euler y de Euler mejorado. Para el c´alculo del error, usar que la soluci´on exacta es f (x) = −2/x/(log(x)2 − 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´on 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´etodo Euler mejorado comparando ´este 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´etodo de Euler y Euler mejorado, es y em respectivamente, para el problema (a) con d = 1 y 24 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´on de este tipo de problemas debe permitirnos llegar a conclusiones experimentales y te´oricas, y adem´as, ganar experiencia en la aplicaci´on de estos m´etodos.

5.4

´ ´ OTROS METODOS NUMERICOS

Podemos dise˜ nar otros m´etodos num´ericos para resolver y 0 = f (x, y), y(x0 ) = y0 , calculando aproximadamente la integral Z xk+1 f (x, y(x))dx, xk

mediante nuevos procedimientos num´ericos. Seg´ un sea el orden del error producido por la f´ormula de integraci´on que utilicemos, as´ı ser´a el orden del m´etodo num´erico que permitir´a calcular aproximadamente la soluci´on del PVI.

5.4.1

Ejemplo de m´ etodo de dos pasos

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

(5.8)

Clase Pr´actica 5

83

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

5.5

CONCLUSIONES

En general podemos postular que 1. El orden de un m´etodo es un indicador que est´a relacionado con el n´ umero de operaciones l´ogico-aritm´eticas que han de ejecutarse para alcanzar una precisi´on dada. Es una referencia vaga ya que en la pr´actica experimental influyen diversos factores en la eficiencia de un m´etodo. No obstante, permite hacer una primera comparaci´on entre dos m´etodos diferentes. 2. No existe un m´etodo num´erico 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´etodo a un problema concreto cuya soluci´on exacta es desconocida. 4. El m´etodo 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´etodos num´ericos que resuelven aproximadamente problemas de valor inicial. El m´etodo de Euler pertenece a una clase de m´etodos llamados de “un paso”, y que en general se caracterizan por la siguiente relaci´on recurrente: yk+1 = yk + φ(h, xk , yk ).

Clase Pr´actica 5

84

5.6

EJERCICIOS ADICIONALES

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

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

(5.9)

Detectar una anomal´ıa, en este caso la presencia de una “singularidad” en x = 3, calculando la soluci´on 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 u ´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´on exacta del problema es y = 9/(9 − x2 ).

5.6.1

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´on 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. Ejercicio 5.6 Utilizar el c´odigo EULER3 para obtener soluciones aproximadas de xy 0 = y b(x), y(1) = 1, en los dominios siguientes (I) 1 ≤ x ≤ 2, (II) 1 ≤ x ≤ 2 × 102 , (III) 1 ≤ x ≤ 2 × 1012 , 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´erica? 3. ¿Se justifica aplicar un m´etodo num´erico cuando b(x) es constante (apartado 1)?

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

OBJETIVOS 1. Aplicar los m´etodos num´ericos 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´aficos. En esta clase el alumno comprobar´a la importancia que tiene transformar ecuaciones y sistemas diferenciales de orden superior en sistemas de primer orden. Al igual que la clase anterior, s´olo utilizaremos herramientas gr´aficas y num´ericas 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´ognitas x = x(t) e y = y(t).  0  x = f (t, x, y) y 0 = g(t, x, y) S:  x(t0 ) = x0 , y(t0 ) = y0 Notemos que el anterior sistema S puede ser no-lineal. Si los t´erminos de la derecha no dependiesen de t dir´ıamos que S es un “sistema aut´onomo”.

6.1.1

M´ etodo de Euler mejorado

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

Euler mejorado pertenece a la familia de m´etodos de Runge-Kutta que m´as adelante estudiaremos.

87

Clase Pr´actica 6

88

1. En cada etapa debe ser aplicado el m´etodo a todas las ecuaciones antes de pasar a la siguiente etapa. Adem´as, 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´on num´erica de una ecuaci´on o sistema diferencial de orden mayor que uno debe ser primero transformado en un sistema de orden uno. A ´este u ´ltimo se le aplica directamente el m´etodo num´erico. 3. Cualquiera de los m´etodos de diferencias finitas que estudiemos para resolver el PVI escalar y 0 = F (x, y), y(x0 ) = y0 , puede tambi´en ser aplicado a un sistema diferencial. A continuaci´on recordamos c´omo se transforma una ED de orden dos y 00 = f (t, y 0 , y)

(6.1)

en un sistema de EEDD de orden uno. Si hacemos x = y 0 , entonces y 00 = x0 . Sustituyendo en (6.1) y agregando la nueva ecuaci´on y 0 = x, originada por el propio cambio, obtenemos el sistema de primer orden x0 = f (t, x, y) y0 = x De su soluci´on (x(t), y(t)) s´olo necesitamos la componente y = y(t) que es la soluci´on del problema original. El algoritmo de Euler mejorado aplicado a un sistema diferencial gen´erico S queda formulado tal como se muestra a continuaci´on Euler mejorado aplicado a un sistema xk+1 = xk +


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

(6.2)


h ∗ yk+1 = yk + g(tk , xk , yk ) + g(t∗k+1 , x∗k+1 , yk+1 ) 2 donde las predicciones para k + 1 est´an dadas por t∗k+1 = tk + h x∗k+1 = xk + hf (tk , xk , yk )

(6.3)

∗ yk+1 = yk + hg(tk , xk , yk )

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´odulos. El primero, a modo de interfaz logra que el c´odigo sea transparente al usuario. El alumno s´olo tiene que saber qu´e datos debe dar como entrada, e interpretar los resultados.

Clase Pr´actica 6

89

Los datos relativos al campo de direcciones o a la entrada opcional de la soluci´on exacta se entregan cuando lo solicita el programa, con la sintaxis o notaci´on usual, sin encerrar entre ap´ostrofes ni poner puntos delante de operadores. Esto u ´ltimo es necesario para la ejecuci´on pero el propio programa se encarga de hacerlo autom´aticamente. Como sintaxis se admite una cualquiera de las siguientes:(2 s2_euler y=s2_euler [t,y]=s2_euler [t,x,y]=s2_euler El segundo m´odulo est´a constituido por dos sub-programas situados al final del c´odigo. Uno de ellos, tal como se se˜ nal´o anteriormente, se encarga de “escribir” el punto delante del operador divisi´on \, la multiplicaci´on *, y la exponenciaci´on ^. La otra sub-funci´on, de mayor inter´es para el alumno, permite simular el m´etodo de Euler mejorado aplicado a un sistema diferencial de dos ecuaciones, y se corresponde con la formulaci´on (6.2-6.3). Los c´alculos se realizan atendiendo al paso h dado por el usuario, y para h/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´aficos de las soluciones aproximadas aparecen en pantalla superpuestos y s´olo corresponden al u ´ltimo paso. El programa RK CELL, que tambi´en podemos obtener en la Web de la 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 de cualquier n´ umero de dimensiones. El propio programa crea autom´aticamente el fichero FUNCTION correspondiente al campo vectorial, con las caracter´ısticas exigidas por ODE45. Este comando Matlab ejecuta algoritmos correspondientes a los m´etodos de Runge-Kutta que estudiaremos en la clase pr´actica 7.

6.1.3

Ejercicios con indicaciones para ser resueltos

Ejercicio 6.1 Dado el sistema diferencial lineal siguiente x0 = x + y, y 0 = x − y, √ x(0) = 1 + 2, y(0) = 1,

(6.4)

donde 0 ≤ t ≤ 1, emplear el programa S2 EULER con los pasos h = 0.05/2k , k = 0, 1, 2, 3, 4, para resolverlo. Indicaciones para el ejercicio 6.1 Escribir 2

Para m´ as informaci´ on ejecutar help s2 euler

Clase Pr´actica 6

90

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

√ √ t√2 2)e , y(t) = et 2

El problema (6.4) es lineal y homog´eneo. El alumno deber´ıa intentar resolverlo artesanalmente mediante los m´etodos estudiados en clases. Ver la figura 6.1.

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

Ejercicio 6.2 Resolver num´ericamente el sistema no lineal: x0 = x2 + x + 5y y 0 = y 2 + sen(x) + 3x − y x(0) = 0.001, y(0) = 0.001,

0 ≤ t ≤ 1.

Obtener una estimaci´on del error para h = 0.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´es obtener los resultados. No cerramos el gr´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.

Clase Pr´actica 6

91

Fig. 6.2: Gr´afico producido por S2 EULER al resolver el ejercicio 6.2 Ejercicio 6.3 Resolver num´ericamente el sistema lineal: x0 = x + 5y y 0 = 4x − y x(0) = 0.001, y(0) = 0.001,

0 ≤ t ≤ 1.

Obtener una estimaci´on del error para h = 0.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´afico anterior, del ejercicio 6.2, quede superpuesto al gr´afico del ejercicio 6.3. A continuaci´ on se muestra c´ omo se ha “perturbado localmente” el problema 6.2 para transformarle en el problema 6.3 Los problemas planteados en los ejercicios 6.2 y 6.3 est´an relacionados entre 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. Mostraremos c´omo el campo de direcciones del sistema diferencial del ejercicio 6.2 se aproxima localmente mediante su diferencial, y c´omo a partir de ´este u ´ltimo se producen las ecuaciones lineales del problema 6.3. Se dice que el sistema del ejercicio 6.3 es una aproximaci´on lineal del sistema del ejercicio 6.2. Las expectativas son: que la soluci´on de 6.3 sea una aproximaci´on aceptable de la correspondiente a 6.2 en un intervalo [0, ε], con ε > 0 “suficientemente peque˜ no”, y que el m´etodo num´erico seleccionado sea m´as efectivo resolviendo el problema lineal 6.3, que el problema no lineal 6.2.(3 3

La frase “suficientemente peque˜ no” encierra otro problema cuya soluci´on no abordaremos en este material.

Clase Pr´actica 6

92

Fig. 6.3: Gr´afico producido por S2 EULER al resolver el ejercicio 6.3 El c´alculo de los siguientes diferenciales en (0, 0) D(x2 + x + 5y)(0, 0) = dx + 5dy D(y 2 + sen(x) + 3x − y)(0, 0) = 4dx − dy y la conocida f´ormula f (x, y) ≈ f (0, 0) + Df (0, 0)(x − 0, y − 0) v´alida para x e y “peque˜ nos”, nos muestran que x2 + x + 5y ≈ x + 5y y 2 + 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´enea de orden dos con coeficientes variables. y 00 − ty 0 + y = 0, y(0) = 2, y 0 (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´etodo 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´afico obtenido al resolver el ejemplo 6.1 Como comprobaci´on de los resultados apreciaremos que la gr´afica de la soluci´on num´erica, obtenida seg´ un Euler, es muy similar a la del polinomio y(t) = 2 + t − t2 −

t6 t4 − , 12 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´on de la soluci´on exacta y los gr´aficos de las soluciones aproximadas. Comencemos viendo que la soluci´on exacta del problema (6.5), expresada en t´erminos de una serie de potencias en torno al punto x = 0, es ∞ X

an tn ,

n=0

siendo los an definidos por la recurrencia a tres t´erminos an+2 =

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

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

Observar que siendo t ∈ [0, 1], y atendiendo a la precisi´on deseada, el t´ermino con la potencia t6 tambi´en podr´ıa despreciarse en la aproximaci´ on (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´on de la soluci´on del problema (6.5) 6 X

an tn = a6 t6 + · · · + a1 t + a0 .

n=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´etodo Euler mejorado aplicado al sistema asociado, ejecutar S2 EULER optando por la salida gr´afica 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, ´estos deben ser convertidos previamente en un sistema de orden uno, y aplicarle a este u ´ltimo uno de los m´etodos estudiados para el caso escalar. La serie de potencias que representa a la soluci´on de un problema con datos anal´ıticos, puede ser simulada num´ericamente utilizando una suma parcial seleccionada convenientemente, lo que, en principio, constituye otro m´etodo para aproximar la soluci´on.

6.2

EJERCICIOS ADICIONALES

Ejercicio 6.4 Aproximar los valores de la soluci´on del PVI y 00 + 2t2 y 0 + 5y = 0, y(0) = 1, y 0 (0) = −1, t ∈ [0, 2]. en los puntos t = 0.21, 0.95, 1.31, 1.84. Ejercicio 6.5 Aproximar los valores de la inc´ognita y = y(t) x0 = tx + 5y + t2 , y 0 = −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

95

Ejercicio 6.6 x0 = tx − 10y − 4t − 2t2 , y 0 = x + 6y, x(0) = 1, y(0) = −1, t ∈ [0, 1]. Ejercicio 6.7 y 0 = x + 6y + 3, x0 = 2x + y, x(0) = 2, y(0) = 4, t ∈ [0, 2].

6.2.1

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´on por la v´ıa num´erica. Para ello escribimos la EDLH de 2do orden del ejercicio 6.4 como un sistema de primer orden x0 = −2t2 x − 5y, y 0 = x, x(0) = −1, y(0) = 1, para aplicar S2 EULER.M o RK CELL.M. Hagamos >>[t,y]=s2_euler porque s´olo nos interesa almacenar en memoria los valores aproximados de y = y(t), la inc´ognita del problema original. Entrar h = 0.2 Desplegamos el gr´afico para contrastarlo con el de la figura 6.5. Los puntos t donde se nos pide el valor aproximado de la soluci´on no est´a 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´on aparecen las listas de los puntos interpolados (T1 , Y1 ) y el gr´afico que muestra su posici´on (ver figura 6.5). La leyenda que aparece en la esquina superior izquierda puede moverse arrastr´andola con el cursor.

Clase Pr´actica 6

96

Fig. 6.5: Gr´aficas 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´ermino tx de la primera ecuaci´on provoque una complejidad excesiva o no existe, de hecho, una soluci´on representable con el cat´alogo de funciones MATLAB. Nos queda la alternativa de aplicar S2 EULER o RK CELL al sistema anterior. Los resultados gr´aficos que se muestran en la figura 6.6 corresponden al m´etodo 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

Clase Pr´actica 6

97

>>hold on Podemos comprobar experimentalmente la coincidencia de dos m´etodos vistos en clases. Si transformamos el sistema del problema 6.6 en una EDL de 2do orden obtenemos y 00 − (t + 6)y 0 + (6t + 10)y = −4t − 2t2 , y(0) = −1, y 0 (0) = −5. La anterior ecuaci´on es excesivamente complicada para el n´ ucleo simb´olico de MATLAB. La ecuaci´on homog´enea asociada tiene soluci´on mediante el comando DSOLVE en t´erminos de la funci´on especial WhittakerM.(6 Notemos que la ED lineal y 00 − (t + 6)y 0 + (6t + 10)y = −4t − 2t2 , puede ser resuelta mediante series de potencias aunque su proceso de resoluci´on es algo m´as complicado que los problemas vistos en clases. Sin embargo, lo aprendido nos sobra para abordar su soluci´on en torno a t = 0. Despu´es de asumir que y(t) es anal´ıtica en x0 = 0, llegamos a que

y(t) =

∞ X

an tn

n=0

donde

Fig. 6.7: Soluci´on gr´afica del ejercicio 6.6 6

Usar mhelp en lugar de help para ver m´as detalles sobre esta funci´on especial.

Clase Pr´actica 6

98

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

−1, −5, 3a1 − 5a0 = −10, (12a2 − 9a1 − 6a0 − 4)/6, (18a3 − 8a2 − 6a1 − 2)/12, ··· 6(n + 1)an+1 + (n − 10)an − 6an−1 , n = 3, 4, 5, ... = (n + 2)(n + 1)

= = = = =

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´ermino N = 12, procedemos como sigue. >>t=linspace(0,1,100); >>yy=truncada(t) >>hold on >>plot(t,yy,’+’) En la figura 6.7 se comprueba gr´aficamente el nivel de coincidencia entre las dos soluciones aproximadas. Indicaciones para el ejercicio 6.7 El ejercicio 6.7 tiene soluci´on calculable con DSOLVE pero con una expresi´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 las obtenidas usando S2 EULER o RK CELL. Para hacer los gr´aficos se recomienda usar EZPLOT. 7

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

Clase Pr´actica 6

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(S.x,[0,2])%cambiar el color de la grafica >>hold on >>ezplot(S.y,[0,2]) >>title(’ ’)% titulo vacio >>legend(’S.x’,’S.y’)

6.3 6.3.1

SERIES DE FOURIER Programas

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

6.3.2

Ejercicios

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

−1 si −π ≤ x < 0, 1 si 0 ≤ x < π

(6.7)

representar gr´aficamente la funci´on f y las aproximaciones de su serie de Fourier trigonom´etrica mediante las primeras sumas parciales. Indicaciones para resolver el ejercicio 6.8 La serie de Fourier de f viene dada por ∞ 4 X sen((2n − 1)x) . π n=1 2n − 1

Basta ejecutar los programas SERIE y SERIE1. Ejercicio 6.9 Representar la funci´on 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

∞ X sen(nx) 2 (−1)n+1 . n n=1

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´ umero arbitrario de ecuaciones, utilizando el m´etodo Euler mejorado. Este c´odigo 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´erminos vectoriales. La otra subrutina FUNCTION, llamada CEVAL, permite la evaluaci´on num´erica 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´ognitas escribimos >>x=sn_euler; Si adem´as 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´ umero de filas de x se corresponde con el n´ umero de puntos t0 ,...,tm . Para visualizar una breve descripci´on 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 0 = F (x, y), y(x0 ) = y0 . Una vez calculados los valores aproximados de la soluci´on y(x) de un PVI, sean estos yk ≈ y(xk ), k = 1, ..., n nos puede interesar aproximar el valor de y(t) para t ∈ (xk−1 , xk ). Seguiremos el procedimiento m´as sencillo consistente en “interpolar linealmente”, es decir uniremos los puntos Pk−1 = (xk−1 , yk−1 ), Pk = (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´etodos 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´on mediante poligonales. Para realizar eficientemente los c´alculos 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´etodos de diferencias finitas nos han ofrecido la soluci´on como un vector de Rn , y que seg´ un el procedimiento basado en interpolaci´on estamos hallando una soluci´on en la forma de una funci´on que, aunque continua, desafortunadamente no es derivable en todos los puntos. A continuaci´on se muestra el c´odigo de INTERPOLA.

Clase Pr´actica 6

103

Fig. 6.11: Ilustraci´on gr´afica del m´etodo de interpolaci´on function interpola disp(’Se exige que x1>[T,Y]=rungekutta

dependiendo de la informaci´on 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´on. En cualquier caso el programa permite opcionalmente ver en la pantalla los valores T (k) e Y (k) (hasta 50 como m´aximo), pero no los almacena en memoria.

7.1.4

Ejercicios

Ejercicio 7.1 Aproximar la soluci´on del PVI y0 = 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. Ejercicio 7.2 Aproximar la soluci´on del PVI 2y 0 = (x − y), y(0) = 1 para x ∈ [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 − 2 + 3e−x/2 calcular el error exacto cometido y compararlo con el estimado a posteriori.

Clase Pr´actica 7

108 Ejercicio 7.3 Aproximar la soluci´on del PVI (ex + 1)y 0 + y(ex − 1) = 0, y(0) = 3 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 =

12ex calcular el error exacto cometido y (ex + 1)2

compararlo con el estimado a posteriori.

7.1.5

Estabilidad de un m´ etodo num´ erico

Un m´etodo num´erico para resolver una EDO de primer orden, est´a dado por una ecuaci´on en diferencias finitas (EDF). Decimos que el m´etodo es estable si su EDF produce una soluci´on acotada cuando la soluci´on exacta es acotada. Si la soluci´on obtenida no es acotada, siendo la exacta acotada, decimos que el m´etodo es inestable. S´olo analizamos la estabilidad cuando tratamos con EEDD ordinarias lineales y 0 + py = q(x), p ∈ R, en las que el t´ermino q(x) no influye. Los siguientes enunciados pueden verse en Hoffman [3], secciones 7.8-7.9. • Los m´etodos de Euler y Euler mejorado son estables cuando 0 < ph ≤ 2. • Runge-Kutta es estable cuando 0 < ph ≤ 2.784. Ejercicio 7.4 Estudiar num´ericamente la soluci´on del siguiente problema de valor inicial (Hoffman [3], 7.14) aplicando el programa EJE7.4 (Web de la asignatura [15]). y 0 = p(y − (x + 2)) + 1, y(0) = 1, x ∈ [0, +∞). Considere los valores p = ±1000 considerando los intervalos [0, b], b > 0, y los pasos h seg´ un se indica en la tabla 7.1.5. 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´enea asociada. Seg´ un se aprecia en la tabla 7.1.5, el tama˜ no de h se selecciona m´as peque˜ no para mostrar: (A-B) c´omo cambia la precisi´on (p = 1000)

Clase Pr´actica 7

109

b A 0.01 B 0.01 C 0.01 D 1 E 1

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

RK-4 ODE45

Tabla 7.1: Informe de resultados seg´un los valores de b, h y p

(D-E) estabilidad condicional (p = −1000) (Se detectan en los valores num´ericos del error) (ODE45 no se afecta pues selecciona autom´aticamente 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´aficamente 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

hE hA y que hD = . 2 5

Descripci´ on del programa EJE7.4

Tal como sugiere su nombre, este c´odigo est´a especialmente dise˜ nado para resolver el ejercicio 7.4 de esta pr´actica. Para facilitar el trabajo se cuenta con una interfaz de usuario que permite seleccionar uno de los dos posibles PVI mediante una ventana-men´ u, atendiendo a que la ecuaci´on que le define depende del valor de un par´ametro p. El programa se encarga de aplicar los m´etodos de Euler, Euler mejorado y Runge-Kutta 4, mostrando al final los respectivos errores y los gr´aficos (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´aximo permitido), al menor de los errores producidos por los m´etodos anteriores. Los principales objetivos de este experimento consisten en la detecci´on experimental de inestabilidad y escalas relevantes en el dominio.

Clase Pr´actica 7

110

Fig. 7.1: Gr´aficas 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 00 + b(x)y 0 + 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 0 (x), B2 [y](x) = β1 y(x) + β2 y 0 (x). Al siguiente problema Λ[y](x) = f (x), x ∈ [a, b], B1 [y](a) = 0, B2 [y](b) = 0,

(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´on asociado a Λ[y], es decir L[y] = w(x)Λ[y] = (p(x)y 0 )0 + 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´ uan longitudinalmente sobre una columna empotrada y las deformaciones correspondientes disponemos del modelo matem´atico siguiente λ y(x) = 0 y 00 (x) + 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 00 (x) + y(x) = 0 y diferentes condiciones homog´eneas de frontera, para mostrar gr´aficamente las EI deformaciones (autofunciones) y calcular las cargas que producen deformaci´on (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

Ejemplo 7.2 Representar gr´aficamente la soluci´on del siguiente problema: 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´afica producida por CALOR.M Resoluci´ on del ejemplo 7.2 La soluci´on exacta u(x, t) de (7.5) es u(x, t) = 5 sen(4πx) exp(−5(122 )π 2 t/9) − 3 sen(8πx) exp(−5(242 )π 2 t/9) + 2 sen(10πx) exp(−5(302 )π 2 t/9), y se ha obtenido aplicando el m´etodo de separaci´on de variables. El programa CALOR (se descarga en [15]) utiliza la soluci´on de (7.5) para hacer los gr´aficos del ejemplo 7.2. Ejemplo 7.3 El flujo de calor a trav´es de una barra cil´ındrica de longitud π, es soluci´on de: ut = 2uxx ; 0 < x < π, t > 0, u(0, t) = u(π, t) = 0,  x 0 < x ≤ π/2, u(x, 0) = π−x π/2 < x < π. 1. Obtener una aproximaci´on de la soluci´on (serie (7.8)) en el instante t = 0, en x = π/2.

(7.6) (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´etodo de separaci´on de variables permite obtener la soluci´on como la suma de la serie u(x, t) =

∞ exp(−2(2k + 1)2 t) sen((2k + 1)x) 4X (−1)k π k=1 (2k + 1)2

(7.8)

Usar la funci´on TEMP.M para aproximar la temperatura u(x, t) cualquiera sea el instante t y la posici´on 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´odigos se muestran a continuaci´on. 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>x=pi/6;t=linspace(0,10,50);Z=eval(T);plot(t,Z) Apreciar el punto donde se sospecha se alcanza el m´aximo, 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´aximo. 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

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

(7.10)

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´afica 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´es 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´on del ejercicio puede obtenerse aplicando el m´etodo de separaci´on de variables y viene dada por   ∞ X
(2n + 1)x 8(−1)n y(x, t) = cos sen (2n + 1)t . (7.11) 2 π(2n + 1) 2 n=0 La serie (7.11) converge para todo x, t ∈ R. Un an´alisis 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˜ nado atendiendo al anterior an´alisis y permite aproximar hasta con dos d´ıgitos correctos a la serie (7.11). A continuaci´on se muestra el c´odigo 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´aficos 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´ericas 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´es del n´ ucleo MAPLE. >>maple(’pdesolve(ecuacion,func(var1,...,varn))’) resuelve la ecuaci´on en derivadas parciales identificada por ecuacion, para la funci´on f = f unc(x1 , ..., xn ) 2

La precisi´ on 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) fxx + 5fxy = 3 b) 3gx + 7gxy = xy c) htt − hxx = 0 d) uxx + 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) Ejercicio 7.7 Resolver la ecuaci´on en derivadas parciales del Ejemplo 7.5, apartado d), aplicando herramientas simb´olicas MAPLE-MATLAB. d) uxx + uyy = x + y 2

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

BIBLIOGRAF´IA Y SITIOS WEB [1] ATKINSON, K. E., Elementary numerical analysis, J. Wiley and Sons, 1999. [2] COOMBES, K. R. et al., Differential equation with MATLAB. John Wiley and Sons, 1999. [3] HOFFMAN, J. D., Numerical methods for engineers and scientists, McGraw-Hill, 1992. [4] MALEK-MADANI, R., Advanced Engineering Mathematics with Mathematica and MATLAB, vol I, Addison-Wesley, 1998. ´ ´ L. & ZARZO, A., Ecuaciones Diferenciales. Pro[5] MARCELLAN, F., CASASUS, blemas lineales y aplicaciones, McGraw-Hill, 1990. [6] MATHEWS J. H. y FINK, K. D., M´etodos num´ericos con MATLAB. Prentice Hall, 2000. [7] MATLAB, edici´ on estudiante. Prentice Hall, 1995. [8] NAGLE, R. Y SAFF,E. Fundamentals of differential equations. Benjamin Cummings, 1986. [9] NAGLE, R., SAFF, E. & SNIDER, A. D. Ecuaciones Diferenciales y Problemas con valores en la frontera. Pearson Educaci´on, M´exico, 2001 [10] HANSELMAN, D. y LITTLEFIELD, B., The Student edition of MATLAB: version 5 user’s guide, Prentice Hall, 1997. [11] ZILL, D. G. Ecuaciones diferenciales con aplicaciones de modelado, Internat. Thomsom Editores, 1997. [12] Web de MATLAB: http://www.mathworks.com/products/matlab/ [13] El rinc´on de MATLAB: http://matlab.universas.com/ [14] Wikipedia: http://es.wikipedia.org/wiki/MATLAB [15] Web de la asignatura: http://webs.uvigo.es/ecuacionesdiferenciales/

119