Octave Practicas

Octave Pr´ acticas con Octave

´ Indice general

1. Un lenguaje lenguaje de programac programaci´ i´ on on

4

1.1. Introducci´ Introducci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2. Operadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2.1. Operacionales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2. 1.2.2. 2. L´ ogicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3. Comandos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3.1. Generales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3.2. Aritméticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3.3. Algebraicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2. Matrices, submatrices, determinantes, sistemas de ecuaciones

9

2.1. 2.1. C´ alculo matricial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.1.1. Comandos generales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.1.2. Tipos de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.1.3. Operaciones con matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2. .2. Sist Sistem emaas de ecua ecuaccione iones. s. Método todoss dir directo ectoss . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.2.1. 2.2 .1. Soluci Soluci´ón o n de sist sistem emaas Com Compati patibl bles es Deter etermi mina nado doss . . . . . . . . . . . . . . . . . . .

12

2.2.2. 2.2 .2. Soluci Soluci´ón o n de sist sistem emaas Com Compati patibl bles es Inde Indete term rmin inad ados os . . . . . . . . . . . . . . . . . .

12

2.2.3. 2.2.3. Algoritmo Algoritmo de eliminaci´ eliminaci´ on de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.2. 2.2.44. Consi onside dera raccione ioness impo imporrtan tantes tes sobr sobree sis sistema temass . . . . . . . . . . . . . . . . . . . . . .

14

3. Ecuaciones no Lineales

20

3.1. Funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.2. Represent Representaci´ aci´ on on gráfica de funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

1

´ Indice general

1. Un lenguaje lenguaje de programac programaci´ i´ on on

4

1.1. Introducci´ Introducci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2. Operadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2.1. Operacionales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2. 1.2.2. 2. L´ ogicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3. Comandos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3.1. Generales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3.2. Aritméticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.3.3. Algebraicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2. Matrices, submatrices, determinantes, sistemas de ecuaciones

9

2.1. 2.1. C´ alculo matricial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.1.1. Comandos generales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.1.2. Tipos de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.1.3. Operaciones con matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2. .2. Sist Sistem emaas de ecua ecuaccione iones. s. Método todoss dir directo ectoss . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.2.1. 2.2 .1. Soluci Soluci´ón o n de sist sistem emaas Com Compati patibl bles es Deter etermi mina nado doss . . . . . . . . . . . . . . . . . . .

12

2.2.2. 2.2 .2. Soluci Soluci´ón o n de sist sistem emaas Com Compati patibl bles es Inde Indete term rmin inad ados os . . . . . . . . . . . . . . . . . .

12

2.2.3. 2.2.3. Algoritmo Algoritmo de eliminaci´ eliminaci´ on de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.2. 2.2.44. Consi onside dera raccione ioness impo imporrtan tantes tes sobr sobree sis sistema temass . . . . . . . . . . . . . . . . . . . . . .

14

3. Ecuaciones no Lineales

20

3.1. Funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.2. Represent Representaci´ aci´ on on gráfica de funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

1

Pr´ acticas acticas con Octave

2

3.3. Polinomios

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

22

3.4. Interpolaci Interpolaci´ on oń Polinómica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

3.5. Ecuaciones no lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.6. 3.6. Sist Sistem emas as de Ecua Ecuaci cion ones es no line lineal ales es.. M´ Métod e todoo de Newt Newton on.. . . . . . . . . . . . . . . . . . . .

28

4. Sistemas de ecuaciones lineales

4.1.. Soluci 4.1 Solución ó n de sist sistem emas as de ecua ecuaci cion ones es con con m´ métod e todos os iter iterat ativ ivos os . . . . . . . . . . . . . . . . . . 4.1.1. Método de Jacobi

30

30

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

30

4.1.2. Método odo de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

5. Autovalores, autovectores de una matriz

34

5.1. 5.1. C´ alculo por por método odos directos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

5.2. 5.2. C´ alculo por por método odo iterativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

6. Derivaci´ Derivaci´ on on e Integraci´ on

39

6.1. Derivaci´ Derivaci´ on on e Integración con Symbolic . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

6.2. Derivadas numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

6.3. Integraci Integraci´ón Numérica

43

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

7. Ecuaciones Diferenciales

44

7.1. Campo de direccio direcciones nes de una ecuaci´ ecuaci´ on diferencial . . . . . . . . . . . . . . . . . . . . . . .

44

7.2. 7.2. Prob Proble lema mass de Valor alor Inic Inicia ial. l. Ecua Ecuaci cion ones es Dife Difere renc ncia iale less . . . . . . . . . . . . . . . . . . . . .

46

7.2.1. Usando Symbol bolic. Con Dsolve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

7.2.2. Usando la orden Lsode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

7.2.3. Usando Métod todos numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

7.3. Sistemas de Ecuaciones Diferenciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

7.3.1. 7.3 .1. Ecuaci Ecuaci´ón de Orden Superior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

7.3.2. Usando Symbol bolic. Con Dsolve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

7.3.3. 7.3 .3. Resolu Resoluci´ ci´ on con lsode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

7.3.4. 7.3 .4. Resolu Resoluci´ ci´ on o n por por Métod e todos os Num´ uméric e ricoos . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

7.3.5. Método odos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76

7.4. Problemas de Valor Inicial y de Contor torno . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

8.

7.4.1. Problemas de transporte estacionario 1-dimensionales . . . . . . . . . . . . . . . .

80

7.4.2. Problemas de difusión evolutiva 1-dimensionales . . . . . . . . . . . . . . . . . . .

80

Otros comandos y funciones u ´ tiles

81

8.0.1. Los comandos max y min . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

8.0.2. C´ alculo de posiciones. El comando find. . . . . . . . . . . . . . . . . . . . . . . . .

81

9. Transformada de Laplace. Funciones definidas a trozos. La funci´ on Heaviside

83

9.0.1. Transformada de Laplace. Funciones definidas a trozos. La función Heaviside . . .

83

9.1. Transformada de Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

10.An´ alisis Vectorial

85

10.0.1. Representación de Campos Vectoriales . . . . . . . . . . . . . . . . . . . . . . . . .

85

10.0.2. Representación de curvas y superficies . . . . . . . . . . . . . . . . . . . . . . . . .

87

10.0.3. Cálculos de gradiente, divergencia, rotaciona. Representaciones. . . . . . . . . . . .

91

Cap´ıtulo 1

Un lenguaje de programaci´ on 1.1.

Introducci´ on

Octave es un lenguaje de programación de alto nivel especializado en cálculos numéricos. Al iniciar el programa, ya sea en un entorno gráfico o por l´ınea de comando, nos encontraremos con alg´ un mensaje de bienvenida y la l´ınea de comando en blanco para comenzar a trabajar. Existen dos maneras de trabajar con Octave: de forma directa, ingresando comandos por la l´ınea de comandos, o bien generando un programa en el editor. Un programa es un archivo de texto plano que contiene una serie de instrucciones que octave puede interpretar y ejecutar, de extensión .m.

1.2. 1.2.1.

Operadores Operacionales

1. Suma, resta, multiplicaci´ on, división >> 8 + 2 (13 ans = 31

∗

− 3) + 15/5

2. Ra´ız cuadrada >> sqrt(25) ans = 5 3. Potencia >> 3 ∧ 2 ans = 9 >> 5

∗ ∗2

ans = 25

4. Exponencial y logaritmo neperiano >> exp(log(2)) ans = 2 4

5

Pr´ acticas con Octave

5. Razones trigonométricas, sin, cos, tan, asin, acos, atan >> sin(pi) ans = 0 6. Valor absoluto >> abs(-20.14) ans = 20,14

1.2.2.

L´ ogicos

En Octave el valor verdadero viene representado por un 1 y el valor falso por un 0 1. Asignaci´ on de igualdad >> a = 3 a=3 2. Comparaci´ on de igualdad >> 7 == 3 ans = 0 >> 7 == 7 ans = 1 3. Comparaci´ on de desigualdades; menor <, menor igual <=, mayor >, mayor igual >= >> 9 >= 7 ans = 1 4. Distinto >> 5! = 5 ans = 0 5. Conjunci´ on >> a = 6; (a < 9)&(a! = 4) ans = 1 6. Disyunci´ on >> a = 5; (a > 7) (a < 6) ans = 1

|

7. Negaci´ on >> a = 5;!((a < 3) (a = 1)) ans = 1

|


1.3. 1.3.1.

Comandos Generales

1. Borrar la ventana de comandos >> clc 2. Borrar la memoria >> clear 3. Salir del programa >> exit 4. Muestra m´ as cifras decimales >> format long 5. Muestra menos cifras decimales >> format short 6. No muestra por pantalla el resultado >> 5;

1.3.2.

Aritm´ eticos

7. Mayor entero que sea menor o igual a x >> floor(3,254) ans = 3 8. Menor entero que sea mayor o igual a x >> ceil(3,254) ans = 4 9. Parte entera de x >> fix(7,254) ans = 7 10. Redondea al entero más próximo a x >> floor(3,854) ans = 4 11. Devuelve el resto de la divisi´ on >> rem(8, 5) ans = 3 12. Devuelve la descomposici´ on en factores primos >> factor(306) ans = 2 3

3 17

6

7


13. Devuelve el m´ aximo com´ un divisor >> gcd(36, 14) ans = 2 14. Devuelve el m´ınimo común m´ ultiplo >> lcm(36, 14) ans = 252 15. Calcula el factorial de un número >> factorial(6) ans = 720 16. Dice si es un n´ umero primo >> isprime(29) ans = 1 17. Devuelve los números primos menores de 11 >> primes(11) ans = 2 3

1.3.3.

5 7 11

Algebraicos

1. Vectores >> [8, 5, 1] ans = 8 5 1 2. Vectores, mediante rangos >> [8 : 12] ans = 8 9 10 11 12 >> [8 : 2 : 12] ans = 8 10 12 >> linspace(1 : 4 : 6) ans = 1 1,6 2,2 2,8 3,4 4,0 3. Matrices >> [8, 5, 1; 3, 7, 4] ans =

8 5 1 3 7 4

4. Genera una matriz aleatoria 3x2 con n´ umeros entre >> randi([ 7, 10], 3, 2) ans =

−

8 1 3

−

−5

9 6

−7 y 10


5. Genera una matriz aleatoria 3x2 con n´ umeros entre 0 y 1 >> rand(3, 2) 0,12 0,01 ans = 0,20 0,06 0,70 0,12 6. Construye una matriz cuyo polinomio caracter´ıstico tiene de coeficientes los del polinomio [1,2,3] >> compan[(1, 2, 3)]

−2 −3

ans =

1

0

7. Longitud de un vector o de una matriz, n´ umero de columnas >> length([2, 5, 1]) ans = 3 8. Dimensi´ on de un vector o de una matriz, número de filas y columnas >> size([2, 5, 1]) ans = 1

3

9. Suma los elementos de un vector >> sum([8, 5, 1]) ans = 14 10. Suma los elementos de cada columna de una matriz >> sum([8, 5; 2, 4]) ans = 10

9

11. Suma los elementos de cada fila de una matriz >> sum([8, 5; 2, 4], 2) ans =

13 6

12. Producto escalar de dos vectores >> dot([8, 5, 1], [ 1, 1, 0]) ans =

−3

−

13. Norma de un vector >> norm([ 1, 1, 0])

− √ ans = 2

8

Cap´ıtulo 2

Matrices, submatrices, determinantes, sistemas de ecuaciones 2.1. 2.1.1.

C´ alculo matricial Comandos generales

1. Extraer el elemento a12 de la matriz A >> A = [8, 5; 2, 4]; ans = 5

A(1, 2)

9

2. Extraer la segunda fila >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3];

A(2, :)

ans = 2 4 4 3. Extraer la tercera columna >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3];

A(:, 3)

1 ans = 4 3 4. Submatriz fila de 2 a 3 y columna de 2 a 4 >> A = [8, 5, 1, 7; 2, 4, 4, 3; 7, 1, 3, 1]; ans =

A(2 : 3, 2 : 4)

4 4 3 1 3 1

5. Submatriz fila 1 y 3 y columna 2 y 4 >> A = [8, 5, 1, 7; 2, 4, 4, 3; 7, 1, 3, 1]; ans =

A([1, 3], [2, 4])

5 7 1 1 9

10


6. Borra la 3a fila >> A = [8, 5, 1, 7; 2, 4, 4, 3; 7, 1, 3, 1]; ans =

A(3, :) = []

8 5 1 7 2 4 4 3

7. Unir matrices en filas >> [A, B] 8. Unir matrices en columnas >> [A; B] 9. Unir matrices por bloques >> blkdiag(A, B)

2.1.2.

Tipos de matrices

1. Matriz nula >> zeros (2,3) ó zeros(3) si es cuadrada ans =

0 0 0 0 0 0

2. Matriz de unos >> ones (2,3) ó ones(3) si es cuadrada ans =

1 1 1 1 1 1

3. Matriz unidad >> eye (3) 1 0 0 ans = 0 1 0 0 0 1 4. Matriz daigonal >> diag ([3,-1,5) 3 ans = 0 0

0 0 1 0 0 5

−

5. Matriz traspuesta >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3];A. 8 2 7 ans = 5 4 1 1 4 3 6. Matriz inversa >> A = [2, 1; 1, 1]; inv(A) ans =

−

1 1

−1

2

11


7. Matriz triangular superior >> triu(A) 2 1 0 1

ans =

8. Matriz triangular inferior >> tril(A) 2 0 1 1

ans =

2.1.3.

Operaciones con matrices

1. Multiplicar >> A B

∗

2. Potencia >> A∧ 3 3. Determinante de una matriz >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3]; det(A) ans = 148 4. Rango de una matriz >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3]; rank(A) ans = 3 5. Traza de una matriz >> trace(A) ans = 15 6. Sumar una constante a todos los elementos de la matriz >> A+10 7. Multiplica por 6 la Fila 1 >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3]; A(1, :) = 6 A(1, :) 48 30 6 ans = 2 4 4 7 1 3

∗

8. Cambia la Fila 1 y la Fila 2 >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3]; A([1, 2], :) = A([2, 1], :) 2 4 4 ans = 8 5 1 7 1 3 9. Multiplica por 6 la Fila 2 y le suma 3 veces la Fila 1, dej´ andola en la Fila 2 >> A = [8, 5, 1; 2, 4, 4; 7, 1, 3]; A(2, :) = 6 A(2, :) + 3 A(1, :) 8 5 1 ans = 36 39 27 7 1 3

∗

∗

12


2.2.

Sistemas de ecuaciones. M´ etodos directos

2.2.1.

Soluci´ on de sistemas Compatibles Determinados

Sea A =

 −

1 0 1 1 0 1

independiente.

−1 2 −1

 

la matriz de los coeficiente y B =

  −  8 0 1

, la matriz de los términos

Viendo que rank(A) = 3 resolvemos el sistema mediante Octave. >> inv(A) B

∗

13 3 5

ans =

Tambi´ en podemos unir las matrices A y B y aplicar el método de GAuss-Jordan con el comando rref (reduce row echelon from) >> Ampliada = [A, B] >> rref (Ampliada) ans =

2.2.2.

 

1 0 0

 

0 0 13 1 0 3 0 1 5

Soluci´ on de sistemas Compatibles Indeterminados

Estudiar y resolver el sistema

− −

 − −

3 4 1 1 Sea A = 1 2 términos independiente.

3x + 4y z + 2t + 4s = x y 2z t s = x + 2y 5z + 4t + 2s =

−1 2 4 − −2 −1 −1 −5 4 2

− − − − − −

 

12 6 8

− −

la matriz de los coeficiente y B =

   −−  12 6 8

>> Ampliada = [A, B] >> rref (Ampliada) Obtenemos la matriz escalonada reducida

 

1 0 0

0 1 0

−9 −7

0 0 0 1 0 1 0

−16 −8 −2

 

Las inc´ ognitas son x, y,t la z y la s son parámetros y la solución viene dada por:

, la matriz de los

13


      x y t

2.2.3.

=

9z 16 7z s 8 2

− − − −

 

∀ z, s ∈ R

Algoritmo de eliminaci´ on de Gauss

Construimos una funci´ on que calcule un sistema de ecuaciones empleando el m´ etodo de eliminación de Gauss

o bien, esta otra, empleando la matriz ampliada


2.2.4.

14

Consideraciones importantes sobre sistemas

Para calcular la solución de un sistema compatible determinado A x = b, podemos usar la orden sol=inv(A)*b que equivale a escribir A  b Pero hay que tener cuidado con el operador  ya que no siempre coincide con lo dicho. Si tenemos el sistema sobredeterminado

 

x + y = 1 2 x y = 3 , 2 x + 4 y = 1

−

Estudiando los rangos o bien con la orden rref podemos ver que es incompatible, pero si hacemos: A=[1 1;2 -1;2 4] b=[1;3;1] A  b nos da una solución ans = 1.31429 -0.38571 Que evidentemente no puede ser correcta, en efecto, A*ans nos da ans = 0.92857 3.01429 0.54286 No coincide con b, aunque son valores aproximados. ¿qué ha pasado? En primer lugar parece haber calculado la inversa de una matriz 3x2 lo cual no tiene sentido matem´ atico. Si escribimos inv(A) nos dirá que no puede calcularla porque no es ua matriz cuadrada, entonces ¿qué operación ha hecho? Realmente el operador  , cuando la matriz no es cuadrada calcula la pseudoinversa, que en realidad consiste en hacer las siguientes operaciones Amxn xnx1 = b mx1

15


Multiplicando por la traspuesta de A por la izquierda At A x = At b At A es siempre un matriz cuadrada, en nuestro caso nxn, si tiene inversa, podemos hacer x = (At A)−1 At b A la matriz (At A)−1 At , que existirá cuando (At A) tenga inversa, la llamamos pseudoinversa por la izquierda. Es una matriz que cumple que M.A = I , pero no A.M = I , ya que este u ´ ltimo producto ni siquiera tiene sentido si estudiamos las dimensiones De forma similar se puede calcular la pseudoinversa por la izquieda. La pseudoinversa de una matriz se puede calcular con el comando pinv, aunque en Octave, cuando la matriz A no es cuadrada la orden A  b calcula s=pinv(A)*b que no es la solución, pero puede ser una aproximación a la solución, ya que lo que realmente estamos haciendo es aproximar por m´ınimos cuadrados. Gr´ aficamente, en nuestro ejemplo anterior, son tres rectas que no se cortan, pero los puntos de corte de cada par de ellas está próximo, de ah´ı que consigamos una aproximación de la solución relativamente aceptable.

Si repetimos lo anterior con el sistema

 

x + y = 1 2 x y = 3 , x y = 5

− −

las soluciones serán peores aproximaciones, ya que los puntos de corte de las tres rectas están muy separados

16


A=[1 1;2 -1;1 -1];b=[1; 3 ;5] A b ans = 1.5714 -1.2857 A*ans ans = 0.28571 4.42857 2.85714

Tambi´ en es importante ver qué ocurre con sistemas compatibles indeterminados, donde aunque la matriz de los coeficientes fuese cuadrada, no tiene inversa



x + y = 1 , 2 x + 2 y = 2

es un sistema compatible indeterminado de infinitas soluciones, (α, 1

− α). Sin embargo la orden

Ab en este caso nos dará una u ńica solución ans = 0.50000 0.50000 Pero, no nos ha dado una solución cualquiera, nos ha dado la solución s que cumple que

 s  sera

17


m´ınimo. La pseudo inversa puede utilizarse para hacer una aproximación por m´ınimos cuadrados, supongamos que tenemos unos datos sobre unas tensiones ejercicidas sobre un cuerpo y el aumento en la longitud que producen



Tensiones 0 0,06 0,14 0,25 0,31 0,47 0,6 0,7 , Alargamientos 0 0,08 0,14 0,20 0,23 0,25 0,28 0,29

Buscamos si existe una relación lineal entre las variables y = x1 T + x2 siendo T las tensiones y y los alargamientos, para ello construimos las matrices

∗

A=[0 1; 0.06 1;0.14 1;0.25 1;0.31 1;0.47,1;0.6 1;0.7 1] b=[0;0.08;0.14;0.2;0.23;0.25;0.28;0.29] y planteamos el sistema A x = b que resolvemos con

∗

xs=pinv(A)*b Obtenemos xs = 0.374099 0.065441 El mismo resultado hubiesemos conseguido con la orden de interpolación polin´ omica de octave polyfit polyfit(A(:,1),recta(A(:,1),xs(1),xs(2)),1) Por tanto la recta que mejor aproxime esos puntos, recta regresión, será Y = 0,374099 X + 0,065441, la representamos tomando como puntos sobre el eje x los valores de las tensiones. figure 1 Primero dibujamos los puntos en rojo con un asterisco con la orden plot(A(:,1),b,’r*’) y en la misma gráfica, la recta, que definimos con parámetros, as´ı hacemos el programa más general para cuando quiera aplicarla a otros conjuntos de puntos recta=inline(’p1.*x+p2’,’x’,’p1’,’p2’) hold on plot(A(:,1),recta(A(:,1),xs(1),xs(2)),’b’)

18


Resultar´ a Observando la gráfica se aprecia que una aproximación por una parábola puede ser más exacta, para hacerlo tomamos una nueva matriz 8x3 donde la primera columna sean los valores de las tensiones al cuadrado. Para no hacerlo manualmente podemos teclear A1=[(A(:,1)).2 ]; B=[A1 A]; La resolvemos como antes ys=pinv(B)*b obteniendo los coficientes de la parábola ys = -0.739870 0.889101 0.018426 Hubiesemos conseguido lo mismo con la orden de interpolación polin´ omica de Octave polyfit(B(:,2), parabola(B(:,2),ys(1),ys(2),ys(3)),2) La par´ abola será Y =

−0,739870 X 2 + 0,889101 X + 0,018426

La representamos, los puntos como antes, la parábola en verde usando los puntos de las tensiones y finalmente dibujamos la parábola en azul tomando m´ as valores para que se aprecie que es una parábola figure 2 parabola=inline(’p1.*x.2 +p2.*x+p3’,’x’,’p1’,’p2’,’p3’); plot(B(:,2),b,’r*’)


hold on plot(B(:,2), parabola(B(:,2),ys(1),ys(2),ys(3)),’g’) t=[0:0.01:0.7]; yt=parabola(t,ys(1),ys(2),ys(3)); plot(t,yt,’b’) Obteniendo

19

Cap´ıtulo 3

Ecuaciones no Lineales 3.1.

Funciones

En Octave existen tres tipos de funciones: las funciones anónimas, funciones “en linea” y las funciones definidas mediante el comando function, utilizaremos ésta última para definir expresiones más complejas. En su forma más sencilla tiene la siguiente sintaxis: function nombre comandos

endfunction Donde nombre es el nombre con el que llamaremos a la función 1 . La guardaremos en un fichero con la extensión .m y con el mismo nombre que el nombre de la función.

3.2.

Representaci´ on gr´ afica de funciones

Para realizar una gráfica Octave necesita dos vectores de la misma longitud, uno x = (x1 , x2 , . . . , xn ) con las coordenadas del eje x y otro y = (y1 , y2 , . . . , yn ) con las coordenadas del eje y. De esta manera puede construir los pares (xi , yi ) y dibujar mediante interpolaci´ on de los pares (xi , yi ). Para dibujar una colecci´ on de puntos, por ejemplo tomando los puntos ( 1, 2), (1, 4), (2, 8), (3, 1) escribimos

−

x = [-1, 1, 2, 3] y = [2, 4, 8, 1] plot(x,y)

Despu´ es de cada ”pareja de vectores”(x, y) podemos a˜ nadirle atributos o detalles de color, forma, etc. Por ejemplo: 1

No utilizar tildes, espacios, la letra “˜ n” y caracteres extra˜ nos o reservados.

20

21


’b’ - blue ’r’ - red ’g’ - green ’k’ - black o’ - dibuja c´ırculos ´ ’x’ - dibuja cruces ’-’ - dibuja un recta entre dos puntos sucesivos.

Prueba como ejemplo con los siguientes comandos plot(x,y,’r’) plot(x,y,´ or’) plot(x,y,’g-’)

Para dibujar la gr´ a fica de un funci´ on en un intervalo [a,b] creamos un vector con numerosos puntos entre a y b con los comandos [a:n:b] linspace(a,b,n)

Crea un vector con entradas desda a hasta b, con pasos de n en n Crea un vector con n entradas equiespaciadas entre a y b

afica de f (x) = x 2 + x + 1 en el intervalo [ 2, 2] Ejemplo. La gr´

−

f = inline(’x.∧ 2+x+1’,’x’) I = [-2:0.01:2] plot(I,f(I))

Prueba a dibujar f (x) usando el intervalo J=[-2:1:2] ¿Qué ocurre?

Presentaci´ on de gr´ aficas Para dibujar varias funciones en la misma ventana se van a˜ nadiendo al comando plot de dos en dos: f = inline(’x.*sin(x)’,’x’); g = inline(’sin(x)./x’,’x’); I = [-5:0.01:5]; plot(I,f(I),’b’,I,g(I)),´ ok’)

Para mostrar distintas gr´ aficas en distintas ventanas se usa el comando figure antes de cada plot: figure(1) f = inline(’x.*sin(x)’,’x’); figure(2) g = inline(’sin(x)./x’,’x’);

Para mostrar varias subventanas en la misma imagen utilizamos el comando subplot(filas,cols,index):

22


f = inline(’x.*sin(x)’,’x’); g = inline(’sin(x)./x’,’x’); I = [-5:0.001:5]; subplot(2, 1, 1) plot(I,f(I)); subplot(2, 1, 2) plot(I,g(I));

Crea una matriz de 2 filas y 1 columna Dibuja f en la primera casilla Activa la segunda casilla Dibuja g en la segunda casilla

Podemos a˜ nadir después de cada plot comandos para incluir un t´ıtulo, nombres en los ejes, leyendas, etc. f1 = inline(’cos(x)’,’x’); f2 = inline(’sin(x)’,’x’); I = [-5:0.01:5]; plot(I,f1(I),I,f2(I)) title(’Funciones seno y coseno’) legend(’cos(x)’,’sin(x)’) xlabel(´ eje x’) ylabel(´ eje y’)

El comando fplot puede dibujar directamente una función conocida sin definir ningún vector, fplot (fun,limites,n´ umero de puntos) probamos fplot(’1/(1+x**2),[-10,10],10); fplot(’1/(1+x**2),[-10,10]); fplot(@cos,[-pi,pi]); fplot(’[cos(x),sin(x)]’,[0,2*pi]); fplot(’[x**3-3*x+2,1-x**2]’,[-4,4]);

o bien definiendo la función con @(x) f1=@(x) 1./x; fplot(f1,[-1,1]); ezmesh (@(x, y) x.**2 - y.**2 - 1);

3.3.

Polinomios

Un polinomio en Octave se puede representar como un vector, p(x) = x5 + 6 x3 + 2 x2 2 x + 1 se puede escribir como una función como ya hemos visto, pero cuando se trata de un polinomio tambi´ en se puede escribir como un vector utilizando los coeficientes de mayor a menor

−

p=[1 0 6 2 -2 1] De esa forma podemos hacer distintas operaciones con vectores. Para comprobar que la expresión escrita es correcta, podemos poner polyout(p,’x’) Para sumar y multiplicar por un escalar, si q=[3 2 1 9]

23


p+2*q dar´ıa un error, ya que los vectores no son de la misma dimensión, luego habra que escribir q=[0 0 3 2 1 9] Para multiplicar prod=conv(p,q) Ser´ıa conveniente quitar los ceros que pueden aparecer a la derecha, con la orden reduce(prod) Para dividir, el comando deconv nos indica el cociente y el resto de la división, para ello q no puede empezar por 0, luego antes hemos de quitarlos con polyreduce q=polyreduce(q) [q,r]=deconv(p,q) Para obtener una descomposición en fracciones simples usamos la orden residue, con la forma [coef,den,coc,exp]=residue(p,s) donde coc son los exponentes del polinomio cociente, y la expresión del resto en fracciones simples nos la dan los otros datos, coef son lo coeficientes de la expresión racional, den dará unos valores x k que indican que los denominadores son x xk , elevados a los exponentes que indique exp, todos ordenados respectivamente. Por ejemplo si tomamos

−

s=[1,-2,1,0] [coef,den,coc,exp]=residue(p,s) nos dará coef = 17 8 1 den = 1 1 0 k= 129 exp = 1 2

24


1 luego podemos decir que x2 + 2 x + 9 +

17 x

−1

+

8 (x

2

− 1)

+

1 . x

Para dar valores a un polinomio la orden es polyval t=linspace(-2,2,20); p1=polyval(p,t) Podemos representar el polinomio con plot(t,p1) Y calcular las raices del polinomio con roots(s) roots(p) Aunque hay que precaución con esa orden, puede inducir a errores cuando el polinomio tiene ra´ıces multiples, para un polinomio con 1 raiz de multiplicidad 8, no dar´ıa p8=[1 -8 28 -56 70 -56 28 -8 1] roots(p8) ans = 1.02171 + 0.00000i 1.01523 + 0.01540i 1.01523 - 0.01540i 0.99983 + 0.02154i 0.99983 - 0.02154i 0.98477 + 0.01506i 0.98477 - 0.01506i 0.97863 + 0.00000i Tambi´ en permite calcular derivadas e integrales indefinidas de polinomios con las ordenes polyder(p) polyint(p) El resultado no dará el término independiente 0, pero si se cambia por C se tiene la integral indefinida.


3.4.

25

Interpolaci´ on Polin´ omica

La orden polyfit(x,y,n) permite calcular un polinomio de grado n que pase por los puntos determinados por los vectores x e y. Si x=[ 1 2 3 4 5 6 7 8 9 10 11 12] e y=[ 2.7 5.1 7.2 9 10.7 12 13 13.7 14 14.2 14.1 13.7] Como hay 12 valores, probamos con un polinomio de grado 11 p11=polyfit(x,y,11) x11=linspace(1,12,30) y11=polyval(p11,x11) Y las comparamos gráficamente plot(x,y, o ,x11,p11, r )

¿c´ omo funciona realmente esta interpolación? Para conseguir un polinomio de grado n que aproxime un conjunto de puntos, o una función f (x), buscamos una polinomio que pase por los puntos (xi , yi ) para i = 1,..n, en el caso de una función yi = f (xi ).

26

Pr´ acticas con Octave n

Vamos a utilizar la función auxiliar φ k (x) =



j = 0, j <> k .

x xk

−x −x

j

, que vale 1 cuando x = x k y 0 en los

j

demás x i , i <> k. n

De esta forma el polinomio de interpolación será p(x) =



φk (x).

k

Por ejemplo, el polinomio que pase por los puntos (1 , 2), ( 1, 3) y (3, 5), será de segundo grado y lo construiremos

−

(x + 1)(x p(x) = (1 + 1)(1

3.5.

− 3) 2 + (x − 1)(x − 3) 3 + (x − 1)(x + 1) 5 = x2 − x + 4 2 − 3) (−1 − 2)(−1 − 3) (3 − 1)(3 + 1)

Ecuaciones no lineales

Si queremos aproximar la ra´ız de la ecuación f (x) = 0 con una tolerancia tol y en un n´ umero máximo de nmax de iteraciones, se define la función f como una funci´ on inline con el comando f = inline(’...’,’x’)

y se elige uno de los siguientes métodos:

M´ etodo de la bisecci´ on Se eligen a y b de manera que f contenga una u ńica ra´ız en el intervalo [a, b] y se ejecuta el comando: [z,res,niter] = bisection(f,a,b,tol,nmax)

En la salida z contiene la aproximación buscada, res el residuo y niter el n´ umero de iteraciones que ha tardado el método en alcanzar z .

M´ etodo de Newton. Se define la derivada de f como otra función inline df = inline(’...’,’x’)

Se elige el valor inicial x 0 y se ejecuta el comando [z,res,niter] = newton(f,df,x0,tol,nmax)

En la salida z contiene la aproximación buscada, res el residuo y niter el n´ umero de iteraciones que ha tardado el método en alcanzar z .


27

M´ etodo del punto fijo y m´ etodo de Aitken. Si queremos resolver el problema de punto fijo φ(x) = x, primeramente se define φ como una función inline con el comando phi = inline(’...’,’x’)

Se elige el valor inicial x 0 y se ejecuta el comando [z,niter] = puntofijo2(phi,x0,tol,nmax)

En la salida z contiene la aproximación buscada, res el residuo y niter el n´ umero de iteraciones que ha tardado el método en alcanzar z . Si queremos aplicar el método de aceleraci´ on de Aitken para resolver el mismo esquema se ejecuta el comando [z,niter] = aitken(phi,x0,tol,nmax)

Orden fzero. fzero nos da una aproximación a partir de un punto por métodos num´ ericos, pero no sabemos qué método est´ a utilizando, y por tanto no lo podemos controlar fun=@(x) x. 3

− 3x + 1

fplot(fun,[-3,3]) fzero(fun,1) ans = 0.34730 fzero(fun,-1) ans = -1.8794 fzero(fun,0) ans = 0.34730 fzero(fun,1) ans = 0.34730 fzero(fun,3) ans = 1.5321

28


3.6.

Sistemas de Ecuaciones no lineales. M´ etodo de Newton.

Sea el sistema de ecuaciones no lineales

  

f 1 (x1 , x2 , . . . , xn ) = 0 f 2 (x1 , x2 , . . . , xn ) = 0 .. . f n (x1 , x2 , . . . , xn ) = 0

que podemos escribir de manera compacta como f (x) = 0

en donde f = (f 1 , . . . , fn )T y x = (x1 , . . . , xn )T . (0) T Tomando como dato inicial el vector columna x0 = (x(0) Rn , para resolver el sistema 1 , . . . , xn ) mediante el método de Newton necesitaremos crear 3 archivos .m distintos:

∈

1. Archivo Ffun.m que contiene al vector columna f = (f 1 , . . . , fn )T . 2. Archivo JFun.m que contene la matriz Jacobiana J f . 3. Archivo en donde ejecutaremos el método mediante el programa newtonsys.m. Debemos tener en cuenta además las siguientes observaciones: La llamada a las funciones Ffun y Jfun re realiza usando el s´ımboo @. Dentro de los archivos Ffun.m y Jfun.m las variables x 0 , . . . , xn se denotan por x(0), . . . , x(n). No es necesario escribir punto antes de ,

∗

∧

y /.

Ejemplo: Para el sistema de ecuaciones no lineales Sea el sistema de ecuaciones no lineales



x21 10x1 + x22 + 8 = 0 x1 x22 + x1 10x2 + 8 = 0

−

−

y tomando como semilla x 0 = (0,5, 0,5)T tenemos que f =



x21 10x1 + x22 + 8 x1 x22 + x1 10x2 + 8

−

−



y J f (x) =



−

por lo que creamos los archivos Ffun.m y JFun.m de la siguiente forma: Archivo Ffun.m function F=Ffun(x) F(1,1) = x(1)∧ 2-10 x(1)+x(2)∧ 2+8; F(2,1) = x(1) x(2)∧ 2+x(1)-10 x(2)+8; endfunction

∗

∗

∗



2x1 10 2x2 2 x2 + 1 2x1 x2 10

−


Archivo Jfun.m function J=Jfun(x) J(1,1) = 2 x(1)-10; J(1,2) = 2 x(2); J(2,1) = x(2)∧ 2+1; J(2,2) = 2 x(1) x(2)-10; endfunction

∗ ∗ ∗

∗

Ahora el sistema se resuelve mediante los comandos: x0 = [0.5;0.5]; [z,res,niter] = newtonsys(@Ffun,@Jfun,x0,tol,nmax)

29

Cap´ıtulo 4

Sistemas de ecuaciones lineales En el cap´ıtulo anterior hemos utilizado funciones para aplicar el método de eliminaci´ on de Gauss.

4.1.

Soluci´ on de sistemas de ecuaciones con m´ etodos iterativos

Un método iterativo es un método que progresivamente va calculando aproximaciones a la soluci´ on de un problema. En Matemáticas, en un m´ etodo iterativo se repite un mismo proceso de mejora sobre una solución aproximada: se espera que lo obtenido sea una solución m´ as aproximada que la inicial. El proceso se repite sobre esta nueva solución hasta que el resultado m´ as reciente satisfaga ciertos requisitos. A diferencia de los métodos directos, en los cuales se debe terminar el proceso para tener la respuesta, en los m´ etodos iterativos se puede suspender el proceso al termino de una iteraci´ on y se obtiene una aproximaci´ on a la solución. Un método iterativo consta de los siguientes pasos. 1. Inicia con una soluci´ on aproximada (Semilla) 2. Ejecuta una serie de cálculos para obtener o construir una mejor aproximación partiendo de la aproximaci´ on semilla. La f´ ormula que permite construir la aproximación usando otra se conoce como ecuación de recurrencia. 3. Se repite el paso anterior pero usando como semilla la aproximación obtenida.

4.1.1.

M´ etodo de Jacobi

El m´ etodo Jacobi es el método iterativo para resolver sistemas de ecuaciones lineales más simple y se aplica solo a sistemas cuadrados, es decir a sistemas con tantas incógnitas como ecuaciones. CX = B

(4.1)

1. Primero se determina la ecuación de recurrencia. Para ello se ordenan las ecuaciones y las incógnitas. De la ecuación (4.1) se despeja la incógnita X . En notación matricial se escribirse como: X = B + AX 30

31


2. Se toma una aproximaci´ on para las soluciones y a ésta se le designa por X 0 3. Se itera en el ciclo que cambia la aproximación X n+1 = B + AX n Uno de los principales problemas de los métodos iterativos es la garant´ıa de que el método va a converger, es decir, va a producir una sucesión de aproximaciones cada vez efectivamente más próximas a la solución. En el caso del m´ etodo de Jacobi no existe una condici´ on exacta para la convergencia. Lo mejor es una condición que garantiza la convergencia, pero en caso de no cumplirse puede o no haberla es la siguiente: Si la matriz de coeficientes original del sistema de ecuaciones es diagonalmente dominante, el método de Jacobi seguro converge. Una matriz se dice matriz diagonalmente dominante, si en cada uno de los renglones, el valor absoluto del elemento de la diagonal principal es mayor que la suma de los valores absolutos de los elementos restantes del mismo renglón. A veces la matriz de un sistema de ecuaciones no es diagonalmente dominante pero cuando se cambian el orden de las ecuaciones y las incógnitas el nuevo sistema puede tener matriz de coeficientes diagonalmente dominante. Con octave construimos una función itermeth que realiza el método iterativo general. functuion [x, iter] = itermeth(A,b,x0 ,nmax,tol,P ) Donde A, B son las matrices de los coeficientes y la matriz de los términos independientes, x0 es la condición inicial, nmax el n´ umero máximo de iteraciones, tol la tolerancia y P la opción de elegir el método de Jacobi o del Gauss-Seidel.

4.1.2.

M´ etodo de Gauss-Seidel

El método de Gauss-Seidel es muy semejante al método de Jacobi. Mientras que en el de Jacobi se utiliza el valor de las incógnitas para determinar una nueva aproximación, en el de Gauss-Seidel se va utilizando los valores de las incógnitas recien calculados en la misma iteración, y no en la siguiente. Por ejemplo, en el método de Jacobi se obtiene en el primer cálculo X i+1 , pero este valor de X no se utiliza sino hasta la siguiente iteración. En el m´ etodo de Gauss-Seidel en lugar de eso se utiliza de X i+1 en lugar de X i en forma inmediata para calcular el valor de Y i+1 de igual manera procede con las siguientes variables; siempre se utilizan las variables recien calculadas.


32

Con esta función obtenemos algunos ejemplos para obtener la solución de un sistema por el método de Jacobi.

Obteniendo en la iteración veinticuatro una aproximación de acuerdo con la tolerancia


El mismo ejemplo lo hacemos por Gauss-Seidel

Obteniendo en la iteración trece una aproximación de acuerdo con la tolerancia

33

Cap´ıtulo 5

Autovalores, autovectores de una matriz 5.1.

C´ alculo por m´ etodos directos

Veamos los comamdos para obtener los valores que diagonalizan una matriz. 1. Autovalores >> A = [ 2, 1, 2; 1, 0, 2; 1, 1, 1]; lambda = eig(A)

−

lambda =

1 1 1

−

−

− −

2. Coeficiente del polinomio caracter´ıstico >> A = [ 2, 1, 2; 1, 0, 2; 1, 1, 1]; poly(A)

−

ans = 1 1

−

−1 −1

−

3. Otra forma de obtener los autovalores >> A = [ 2, 1, 2; 1, 0, 2; 1, 1, 1]; roots(poly(A))

−

ans =

− −

1 + 0i 1 + 0i 1 + 0i

−

−

4. Salida de matriz con autovectores (matriz de paso) y matriz diagonal de autovalores (matriz diagonal) >> A = [ 2, 1, 2; 1, 0, 2; 1, 1, 1]; [vectorpropio, valorpropio]=eig(A)

−

−

vectorpropio = 0,904534 0,577350 0,301511 0,577350 0,301511 0,577350

− − −

−

−0,636162 −0,768704 0,066271

valorpropio = Diagonal Matrix 34

35


−1,00000

0 0 1,00000 0 0

0 0 1,00000

−

5. Un programa para obtener la diagonalizaci´ on de una matriz

5.2.

C´ alculo por m´ etodo iterativo

Existen métodos num´ ericos que nos conducen, a partir de la matriz del endomorfismo, a los autovectores y autovalores con una cierta aproximación. Uno de estos métodos es el de las potencias, aplicable cuando la matriz es de diagonal estrictamente dominante (los elementos de la diagonal principal son mayores en valor absoluto, que la suma de los valores absolutos de los demás elementos de la fila correspondiente) Sean λ1 , λ2 , λn los autovalores (iguales o no) de la matriz Anxm . Se dice que λi es dominante si λi > λj i = i. Sea xi su autovector correspondiente a λ i , se dice que xi es autovector dominante.

··· | | | | ∀

M´ etodo de las potencias

Sea A nxm una matriz diagonalizable =

⇒ ∃λ /Au = λ u i

i

i i

Sea λ 1 dominante Sea B= ui una base de autovectores asociados

{ }

n

Sea x0

∈ V − {0}

x0 =

n



ci ui tomado al azar

1

Definimos la ecuación recurrente xn+1 = Axn y obtenemos;

k

k

xk = A x0 = λ 1

Y por ser l´ım

k→∞

c1 u1 + c2

k

  λ2 λ1

u2 +

k

  λj λ1



= 0,

obtenemos xk

k

≈ λ1 c1u1

······ + c

n

   λn λ1

k

un

36


Con este procedimiento xk se aproxima al autovector u1 . Para que xk no crezca rápidamente, se normaliza dividiéndolo por la mayor de sus coordenadas en valor absoluto y as´ı hacemos que el autovector u1 su coordenada mayor sea uno. El autovalor dominante asociado lo obtenemos resolviendo la ecuación: A u1 = λ 1 u1 Si el autovector dominante tiene grado de multiplicidad r se converger´ıa a uno de los autovectores. Puede encontrarse los restantes autovectores, repitiendo el proceso con un vector inicial x0 distinto. Para calcular el autovalor de menor valor absoluto aplicamos el método de las potencias a A −1 , pues sabemos, que los autovalores de A−1 son los inversos de A. Para obtener los restantes autovectores y autovalores aplicamos los procedimientos de deflaci´ on. Con Octave hemos construidos el siguiente método iterativo mediante la funci´ on eigpower

Aqu´ı ponemos cuatro ejemplos de aplicaci´ on del método de las potencias para calcular aproximadamente los autovalores y autovectores.


Ejenplo1

Ejemplo2

37


Ejemplo 3

Ejemplo4

38

Cap´ıtulo 6

Derivaci´ on e Integraci´ on 6.1.

Derivaci´ on e Integraci´ on con Symbolic

La orden para derivar con Octave Symbolic es diff, para integrar int. El polinomio de Taylos se puede obtener con la orden taylor(f,variable,punto,órder’,orden). Se pueden realizar tanto de una como de más variables >> pkg load symbolic >> syms x y >> f = x 2 + 2 x exp(x) + 1

∗ ∗

>> subs(df,x,4) Se obtiene el valor de la función en el x = 4 >> df=diff(f) Calcula la derivada >> subs(df,x,4) Se obtiene el valor de la derivada en x = 4 >> db=diff(f,4) Calcula la derivada de cuarto orden >> in=int(f) Calcula la integral indefinida >> in=int(f,[2,5]) Calcula la integral definida entre 2 y 5 >> tayf=taylor(f, x ,2, o´rder’, 6) Calcula el polinomio de Taylor en el punto 2 de orden 6. Para funciones de dos variables

39

40


>> g = x y + x5

∗

− y ∗ sin(x)

Puede calcularse derivalas de distintos órdenes en ambas variables >> diff(g,x) >> diff(g,x,2,y,1) E integrales >> int(g,x) >> int(int(g,y),x) >> int(int(g,x),y) >> int(int(g,y,0,2),x,0,pi) >> int(int(g,x,0,y),y,0,pi/2) O el polinomio de Taylor de la función de 2 variables >> taylor(g, [x,y] , [1,1], órder’, 4)

6.2.

Derivadas num´ ericas

Dada una función continua y derivable, buscamos una aproximación de la primera derivada en un punto de un intervalo [a,b]. Utilizando el desarrollo en serie de Taylor, de orden uno, en torno a x 0 y si h es una cantidad positiva suficiente peque˜ na, buscamos una aproximación de la derivada en x 0 + h f (x0 + h) = f (x0 ) + h.f  (x0 ) + Despejando f  (x0 ) tenemos

h2  f (ξ 1 ) 2

x0 < ξ 1 < x0 + h

(6.1)

f (x0 + h) f (x0 ) Θ(h) h donde Θ(h) es el error de truncamiento, en ese caso de orden 1. De manera que una primera aproximación de la derivada la obtenemos con f (x0 + h) f (x0 ) f  (x0 ) , h que llamamos diferencia finita progresiva. Si repetimos el proceso en x 0 h obtenemos

−

f  (x0 ) =

−

−

≈

−

f (x0

2

− h) = f (x0) − h.f (x0) + h2 

f  (ξ 2 )

x0

− h < ξ 2 < x0

(6.2)

y llegamos a f  (x0 )

0 − h) ≈ f (x0) − f (x h

que llamamos diferencia finita regresiva, otra posible aproximación de la derivada. Restando las ecuaciones 6.1 y 6.2 obtenemos f (x0 + h) f (x0 h) = 2h.f  (x0 ) + Θ(h2 ) (6.3)

−

−

41


f  (x0 )

− f (x0 − h) ≈ f (x0 + h) 2h

que denominamos diferencia finita centrada, y donde el error es menor, ya que es de orden 2. Con el programa siguiente definimos una función que calcula la derivada utilizando la ecuació n de derivada finita centrada, salvo en los extremos, donde utiliza la progresiva y regresiva, respectivamente. Dado un vector x, y calculamos y = f (x), esa función nos dará un vector d, donde cada componente es una aproximación de la derivada en el punto correspondiente del vector x de la función f.

Si tomamos >> x=[0:1:10] >> y = x. 2 >> d=der2p(x,y) Obtendremos ans = 1 2 4 6 8 10 12 14 16 18 19 Hay errores sólo en los extremos, ya que se ha tomado h muy grande Si hubieses tomado >> x=[2:0.01:2] el error en los extremos hubiese sido también mayor, pero de 2, siendo exactos el resto de valores. También se pueden obtener fórmulas con más puntos. Muy usadas son las de tres y cinco puntos. Para deducir la de tres puntos, se parte del desarrollo de Taylor de segundo grado en f (x0 + h) y f (x0 + 2h),

42


y se obtendr´ıa las fórmulas progresivas, regresivas y centradas. Ser´ıan respectivamente + h) − f (x0 + 2h) , ≈ −3f (x0) + 4f (x02h f (x0 − 2h) − 4f (x0 − h) + 3f (x0 ) f (x0 ) ≈ ,

f  (x0 ) 

f  (x0 )

≈ f (x0 − 2h) − 8f (x0 −

2h h) + 8f (x0 + h) + f (x0 + 2h) , 12h

la centrada de orden 4 y las otras dos de orden 2. También se pueden deducir para derivadas de orden mayor.


6.3.

Integraci´ on Num´ erica

43

Cap´ıtulo 7

Ecuaciones Diferenciales 7.1.

Campo de direcciones de una ecuaci´ on diferencial 

x

Vamos a representar un campo de direcciones, por ejemplo de la ecuación diferencial y =ex −y . Para ello, definimos la función, el n´ umero de puntos que vamos a tomar y el dominio dónde queremos dibujar el campo. >> clear >> f = inline( (exp(x)

 

 



− y)./x , x , y )

>> N = 15; >> xmin =

−20;

>> xmax = 20; >> ymin =

−40;

>> ymax = 30; Ahora con la orden meshgrid creamos la red de puntos en el plano >> [x, y] = meshgrid(xmin : (xmax

− xmin)/N : xmax, ymin : (ymax − ymin)/N : ymax);

Calculamos un vector en cada uno de esos puntos, partimos de un valor 1 para x, calculamos y y luego normalizamos >> dx = ones(size(x)); >> dy = f (x, y); >> L = sqrt(dx.2 + dy. 2 ); Para evitar el error que pueda producir el caso L = 0 hacemos: >> L = L + 1e

− 10;

>> dx = dx./L;

44

45


>> dy = dy./L; Y con la orden quiver representamos esos vectores unitarios en cada uno de los puntos >> quiver(x,y,dx,dy); Y representamos >> title(’Campo de direcciones’); >> xlabel( x ); >> ylabel( y );

Si representamos del mismo modo el campo de direcciones de y =

−y2 + x resulta

46


7.2.

Problemas de Valor Inicial. Ecuaciones Diferenciales

Para resolver un Problema de Valor Inicial (PVI) de la forma



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

t

∈ [t0, t

N ]

Vamos a resolverlo desde tres puntos de vista, usando Octave Symbolic, con la orden Lsode y por métodos numéricos programados.

7.2.1.

Usando Symbolic. Con Dsolve

La orden dsolve del paquete de Octave Symbolic resuelve algunas ecuaciones sencillas. Puede utilizarse instalando el programa o bien de manera Online. Ejemplo1. Resolvamos una ecuación conocida y’+2y=-cos(t) Con el paquete simb´ olico cargando definimos las variables con las dependencias con syms >> pkg load symbolic >> syms y(x) Escribimos la ecuación y aplicamos la orden dsolve, obtenemos la solución que simplificada con simplify podemos comprobar que es la misma que calcular´ıamos sin usar ordenador, como una lineal que es. >> de3 = diff (y) + 2 y ==

∗

−cos(t);

>> sg = dsolve(de3)

Con condiciones iniciales, por ejemplo y(0)=2,

47


Para poder representar la curva, podemos convertir la solución en función con la orden function handle, proceder´ıamos:

Obteniendo la gráfica

√

Ejemplo 2. Resolver y”’+ 2y”+y’=0 con y(0)=y’(0)=0 y”(0)=1. Hallar la ecuación general, la particular y representar. Veremos en este ejercicio que a pesar de usar Octave Symbolic, nos dará valores aproximados. Si lo resolviésemos anal´ıticamente la ecuación, las soluciones de la ecuación caracter´ıstica son 0 y 2/2(1 i), con lo cual los valores que da Octave no son exactos.

−√

±

48


7.2.2.

Usando la orden Lsode

Octave dispone de un driver para resolver problemas de Cauchy llamado lsode. Las opciones de integraci´ on se modifican utilizando la función lsode options. Octave resuelve ecuaciones diferenciales no lineales, para ello utiliza la función lsode (the Livermore Solver for Ordinary Differential Equations). Planteamos el Ptoblema de Cauchy (P.C.)



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

t

∈ [t0, t

N ]

Cabe mencionar que la solución que proporciona Octave es puramente num´ erica, es una matriz x, con cada fila correspondiente a un elemento del vector t. El primer elemento de t debe ser t 0 (t sub´ındice 0 o t inicial) y debe corresponder al estado inicial del sistema x 0 , de modo que la primera fila de la salida es x0 y utilizando dicha matriz puedes obtener con plot una gráfica demostrativa de cómo se comporta la ecuación en un intervalo determinado. As´ı que no esperemos resultados con variables, diferenciales y/o

49


constantes como te lo dar´ıa un método convencional de solución. Hay 4 pasos para resolver una ecuación diferencial en octave con lsode: 1. El primer paso es definir una función que represente a la ecuación diferencial. 2. El segundo paso es definir un vector que represente los valores que va a tomar la variable independiente (t). 3. El tercer paso es definir la condición inicial para la ecuación. 4. El cuarto y u ´ ltimo paso es usar el comando lsode para integrar y resolver la ecuación diferencial, que debemos representar para estudiar la solución Veamos un ejemplo. Resolver x 2 x = e −t con x(0)=1 en el intervalo [0,20]. Defino la función, se puede definir de distintas formas, utilizo la función an´ onima >> f 0 = @(x, t) exp( t)/x2 ;

−

>> x0=1; >> t=linspace(0,20); >> y=lsode(f0,x0,t); >> plot(t,y)

Ley de Enfriamiento de Newton. Cuando la diferencia de temperaturas entre un cuerpo y su medio ambiente no es demasiado grande, el calor transferido en la unidad de tiempo hacia el cuerpo o desde el cuerpo por conducción, convección y radiación es aproximadamente proporcional a la diferencia de temperatura entre el cuerpo y el medio externo. De manera matemática, si T(t) es la temperatura en el instante t, podemos expresarlo como dT = k (T dt

− T a)

50


donde k es una constante, que llamamos constante de enfriamiento, que en realidad es el producto del coeficiente de intercambio de calor, que depende del medio, por la superficie de contacto del cuerpo. Ejemplo. Una placa de acero se extrae de un horno a 600 0 C y se sumerge en un baño de aceite a 300 C. Se sabe que la constante de enfriamiento de la pieza es 0 ,023. Determinar la temperatura que tendrá después de 92 segundos. Tenemos que resolver el P.C.



T  (t) = 0,023(30 y(0) = 600

− T ),

t

∈ [0, 92]

Lo hacemos con la orden Isode. >> fun = @(y, t) 0,023 (30

∗

− y)

>>t=[0:0.1:92] >>ys=lsode(fun,600,t) >>ys(plot(t, ys)

Se puede ver que la solución es cercana a 100, pero podemos precisar más, como hemos dividido el intervalo [0,92] con un paso de 0.1, el último valor, correspondiente a 92, se encontrará en la posición 921. >> ys(921) ans = 98.692 Ademas de Lsode Hay otras órdenes para tratar ecuaciones diferenciales de forma numérica en un paquete de Octave llamado odepkg, pero no las vamos a tratar. Si resolvemos este problema manualmente podemos ver que el resultado es exacto. Si la orden la escribimos

51


>>[ys,i,m]=lsode(fun,600,t) Si el procedimiento es satisfactorio i nos dará 2, y cualquier otro valor en caso contrario; m nos dirá si es satisfactorio o nos proporcionará información cuando no lo sea.

7.2.3.

Usando M´ etodos num´ ericos

Sea I un intervalo de la recta real, no reducido a un único punto, y sea t0 un punto fijado de él. Supongamos adem´ as que f es una función definida y continua en IxR y con valores en R, y sea y0 un valor dado de R. En estas condiciones nos planteamos el siguiente problema: Hallar una función continua y diferenciable y(t) definida en I y con valores en R verificando:



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

t

∈ [t0, t

N ]

Este tipo de problemas se conoce con el nombre de problema de Cauchy para la ecuación diferencial. La condici´ on asociada se denomina condición de Cauchy. Toda función y(t) verificando la EDO del problema de Cauchy se denomina integral de la EDO o solución particular de la EDO. En general las soluciones de la EDO estarán dadas por la expresión formal: y(t) = C +



f (t, y(t))dt

Donde C es una constante arbitrariamente elegida. A esta familia de soluciones se la denomina solución general de la EDO. En algunos casos puede haber además otras soluciones que no pertenezcan a la familia de la solución general. Dichas soluciones se denominarán soluciones singulares de la EDO. A aquella solución particular que, además de satisfacer la EDO del problema de Cauchy, verifique la condici´ on de Cauchy se la denomina solución del problema de Cauchy. Teorema (de Cauchy-Lipschitz). Suponiendo que la función f(t,y) es una función continua sobre IxR que es lipschitciana respecto a la segunda de sus variables, es decir:

∃k > 0 / |f (t, y) − f (t, z)| ≤ k ||y − z|| ∀tεI, ∀y,zεR entonces el problema admite una única solución. En numerosos problemas f´ısicos la variable t representa el tiempo. Adem´ as el intervalo I se suele considerar de la forma [t0, t0 + T] y por ello al instante t0 se le llama instante inicial y a la condición impuesta en este instante se la denomina condición inicial. En esta situación a los problemas de Cauchy correspondientes se les denomina habitualmente problemas de valor inicial (P.V.I.), extendiéndose esta terminolog´ıa también a los problemas en los que t no se identifique con la variable f´ısica tiempo. Los problemas de Cauchy que no sean problemas de valor inicial por no ser t0 el extremo izquierdo del intervalo I pueden reducirse a problemas de valor inicial mediante sencillos cambios de variable. En efecto si I= [a,b] y t0 es un punto intermedio de I=[a,b] el problema de Cauchy correspondiente puede descomponerse en dos problemas de la forma:

P C 1



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

t

∈ [a, t0]

52


P C 2



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

t

∈ [t0, b]

El problema PC2 es un problema de valor inicial. Y el problema PC1 puede transformarse en un problema de valor inicial si se realiza el cambio de variable independiente: t =-t. Por este motivo los métodos numéricos que se estudiarán en este tema se plantearán sobre problemas de valor inicial. Para ellos, además, teniendo en cuenta la expresión de la solución general de la EDO, se tendrá que la expresión formal de la solución del problema de Cauchy está dada por: t

y(t) = y 0 +



f (x, y(x))dx.

t0

Aunque existen métodos espec´ıficos para abordar este tipo de problemas (como por ejemplo los llamados métodos de tiro), por cuestiones de tiempo disponible no podremos abordarlos. Su tratamiento numérico, no obstante, puede realizarse mediante métodos en diferencias finitas. Existen numerosos tipos de métodos para la resoluci´ on de problemas de valor inicial. Introduciremos los métodos en diferencias finitas. Este tipo de métodos no persiguen la determinación de la expresión de la función solución y(t) de un problema de valor inicial. Tan sólo persiguen calcular, de forma aproximada, el valor de la solución en ciertos puntos t i ( i = 0, 1, 2, ..., N) seleccionados (bien de antemano por el usuario de dichos métodos, bien de forma autom´ atica a trav´ es de algoritmos preparados para su cálculo con el objeto de asegurar ciertas tolerancias de error entre el valor exacto (desconocido) y el valor aproximado por el método). Para su descripción, siendo I = [t0 , t 0 +T] consideremos dados los puntos en los que se va a determinar la solución y que denotaremos por: t0 < t1 < t2 < .... < t n < tn + 1 < .... < tN = t 0 + T Los m´ etodos en diferencias finitas, en general, se basan en, conocidos los valores aproximados de la soluci´ on en los puntos t n−k , tn−k+1 , tn−k+2 ,.....,tn−1 , que denotaremos como: yn−k

≈ y(t

n−k ), yn−k +1

≈ y(t

n−k +1 ),.....,yn−1

≈ y(t

n−1 )

interpolar, habitualmente mediante un polinomio, la función y(t) sobre el soporte tn−k , ....., tn . Esta técnica permite escribir yn = Φ(yn−k , yn−k+1 ,.....,yn ) En este proceso hay muchas opciones posibles que nos conducir´ıan a expresiones diferentes de la función Φ. Pero en todo caso ser´ıan expresiones en las que el valor aproximado yn se estimará en función de los valores aproximados, previamente calculados, yn−k , yn−k+1 ,.....,yn−1 , donde puede aparecer eventualmente el propio valor yn que se desea aproximar. Por ello, en general, se entenderá por método numérico para resolver un problema de valor inicial todo esquema de cálculo (algoritmo) que nos permita evaluar yn mediante una expresión de la forma: yn = Φ(yn−k , yn−k+1 ,.....,yn−1 )(Expl´ıcito) yn = Φ(yn−k , yn−k+1 ,.....,yn ) (Impl´ıcito) En el primer caso diremos que es un método expl´ıcito pues el cálculo de yn tan sólo implica evaluar la función en los valores anteriores. En el segundo caso, se dirá impl´ıcito pues el valor de yn depende del propio yn , y para determinarlo exigirá resolver la ecuación algebraica. En general son métodos más costosos en cada uno de los pasos de integración (aunque puedan presentar también ventajas como el permitir considerar pasos de integración de mayor tama˜ no, como más adelante detallaremos). A analizar el comportamiento de este tipo de métodos num´ ericos dedicaremos algunos ejemplos posteriormente.

53


Adem´ as tales tipos de métodos diremos que son métodos de k pasos en el sentido de que para evaluar ´ ltimos valores previamente calculados. En el caso de que k = 1 los métodos corresyn utilizan los k u pondientes se dicen métodos de pasos libres o m´ etodos unipaso. Por el contrario si k > 1 los métodos correspondientes se dicen métodos de pasos ligados o métodos multipaso. Los métodos multipaso k pasos sólo podrán utilizarse a partir del conocimiento de los k primeros valores, y para la determinación de estos deberán emplearse métodos de un menor n´ umero de pasos. La combinación de unos y otros no es trivial pues interesa no estropear la precisión que pueda conseguirse con el método de k pasos por utilizar métodos menos precisos en las etapas iniciales. Además, conviene advertir ya, desde el comienzo, que la consideraci´ on de un n´ umero mayor de pasos no implica necesariamente la mejor´ıa de las aproximaciones obtenidas. Es relativamente frecuente que un método de un paso bien elegido conduzca a soluciones más baratas de obtener e igual o más precisas que un método multipaso. M´ eto do de Euler expl´ ıcito

A partir del P.V.I.



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

t

∈ [t0, t0 + T ]

y una subdivisi´ on del intervalo [t0, t0 + T] mediante los puntos: t0 < t1 < t2 < .... < tn < tn+1 < .... < tN = t 0 + T Designaremos por hn a los valores (tn+1 tn ) y por y n a las aproximaciones que se vayan obteniendo de y(tn ), (n = 0, 1, ..., N). Si se supone que la función y(t) es al menos de clase C 2 ([t0 , tN ]), puede expresarse el valor de y(t) en el punto tn+1 en función del valor de y(t) en tn mediante el siguiente desarrollo en serie de Taylor: h 2n   y(tn+1 ) = y(tn + hn ) = y(tn ) + hn .y (tn ) + y (tn ) + ... 2

−

⇒

h2n  y(tn+1 ) = y(tn + hn ) = y(tn ) + hn .f (tn , y(tn )) + y (tn ) + ... 2 Si hn es suficientemente peque˜ no, el desarrollo anterior podrá aproximarse despreciando los términos en los que aparezcan potencias de hn superiores o iguales a 2, obteniéndose as´ı nuevamente el esquema de cálculo del método de Euler: yn+1 = y n + hn .f (tn , y(tn )) Este esquema de cálculo de las soluciones aproximadas del Problema de Valor Inicial (P.V.I.) se conoce con el nombre método de Euler expl´ıcito en honor al matemático Leonard Euler. M´ eto do de Euler impl´ ıcito

Si en lugar de expresar el valor de y (tn+1 ) en funci´ on del valor de la función y(t) y de sus derivadas en t n , se hiciera al rev´ es y se expresara el valor en tn en función del valor de y(t) y de sus derivadas en tn+1 se obtendr´ıa: y(tn ) = y(tn+1

−h

n)

= y(tn+1 )

−

h 2n  hn .y (tn+1 ) + y (tn+1 ) + ... 2 

2

y(tn+1 )

−h

n .f (tn+1 , y(tn+1 ))

= y(tn )

− h2 y n



(tn+1 ) + ...

⇒

54


de la que, despreciando los términos cuadráticos y posteriores en hn, se obtiene la fórmula que constituye el esquema de cálculo del método de Euler impl´ıcito (o retrógrado): yn+1 = y n + hn .f (tn+1 , y(tn+1 )) Este esquema de cálculo exige, en cada paso, resolver una ecuación algebraica. Para ello puede utilizarse alguno de los métodos para resolver ecuaciones no lineales. eto dos. M´ etodo de Crank-Nicholsol θ -m´

Si se combinan los desarrollos en serie anteriores, ponderando el segundo por el parámetro θε[0, 1] y el primero por (1-θ) se obtendr´ıa la familia de θ-métodos. yn+1 = y n + hn θ)f (tn , y(tn )) + θhn f (tn+1 , y(tn+1 )) Para el caso en que θ=0 se trata del método de Euler expl´ıcito, para el caso θ=1 el de Euler impl´ıcito y para θ=0.5 se conoce como método de Crank-Nicholson: hn (f (tn , y(tn )) + f (tn+1 , y(tn+1 ))) 2 Por regla general los métodos expl´ıcitos necesitan más cantidad de pasos, pero estos son menos costosos. En los métodos impl´ıcitos, al tener que resolver una ecuaci´ on no lineal en cada paso, éstos son m´ as costosos, pero se pueden realizar en menos pasos. Una buena opción puede ser el m´ etodo de CrankNicholson, aunque no hay una regla fija en este aspecto, pudiendo para una ecuación ser más ventajoso uno y para otra ecuación otro. Otro aspecto a tener en cuenta es la estabilidad. El método explicito tiene pérdida de estabilidad, y eso no ocurre en los métodos impl´ıcitos. En el ejemplo 2, veremos la estabilidad, antes trabajaremos una ecuación con los tres métodos. En el libro de Quarteroni se dan ficheros para estos tres m´ etodos vistos, sus nombres son feuler.m (para Euler expl´ıcito), beuler.m (para el impl´ıcito) y cranknic.m . yn+1 = y n +

Para llevarlo a cabo con Octave, primeramente definimos la función f (t, y) en la variables t, y como una función inline, el intervalo en en el que se presenta el problema [ t0 , tN ] y la condición inicial y 0 como: f = inline(’...’,’t,y’) tspan = [t0,tN] y0 = ...

Después se especifica el número de intervalos Nh del mallado utilizado, es decir, el número de intervalos en los que se divide el intervalo [t0 , tN ] una vez conocido el paso de discretización h: Nh = ...

y resolvemos el PVI por medio de alguno de los siguientes métodos:

M´ etodo de Euler Expl´ ıcito M´ etodo de Euler Impl´ ıcito Crank-Nicholson rk2 rk3

[t,y] [t,y] [t,y] [t,y] [t,y]

= = = = =

feuler(f,tspan,y0,Nh) beuler(f,tspan,y0,Nh) cranknic(f,tspan,y0,Nh) rk2(f,tspan,y0,Nh) rk3(f,tspan,y0,Nh)

Para todos los métodos, la salida consiste en el vector t que contiene los N h + 1 puntos del mallado en los que se ha dividido el intervalo [t0 , tN ], y en el vector y que contiene las aproximaciones a la solución

55


exacta y (t) en los puntos del mallado. Una vez ejecutado el método podemos dibujar la aproximaci´ on mediante el comando plot(t,y) Ejemplo.



y  (t) = (1 + 6t) e1−y y(0) = log(10) + 1

− 6 t ∈ [0, 6]

La solución anal´ıtica es: y(t) = log(10e6t + t) + 1. Aplicar el algoritmo de Euler Expl´ıcito para calcular la solución en el intervalo temporal [0, 6] con pasos de discretización h1 = 0.5, h2 = 0.15 y h3 = 0.25. Discutir el comportamiento de las soluciones obtenidas en términos de estabilidad num´ erica. Aplicar los algoritmos de Euler Impl´ıcito y Crank-Nicholson Dibujar la soluci´ on anal´ıtica junto con las soluciones estables calculadas anteriormente. Calcular anal´ıticamente el instante en el cual se tiene el m´ınimo de concentraci´ on, as´ı como su valor. En primer lugar definimos la función y la solución anal´ıtica dada, as´ı como su valor inicial: >>f=inline(’(1+6*t)*exp(1-y)-6’,’t’,’y’); >>sol=inline(’log(10.*exp(-6.*t)+t)+1’,’t’); >> y0=log(10)+1; A continuación el intervalo y los hi >> tspan=[0,6]; h1=0.5; h2=0.15; h3=0.25; con ellos calculamos el valor de Ni correspondiente >> N1=(tspan(2)-tspan(1))/h1; >> N2=(tspan(2)-tspan(1))/h2; >> N3=(tspan(2)-tspan(1))/h3; y ya podemos aplicar el método de Euler expl´ıcito con feuler >> [t1,uf1]=feuler(f,tspan,y0,N1); >> [t2,uf2]=feuler(f,tspan,y0,N2); >> [t3,uf3]=feuler(f,tspan,y0,N3); Representamos la primera figura con los tres resultados y la solución anal´ıtica >> figure(1) >> plot(t1,uf1,’.-r’,t2,uf2,’b’,t3,uf3,’*g’,t2,sol(t2),’k’) >>legend(’ uf1’,’ uf2’,’ uf3’,’ sol’) obtenemos Comparamos algunos valores


56

Observamos que la solución obtenida para h1 presenta inestabilidad, vemos que las soluciones uf2 y uf2 se mantienen por debajo de la solución (eso ocurre con todas las expl´ıcitas) y parece ser que uf2 es mejor aproximación que uf3. Aplicamos ahora Euler impl´ıcito y Crank Nicholson: >> [t1,ub1]=beuler(f,tspan,y0,N1); >> [t2,ub2]=beuler(f,tspan,y0,N2); >> [t3,ub3]=beuler(f,tspan,y0,N3); >> figure(2); >> plot(t1,ub1,’.-r’,t2,ub2,’ b’,t3,ub3,’*g’,t2,sol(t2),’ k’) ; >> legend(’ ub1’,’ ub2’,’ ub3’,’ sol’) >> [t1,uc1]=cranknic(f,tspan,y0,N1); >> [t2,uc2]=cranknic(f,tspan,y0,N2); >> [t3,uc3]=cranknic(f,tspan,y0,N3); >> figure(3) >> plot(t1,uc1,’.-r’,t2,uc2, b’,t3,uc3,’*g’,t2,sol(t2),’ k’); >> legend(’ uc1’,’ uc2’,’ uc3’,’ sol’)

57


Calculamos el punto donde la función anal´ıtica alcanza el m´ınimo usando el paquete symbolic >> pkg load symbolic >> syms t >> y=log(10*exp(-6*t)+t)+1 >> df=diff(y,t)v >> solve(1-60*exp(-6*t)==0,t) La soluci´ on que nos da es

log(60)

6

>> log(60)/6 Que nos da el valor aproximado, 0.68239 Podemos comprobar que Octave symbolic no resuelve esta ecuación, y si lo intentamos con lsode, calculamos la aproximación y la representamos junto con la solución anal´ıtica. >> figure(5) >> x0=1; >> t=linspace(0,6); >> y=lsode(f,x0,t); >> ys=sol(t); >> plot(t,y,’b’,t,ys,´ or’)


58

Vemos que la solución dada por lsode no es nada buena, toma valores muy grandes, de forma que al representarla junto con la anal´ıtica, ésta no se aprecia. Observamos que ninguna solución presenta inestabilidad y que todas están por encima de la solución anal´ıtica, ambas cuestiones ocurre siempre en un método impl´ıcito. Nuevamente la mejor aproximación es la correspondiente a h2 y la peor la de h1. Como las aproximaciones obtenidas por Euler expl´ıcito estaban por debajo y por Euler impl´ıcito por encima, es probable que la mejor aproximaci´ on se obtenga por Crank-Nicholson. Comparamos la mejor aproximación de cada uno de los métodos, y efectivamente vemos que la mejor opci´ on es Crank-Nicholson: >> figure(4) >> plot(t2,uf2,’.-r’,t2,ub2,’b’,t2,uc2,’g’,t2,sol(t2),’k’); >> legend(’ uf2’,’ ub2’,’uc2’,’sol’)

59


Ejemplo. Problemas de Población. Desintegración radiactiva. Cuando el crecimiento de una población es directamente proporcional al tamaño de ésta en cada momento, el planteamiento matem´ atico a través  de una ecuación diferencial es: y = k y, siendo y(0) la poblaci´ on inicial. Es aplicable al tama˜ no de una poblaci´ on, a interés compuesto en inversiones bancarias, a problemas de desintegraci´ on radiactiva, etc... En problemas de radiación radiactiva como los del Carbono 14, usado para datar f´ osiles, es importante lo que llamamos semivida, que es el tiempo necesario para que la población se reduzca a la mitad. Cuando se trata de un crecimiento negativo, o decrecimiento, k es un valor negativo, aunque en muchas ocasiones se suele plantear la ecuaci´ on y  = k y, con k > 0. Analicemos esta ecuación con y (0) = 1 y en esta ecuación, dada su sencillez, es posible estudiar a priori fácilmente la estabilidad de las soluciones.

−

La solución anal´ıtica es y(t) = e −kt . Aplicaremos el algoritmo de Euler Expl´ıcito, y analizaremos la estabilidad seg´ un el valor de h. Estudiaremos tamibén los m´ etodos de Euler impl´ıcito y θ-métodos. Y particularizaremos a k = 5. Tomaremos h n = h la misma en cada paso. Euler Expl´ıcito yn+1 = y n + h ( ky n ) = (1

−

yn+1 = (1

− k h) y

n =

(1

− k h)2 y

− k h) y

n

n−1 = ... =

Euler Impl´ıcito yn+1 = y n + h ( ky n+1 )

−

yn+1 = (

⇒

yn+1 = (

(1

− k h) +1y0 n

1 ) yn 1+kh

1 )n+1 y0 1+kh

θ-métodos yn+1 = y n

− h ((1 − θ) k y

n +

yn+1 = (

1

θ k yn+1 )

⇒

yn+1 =

− (1 − θ) k h ) +1y0 n

1+θkh

1

− (1 − θ) k h y 1+θkh

n

60


Veamos cuándo estas soluciones son decrecientes y positivas. Para el método de Euler expl´ıcito, llamamos α=(1-k h), la soluci´ on aproximada en cada instante de tiempo t n+1 es la condición inicial (y(0) = y 0 = 1) multiplicada por α n+1 , es decir escribimos el esquema de Euler expl´ıcito en la forma: yn = α n y0 = α n Para que la solución aproximada efectivamente decrezca (al menos en valor absoluto) deberá verificarse que α < 1. Puesto que hemos supuesto k > 0 y la longitud del paso tambi´ en es positiva, un peque˜ no cálculo nos conduce a que el decaimiento en los valores absolutos de la solución se preservará si 0 < k h < 2 Lo anterior nos limita el tama˜ no del paso que podemos utilizar si deseamos que nuestra solución aproximada decrezca con el transcurrir del tiempo, debiendo verificarse para ello que h < k2 En el caso de que h = k2 la solución numérica no explotará con el transcurrir del tiempo pero como en ese caso α = 1 la solución numérica oscilará en cada instante de cálculo pasando del valor -1 al valor 1. No obstante se puede ser un poco más sutil en el examen de la longitud del paso de cálculo a utilizar. En efecto, si se desea que además la solución aproximada sea siempre positiva la condición a imponer es que 0 < α < 1 , lo que nos conduce a que hε(0, k1 ).

| |

−

N´ otese que el intervalo al que deben pertenecer los valores de h es un intervalo abierto. En efecto, en sus extremos la solución numérica no oscilará entre valores negativos y positivos, pero tambi´ en adolecer´ a de graves defectos. As´ı el valor h = 0 no puede tomarse (pues todos los puntos de cálculo se reducir´ıan al instante inicial t0 ). Y si se tomase h = k1 el valor de α ser´ıa nulo y la solución aproximada en instantes positivos será siempre 0. El problema entonces es saber cuál es el mejor valor que podr´ıamos tomar para el paso de integración. El responder a esta pregunta, en este caso concreto, es muy simple. Puesto que la soluci´ on anal´ıtica es conocida, se tendr´ıa un error nulo si: (1

− h k)

n

= e tn = e nh

h

⇒ 1 − h.k = e ⇒ h = 0

Pero es absurdo pensar en tomar h = 0. No obstante, lo anterior nos indica que cuanto menor sea el paso de integración menor será el error cometido. Debe matizarse no obstante que aqu´ı nos estamos refiriendo sólo a un tipo de error: el que se comete con la aplicación del método sin redondear los valores obtenidos. La consideración de que al trabajar con un programa de cálculo se introducen además errores de redondeo al no poderse almacenar infinitos decimales introducirá (como se verá m´ as adelante) un l´ımite inferior a los valores de h que puedan considerarse. Pero olvidándonos por el momento de los errores de redondeo, efectivamente, cuanto menor sea el tamaño del paso de integració n m´ as precisa ser´ a la aproximación obtenida. Ilustremos lo anterior considerando la solución anal´ıtica del problema y las soluciones aproximadas mediante el m´ etodo de Euler para distintos valores del tama˜ no del paso de integraci´ on y tomando como valor de k = 5 (lo que hace que para pasos menores que 1 /5 = 0,2 el método produzca soluciones siempre positivas que decrecen en valor absoluto con el paso del tiempo): Obsérvese como para h = 0,2 la solución obtenida es, efectivamente, siempre nula y como la solución obtenida con h = 0,02 es más precisa que la obtenida para h = 0,1, y esta, a su vez, es más precisa que la obtenida para h = 0,2. Si tom´ asemos mayores valores del paso de integración comenzar´ıan a aparecer soluciones oscilantes pero si no sobrepasamos el valor cr´ıtico h = 2/k = 0,4 la solución decrece en valor absoluto por lo que para tiempos grandes acaba convergiendo hacia la solución exacta. Si a h le asignamos el valor cr´ıtico h = 0,4 la solución permanece acotada, no explota para tiempos grandes, pero tampoco decrece en valor absoluto y oscila entre -1 y 1. Finalmente si a h le asignamos valores superiores a 0.4, por ejemplo h = 0,5 la solución aproximada, en valor absoluto, crecerá por lo que tenderá a explotar para tiempos elevados. Diremos en este caso que el esquema considerado se hace inestable Representamos las gráficas uf1 correspondiente a h = 0,2, uf2 para h = 0,1 y uf3 para h = 0,01 para comprobar lo explicado anteriormente. En otra gráfica representamos uf4 correspondiente a h = 0,4 y uf5 para h = 0,5 para mostrar los casos de inestabilidad. En ambos casos aparece la solución.

61


Examinemos ahora los métodos impl´ıcitos (es decir con θ > 0). En este caso la solución aproximada también puede expresarse como: yn = α n y0 (n = 1, 2, ....,N) siendo ahora: α =

1

− (1 − θ) k h = 1 + θ k h − k h = 1 − 1+θkh

1+θkh

kh 1+θkh

Puesto que 0 < θ 1, y k > 0 se verificará que a siempre será inferior a 1. Para garantizar que tambi´ en es superior a -1 se debe verificar que:

≤

kh < 2 1+θkh

⇒ (1 − 2 θ) k h < 2

En el caso de que θ = 0, 5 la expresión anterior se verificará para cualquier valor que se tome para h el esquema en ese caso es incondicionalmente estable. Y si θ < 0, 5 la condición a imponer sobre el paso de integraci´ on para garantizar la estabilidad, será: h<

2 k (1

− 2 θ)

Estas elecciones del valor del paso de integración nos aseguran que el valor de la solución aproximada que se obtenga tambi´ en va a decaer hacia 0 para tiempos grandes. Si adem´ as se desease asegurar la positividad de las soluciones aproximadas se deber´ıa hilar un poco más fino. En efecto para que esto suceda se debe obligar a que 0 < α < 1. Ello nos conduce a que: kh < 1 1+θkh

⇒ (1 − θ) k h < 1 ⇒

h<

2 k (1 θ)

−

62


Obsérvese ervese que cuanto más as se aproxime ? al valor 1 mayor será el tama˜ no no de paso máximo aximo que puede ser tomado tomado para para garan garantiz tizar ar tanto tanto la estabi estabilid lidad ad de la soluci soluci´ oń aproximada (es decir que no explote on para tiempos grandes) como su positividad. Particularmente para θ = 1 (esquema de Euler impl´ impl´ıcito) cualquier tama˜ no del paso de integración no on nos servir servir´ıa para asegurar asegurar ambas ambas cosas. cosas. Para Para el esquema esquema de Crank Nicholson la estabilidad nos la garantizar´ garantizar´ıa cualquier tama no n ˜ o del paso pero la positividad nos obligar´ıa ıa a considerar consider ar tama˜ tama nos n ˜ os inferiores a 2/k. Ello nos ampl´ ampl´ıa el tamaño no del paso que puede ser tomado respecto al esquema expl´ expl´ıcito (en el que la positividad p ositividad y la estabilidad se aseguraban con tama˜ nos nos inferiores a 1/k). Ilustremos lo anterior con la aplicación on de los esquemas de Euler impl´ impl´ıcito y de Crank-Nicholson al problema:



y (t) = 5 y y(0) = 1

−

t > 0

para distintos valores del paso de integración. on. Comencemos Comencem os con el método etod o de d e Euler E uler impl´ıcito ıcito que, en este caso se formulará como: y0 = 1 1 yn = ( 1+5 )n (n = 1, 2, 3...) ...) h y que será estable y positivo para todos los valores de h que se consideren. Como puede apreciarse las soluciones, efectivamente siempre permanecen p ositivas y decrecen con el tiempo. Pero también en puede advertirse como su precisión o n es más as pobre a medida que se aumenta el paso de integración temporal. Adem´ as as la solución on es aproximada en exceso es decir con valores mayores que la exacta, a diferencia de lo que suced´ıa ıa en el método etod o expl´ıcito. ıcito. Los resultados resulta dos que se obtienen obtien en con el método etod o de Euler Impl´ıcito ıcito y con el método etod o de d e Crank-Nichol Cr ank-Nicholson: son: y0 = 1 −0,5 5 h n yn = ( 11+0 ) ,5 5 h

(n = 1, 2, 3...) ...)

Procediendo como se hizo anteriormente anteriormente para Euler expl´ expl´ıcito, para los mismos valores de hi, se obtendr´ ten dr´ıan ıan las gráficas aficas

63


Se observa que todas las soluciones permanecen positivas y, además, as, m´ etodo etodo no parece tan sensible al paso de integración on en cuanto a su precisión. on. No obstante si sobrepasamos el l´ımite h = 2/k = 2/5 = 0,4, (ub4, uc4) especialmente en el método etodo de Crank Nicholson, para el tama˜ no del paso comenzaremos no a perder la positividad y a notar la influencia del valor de h. M´ etodo etod o de Heun, o M´ etodo etod o de Predicci´ Predi cci´ on on Correcci´ on on

Esta cualidad de los esquemas impl´ıcitos, ıcitos, poder p oder considerar pasos de integración on m´ as as amplios que en los expl´ıcitos ıcitos sin perder la estabilidad, permite pensar en abordar problemas que se planteen en dominios [t0 , t0 + T] en los que T tome valores valores elevados. elevados. Piénsese, ensese, por ejemplo, ejemplo, que la EDO y = ky es la que rige procesos de desintegración on de isótopos otopos radiactivos que, por su peligrosidad hacia el entorno, pueden requerir estudiarlos durante millones de años nos lo l o que har´ har´ıa inviable invi able la l a consideraci cons ideraci´ón on de pasos de integración on excesivamente excesivamente peque˜ peque nos. n ˜ os. Pero, lamentablemente, no todo es la estabilidad y la positividad de los esquemas de cálculo alculo sino que, además, as, habr´ a que prestar atención on a su precisión on se hará un poco m´ as as adelante. Y a otros problemas como el que se plantea en el ejemplo siguiente.

−

Ejemplo



y  (t) = y 3 y (0) = 1

−

t > 0

La soluci´ soluci´ on on anal an al´´ıtica, es de variable separada, s eparada, es y(t) =

1 . 1+t

Las iteraci ite raciones ones de Euler Eule r impl´ im pl´ıcito ıci to ser´ıan ıan y n+1 + h yn2 +1 = y n Para h=0.1, calculamos el valor aproximado para t1 = 0,1. y1 + 0,1 y12 = 1 -10.916 ó 0.916 0.916079 079

⇒

y1 puede ser

Conociendo la solución on anal´ıtica ıtica es fácil acil comprobar que la solución on correcta es la segunda, pero lo normal es no conocer la solución on exacta. ¿Cómo omo resolver este problema? En ausencia de criterios propios

64


del problema no se puede descartar una solución on frente a otra. Por ello, en estos casos, lo aconsejable es acudir a los m´ etodos etodos predictor-corrector en los que, en s´ıntesis, se utiliza un m´ etodo etodo expl´ expl´ıcito para (predecir) obtener el valor de partida con el que aplicar el m´ etodo etodo de resolución on de ecuaciones no lineales a la ecuación on obtenida del esquema impl´ impl´ıcito que nos permite refinar (corregir) el valor aproximado de la soluci´ on finalmente obtenida. La aplicación on on de esta forma de proceder pro ceder a nuestro problema nos conducir´ conducir´ıa en el primer paso de integración on a: Fase de predicción on mediante Euler expl´ıcito: ıcito: y1(0) = y 0

− h y02 = 1 − 0,1(1)2 = 0,9

Fase de corrección on mediante me diante Euler Eu ler impl i mpl´´ıcito (resolviendo (resolv iendo por el método etod o de aproximaciones aproximaci ones sucesivas): s ucesivas): y1(

iter+1)

= y 0

iter 2

− h (y1

)

(iter = iter = 0, 1, 2,...) ,...)

que nos conduce a: y1(1) = 1 0,1(0, 1(0,9)2 = 0, 0 ,919 (2) y1 = 1 0,1 (0, (0,919)2 = 0,9155439 y1(3) = 1 0,1 (0, (0,9155439)2 = 0, 0 ,916177 (4) 2 y1 = 1 0,1 (0, (0,916177) = 0, 0 ,9160617 (5) 2 y1 = 1 0,1 (0, (0,9160617) = 0, 0 ,9160838 (6) 2 y1 = 1 0,1 (0, (0,9160838) = 0, 0 ,9160792 (7) y1 = 1 0,1 (0, (0,9160792)2 = 0, 0 ,9155439

− − − − − −

−

por lo que tomaremos como y 1 = 0,9160798 9160798.. En el segundo paso de integración on se tendr´ıa: ıa: Fase de predicci´ predicci ón on mediante Euler expl´ıcito: ıcito: y2(0) = y 1

(0,9160798)2 = 0,8321595 − h y12 = 0,0 ,9160798 − 0,1 (0,

Fase de corrección on mediante Euler impl´ıcito ıcito (iter+1)

y2

= y 1

iter 2

− h (y2

)

(iter = iter = 0, 1, 2,...) ,...)

que nos conduce a: (1)

y2 (2) y2 (3) y2 (4) y2 (5) y2 (6) y2 (7) y2

= 0,9160798 = 0,9160798 = 0,9160798 = 0,9160798 = 0,9160798 = 0,9160798 = 0,9160798

− 0,1(0, 1(0,8321595)2 = 0,8468308 − 0,1 (0, (0,8468308)2 = 0,8443675 − 0,1 (0, (0,8443675)2 = 0,8447841 − 0,1 (0, (0,8447841)2 = 0,8447137 − 0,1 (0, (0,8447137)2 = 0,8447256 − 0,1 (0, (0,8447256)2 = 0,8447236 − 0,1 (0, (0,8447236)2 = 0,8447239

por lo que tomaremos como y 2 = 0,8447239 8447239.. El tercer paso de integración on nos permitir permi tir´´ıa calcular y3 mediante Euler expl´ıcito: ıcito:

y (0 − 3) consistir consist ir´´ıa en: Fase de predicci´ predicci ón on ≈ y(0

y3(0) = y 2

(0,8447239)2 = 0,7733681 − h y22 = 0,0 ,8447239 − 0,1 (0, Fase de corrección on mediante Euler impl´ıcito ıcito y3( +1) = y 1 − h (y3 )2 (iter = iter = 0, 1, 2,...) ,...) iter

iter

65


de lo que obtendr´ıamos y3(1) = 0,8447239 y3(2) = 0,8447239 y3(3) = 0,8447239 y3(4) = 0,8447239 y3(5) = 0,8447239 y3(6) = 0,8447239 y3(7) = 0,8447239

− 0,1(0,7733681)2 = 0,7849141 − 0,1 (0,7849141)2 = 0,7831149 − 0,1 (0,7831149)2 = 0,783397 − 0,1 (0,783397)2 = 0,7833529 − 0,1 (0,7833529)2 = 0,7833598 − 0,1 (0,7833598)2 = 0,7833587 − 0,1 (0,7833587)2 = 0,7833589

por lo que y 3 = 0,7833589. Esta estrategia puede extenderse a esquemas numéricos muy diferentes, combinando un esquema expl´ıcito de predicción y un esquema impl´ıcito de correcci´ on. Por ejemplo, si se deseara utilizar un θ -método con θ > 0 para la fase correctora y el m´ etodo de Euler expl´ıcito para la fase predictora, el esquema de cálculo ser´ıa: Fase predictora (Euler expl´ıcito): yn(0) +1 = y 0

−h

n

f (tn , yn )

Fase correctora ( θ -método combinado con aproximaciones sucesivas): yn( +1 +1) = y n + (1 iter

− θ) h

f (tn , yn ) + θ hn f (tn+1 , yn( +1 ) iter

n

(iter = 0, 1, 2,...)

As´ı por ejemplo, se deja como ejercicio propuesto verificar que la estrategia de cálculo anterior con θ = 0,5 aplicada al problema de valor inicial del ejemplo 3 conduce a los valores y1 = 0,9087121, y2 = 0, 8327505, y3 = 0, 7685438, y4 = 0, 7135529, y5 = 0, 6659224, .... De hecho algunos métodos de cálculo muy utilizados en la práctica, y que posteriormente nos reencontraremos, podr´ıan interpretarse como un esquema predictor corrector en los que se realizan pocas iteraciones del m´ etodo de aproximaciones sucesivas en la fase de corrección. Es el caso por ejemplo de combinar el método de Euler expl´ıcito en la fase de predicci´ on y realizar una sóla iteración en la fase correctora utilizando el método de Crank-Nicholson, lo que nos conduce a: yn(0) +1 = y n + hn f (tn , yn ) yn+1 = y n +

h n (0) (f (tn , yn ) + f (tn+1 , yn+1 )) 2

por lo que, resumiendo las dos etapas en una única expresión resulta: yn+1 = y : n +

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

Este esquema de cálculo recibe el nombre de método de Heun y su generalizaci´ on ser´ıan los métodos de Runge-Kutta.

7.3.

Sistemas de Ecuaciones Diferenciales

Los P.V.I. pueden formularse de una forma más general para sistemas de ecuaciones diferenciales de primer orden. Más concretamente, siendo I un intervalo de IR de la forma [t0 , t0 + T ] y siendo f una función continua definida en IxRm y con valores en R m .

66


La extensi´ on de los m´ etodos num´ ericos a problemas en los que intervenga un sistema de ecuaciones diferenciales ordinarias es automática y, en cuanto a su análisis, serán v´ alidos los teoremas y propiedades que se presenten con la precaución de emplear la métrica (normas) correspondiente en IRm.

7.3.1.

Ecuaci´ on de Orden Superior

Cabe se˜ nalar además que numerosos problemas de valor inicial se formulan mediante ecuaciones diferenciales ordinarias de orden superior a 1 de la forma:

  

y ( p (t) = f (t, y(t), y (t),...,y ( p−1 ), y(t0 ) = y 0  y (t0 ) = y 0(2) .................. ( p−1) p − y ( 1 (t0 ) = y 0

tε[t0 , t0 + T ]

Tales tipos de P.V.I de orden superior a 1 pueden ser reducidos a sistemas de p ecuaciones diferenciales de primer orden denominando: z1 (t) = y(t), z2 (t) = y  (t), z2 (t) = y  (t),....,z p (t) = y ( p−1 (t)

7.3.2.

Usando Symbolic. Con Dsolve

Con el paquete Symbolic, dependerá de la versión de Phyton que tengamos instalada en el ordenador, por eso voy a hacerlo en octave online. Si lo intentamos directamente en nuestro octave puede que nos dé problemas y sólo nos resuelva los homog´ eneos. https://octave-online.net/ Con los sistemas lineales no tendremos problemas. Por ejemplo, para resolver el sistema



2t

x = 2 x y e2t + cose 3 (t) y  = 17 x + 6 y + 25 t2

− − −

67


Podemos probar en nuestro programa resolver la homogenea y la hará, aunque no la completa >> A = [ 2, 1, 2; 1, 0, 2; 1, 1, 1]; roots(poly(A))

−

−

−

>> pkg load symbolic >> syms x(t) y(t) >> e1= diff(x,t)==-2*x-y >> e2= diff(y,t)==17*x+6*y >> dsolve([e1,e2]) Pero con sistemas no lineales no funciona. Veamos con este ejemplo de una reacción qu´ımica: Sea la ecuación qu´ımica A+B

↔ C → D + E

donde A,B,C,D, y E representan distintos compuestos que reaccionan entre s´ı. En esta ecuación, la concentraci´ on de la especie A, que llamaremos y 1 , y de la especie C , que denotaremos y 2 reaccionan entre

68


si mediante el sistema de ecuaciones diferenciales



dy1 dt dy2 dt

= =

−y12 + y2 −2y2 + y12

,

teniendo en cuenta las condiciones iniciales y 1 (0) = 1 e y 2 (0) = 0 Utilizando Octave Online ser´ıa

Como vemos no la resuelve.

7.3.3.

Resoluci´ on con lsode

Intentamos resolver el problema de la reacción qu´ımica con lsode Definimos la funci´ on, los valores iniciales y el intervalo entre 0 y 60 con h = 0,1 >> function xdot = f(y,t)

−y(1)2 + y(2); >> xdot(2) = y(1)2 − 2 ∗ y(2); >> xdot(1) =

>> endfunction >> y 0 = [0;1];


69

>> t = [0 : 0,1 : 60]; Aplicamos la orden Isode >> y=lsode(’f’, y0 , t);

7.3.4.

Resoluci´ on por M´ etodos Num´ ericos

Resolvamos ahora el sistema anterior aplicando los otros métodos. Empezemos con Euler Expl´ıcito. El planteamiento es el mismo, pero la función y los datos serán vectores. En primer lugar definimos la función, para ellos utilizamos la función function / endfunction : >> function R=REACCION(t,R) K1=1; K2=1; K3=1; K=0; RR(1)=-K1*R(1)*(R(1)-K)+K2*R(2); RR(2)=-(K2+K3)*R(2)+K1*R(1)*(R(1)-K); R=RR >> endfunction Se debe utilizar una variable auxiliar, en este caso RR, ya que sino, si hubiesemos puesto R(1)=-K1*R(1)*(R(1)-K)+K2*R(2);

70


R(2)=-(K2+K3)*R(2)+K1*R(1)*(R(1)-K); al calcular R(2) ya R(1) se habr´ıa actualizado en el c´ alculo de la l´ınea anterior. Definimos los valores iniciales para este problema: >> RA0=1 >> RC0=0; >> R 0 = [RA0 RC 0]; >> t

0=0;

>> t

f=60;

>> h=0.1; >> N=(t f-t

0)/h;

Y aplicamos la orden feuler >> [t,Re]=feuler(’REACCION’,[t

0, t

f],R 0,N);

Pulsamos q para que no salgan todos los ficheros. Podemos aplicar tambi´ en del mismo modo Euler Impl´ıcito y Crank-Nicholson >> [t,Ri]=beuler(’REACCION’,[t 0, t

f],R 0,N);

>> [t,Rc]=cranknic(’REACCION’,[t 0, t

f],R

0,N);

Representamos las tres juntas >> figure1 >> plot(t,Re,’-.’,t,Ri,’-’,t,Rc,’ o’) >> title(’ Euler Expl´ıcito,Impl´ıcito y Crank Nicholson. Concentracion de las especies’) >> legend(’ A exp’,’ C exp’,’ A imp’,’ C imp’,’ A exp’,’ C Cran’);

71


Obtenemos

En el caso de los sistemas no se cumple que el en del caso de que Euler Expl´ıcito e Impl´ıcito sean convergente, Crank-Nicholson estará entre ambos, al menos en todas las variables. Ejemplo. Modelo de Epidemias. Sea la ecuación

 

x = βxy y = βxy γy , donde z  = γy

−

−

x(t) representa a los individuos susceptibles, es decir, aquellos que no han enfermado anteriormente y por lo tanto pueden resultar infectados al entrar en contacto con la enfermedad. El tiempo medido en semanas y(t) representa a los individuos infectados y por lo tanto en condiciones de transmitir la enfermedad a los del grupo anterior x z(t), por u ´ ltimo, representa a los individuos recobrados de la enfermedad, y que ya no están en condiciones ni de enfermar nuevamente ni de transmitir la enfermedad a otros. Resolver para β = 0,2 y γ = 0,7 en el intervalo [0, 5] sabiendo que al principio hab´ıa una persona infectada en una población de 20 habitantes. Qué situación habrá en la semana 3. Con LSODE Definimos la función para esos valores >> function sis = f(y,t) b=0.2,g=0.7; sis(1) = -b*y(1)*y(2); sis(2) = b*y(1)*y(2)-g*y(2); sis(3)=g*y(2) >>endfunction Damos los valores iniciales, en un principio habrá 1 enfermo, 19 sanos y nadie que haya pasado la enfermedad, luego x = 19, y = 1, z = 0

72


>> y0=[19; 1;0]; Tomamos el intervalo y h = 0,1 >> t=[0:0.1:5]; >> y=lsode(”f”, y0 , t); sis = -3.80000 3.10000 0.70000 sis = -3.80000 3.10000 0.70000 —————————————— Pulsamos q para no seguir viendo los valores. Representamos con plot. >> figure(1) >> plot(t,y) Y obtenemos

Observando los valores iniciales, vemos que x está en az´ ul, y en rojo y z en amarillo. Si queremos asignar nosotros los colores: >>figure 2; plot (t, y(:,1),’b’, t, y(:,2), ’r’,t,y(:,3),’y’);legend(’ y En la semana 3, tendremos >> t(31)

1(t)’,’ y

2(t)’,’ y

3(t)’)

73


ans = 3 >> y(:,1)(31) ans = 0.27486 >> y(:,2)(31) ans = 4.8994 >> y(:,3)(31) ans = 14.826 Habr´ a 5 enfermos y 15 que han pasado la enfermedad. Resoluci´ on por Métodos Numéricos: Definimos la función con su variable auxiliar: >>function R=EPIDEMIA3(t,R) b=0.2; g=0.7; RR(1)=-b*R(1)*R(2); RR(2)=b*R(1)*R(2)-g*R(2); RR(3)=g*R(2); R=RR >>endfunction Ahora los valores iniciales >>R

0=[19 1 0];t 0=0;t

f=5; h=0.1;N=t

´ >>[t,Re]=feuler(EPIDEMIA3’,[t

0, t

f],R

f/h; 0,N);

>>figure 2; plot(t,Re);title(’ Euler Expl´ıcito. Epidemia 3x3’) Obtenemos

Probamos con Euler impl´ıcito

74


>>[t,Ri]=beuler(’ EPIDEMIA3’,[t

0, t

f],R 0,N);

>>figure 3; plot(t,Rb); title(’ Euler Impl´ıcito. Epidemia 3x3’)

Como vemos, en este caso Euler Impl´ıcito no da buenos resultados, y por tanto, tampoco lo dar´ a Crank-Nicholson. Resolvamos una ecuación de oden superior no lineal, en el estudio de algunas capas l´ımite de flujos laminares aparece la denominada ecuación de Blasius que es una EDO de tercer orden de la forma: y  (x) + y(x) y  (x) = 0

x

≥0

Se ha llamado a la variable independiente x porque en esta ocasión la variable independiente no es el tiempo sino que es una variable adimensional relacionada con la coordenada espacial. La ecuación de Blasius se acompaña de las tres condiciones iniciales y(0) = y  (0) = 0, y (0) = α donde α es un valor conocido α = 0,47. Con ello se puede plantear el problema de valor inicial de Blasius en forma de sistema, denotando z 1 (x) = y(x), z2 (x) = y  (x), z 3 (x) = y  (x), de la forma siguiente:

 

z1 (x) = z 2 (x) z2 (x) = z 3 (x) z3 (x) = z1 (x) z3 (x) z1 (0) = o z2 (0) = 0 z3 (0) = α

−

x x x

≥0 ≥0 ≥0

o expresado de manera matricial



z  (x) = f (x, z(x)) z(0) = [0, 0, α]

x

≥0

donde z(x) =[z1 (x), z2 (x), z3 (x)] y f (x, z(x)) = [z2 (x), z3 (x), z1 (x)z3 (x)]

−

Realizamos la resolución en fichero resolucionblasius.m, donde se llama a la función blasius, definida en el fichero ejercicioblasius.m


ejercicioblasius.m: function Z=blasius(x,Z) Z(1)= Z(2); Z(2)=Z(3) Z(3)= -Z(1)*Z(3); return resolucionblasius.m: alpha=0.47 z0=[0 0 alpha]; x0=0; xf=6; h=0.1; N=xf/h; [x,Z]=feuler(’fblasius’, [x0,xf],z0,N); umbral=inline(’0.*x+1’,’x’); figure plot(x,Z,x,umbral(x),’.g’); title(’perfil blasius’) Obteniendo

75

76


7.3.5.

M´ etodos de Runge-Kutta

Los m´ etodos de Runge-Kutta representan el ejemplo m´ as clásico de los m´ etodos de pasos libres. Consideramos el problema de valor inicial (P.V.I.) y la subdivisión del intervalo [t0 , tN ] generando los puntos t 0 < t1 < t2 < .... < tn < tn+1 < .... < tN , y con subintervalos [tn , tn+1 ] de longitud h n (n = 0, 1, ..., N - 1). Siendo y(tn ) el valor de la solució n en tn , el valor exacto de la solució n en tn+1 puede estimarse mediante la expresión: tn+1

y(tn+1 ) = y(tn ) +



f (t, y(t))dx.

tn

Una forma de obtener aproximaciones de dicho valor consistirá en aproximar la integral de la expresión anterior mediante una fórmula de integración numérica: p

yn+1 = y n + hn



aj f (tn,j , yn,j ),

j =1

donde p es el n´ umero de puntos de integración usados en la fórmula escogida, aj (j = 1, ..., p) son los pesos de dicha fórmula, t n,j (j = 1, ..., p) son los p puntos pertenecientes al intervalo [tn , tn+1 ] que act´ uan como soporte para la fórmula de integración y, finalmente, yn,j son los valores (o una aproximación de ellos) de y(t) en los puntos t n,j . El problema que se plantea es cómo evaluar y n,j . Para ello, procediendo de forma similar, se considerará que: tn,j

y(tn,j ) = y(tn ) +



f (t, y(t))dx

j = 1, 2,...,p,

tn

por lo que una aproximación de dichos valores puede obtenerse, nuevamente, aproximando la integral anterior. Y en los m´ etodos de Runge-Kutta esta aproximaci´ on se realiza utilizando los mismos puntos que en la expresión anterior, es decir: p

yn,j = y n + hn



bi,j f (tn,i , yn,i ),

j =1

donde bj,k (j, k = 1, ..., p) es el peso otorgado al punto t n,k en la fórmula de integración numérica utilizada para aproximar el valor de y(t) en tn,j . Las expresiones dadas forman un sistema de (p+1) ecuaciones en las que las incógnitas son los p valores y n,j (j = 1, ..., p) y el valor y n+1 . El uso de una u otra fórmula de integración numérica en las expresiones conduce a muy distintos esquemas de tipo Runge-Kutta. Por ello para definir una fórmula de Runge-Kutta se debe concretar cuáles son los puntos tn,j (j =1, ..., p) utilizados en las fórmulas de integración numérica, cuáles son los coeficientes aj y bj,k (j, k = 1, ..., p) usados en las fórmulas. Habitualmente, los puntos t n,j se suelen definir mediante una expresión del tipo: tn,j = t n + cj hn

j = 1,...,p

y el esquema se define mediante la tabla siguiente: c1 c2 . . c p

b1,1 b2,1 . . b p,1 a1

b1,2 b2,2 . . b p,2 a2

... b1,p ... b2,p ... . ... . ... b p,p . a p

77


En dicha tabla designaremos por vector a, matriz B y matriz C al vector y matrices siguientes:

a =

 

a1 a1 . a p

 

B =

  

b1,1 b2,1 . . b p,1

  

b1,2 ... b1,p b2,2 ... b2,p . ... . . ... . b p,2 ... b p,p

, B =

  

c1 0 . . 0

0 c2 . . 0

... 0 ... 0 ... . ... . ... c p

  

.

Ejemplo de Runge-Kutta. Puede dise˜ narse un método de Runge-Kutta en el que la f´ ormula que nos proporcione y n+1 utilice la f´ ormula de Simpson: b b a a+b φ(x)dx φ(a) + 4 φ( ) + φ(b) 6 2 a

≈ −







Para ello, se tomará p = 3 y en el intervalo [tn , tn+1 ] los puntos de integración serán: tn,1 = t n tn,2 = t n + 0,5 hn tn,3 = t n + hn

(

(→ c1 = 0) → c2 = 0,5) (→ c3 = 1)

y los pesos de integración en la f´ ormula considerada serán: a1 =1/6, a 2 =4/6, a 3 =1/6, resultando: yn+1 = y n+1 + hn





1 4 1 f (tn,1 , yn,1 ) + f (tn,2 , yn,2 ) + f (tn,3 , yn,3 ) 2 6 6

Para emplear la fórmula anterior es necesario evaluar: y n,1 (= y n ),yn,2 e y n,3 . Para ello se sabe que: yn,1 = y n yn,2 = y n + yn,3 = y n +

tn,2 tn tn,3 tn

 

f (t, y(t))dt f (t, y(t))dt

y evaluando la primera de las integrales mediante la fórmula del trapecio: tn,2



f (t, y(t))dt

tn

≈ t 2 2− t n,

n

(f (tn , yn ) + f (tn,2 , yn,2 )) =

hn (f (tn , yn ) + f (tn,2 , yn,2 ) 4

y la segunda mediante el método del punto medio: tn,3



tn

f (t, y(t))dt

≈ (t 3 − t n,

n ) f (tn,2 , yn,2 )

= h n f (tn,2 , yn,2 )

resultar´ a: yn,1 = y n yn,1 = y n + h4n (f (tn , yn ) + f (tn,2 , yn,2 )) yn,3 = y n + hn f (tn,2 , yn,2 )

( b1,1 = 0, b1,2 = 0, b1,3 = 0) b2,1 = 1/4, b2,2 = 1/4, b2,3 = 0) ( b3,1 = 0, b3,2 = 1, b3,3 = 0)

→ (→ →

La segunda de las expresiones anteriores es una ecuación (en general no lineal) que deberá resolverse mediante alg´ un m´ etodo numérico como los presentados en el tema anterior. El esquema en forma de tabla será: 0 0 0 0 1/2 1/4 1/4 0 1 0 1 0 1/6 4/6 1/6

78


El m´ etodo de Euler, puede considerarse como un caso particular de los métodos de Runge-Kutta en el que p = 1 y: 0 0 1 El algoritmo rk2 contiene el método de Heun, mientras que el algoritmo rk3 contiene el método de Runge-Kutta de orden 3 dado por la tabla 0 1/2 1

0 0 0 1/2 0 0 1 2 0 1/6 4/6 1/6

−

El método de Euler modificado y el método de Heun. Siendo α un n´ umero real, el método dado por: 0 α

0 α 1

0 0

− 21

1 2α

α

se traduce en las expresiones: yn,1 = y n yn,2 = y n + α hn f (tn , yn ) yn+1 = y n + hn (1 21α ) (f (tn , yn ) +



−

1 f (tn + 2α

α hn , yn,2 )

Para el caso α = 0,5, el método se denomina de Método de Euler modificado:







hn hn yn+1 = y n + hn f tn + , yn + f (tn + yn ) 2 2 Para el caso α = 1, el método se denomina de Método de Heun: yn+1 = y n +

h n (f (tn , yn ) + f (tn + hn , yn + hn f (tn + yn ))) 2

El ejemplo de Runge-Kutta clásico responde a la tabla: 0 1/2 1/2 1

0 0 0 0 1/2 0 0 0 0 1/2 0 0 0 0 1 0 1/6 1/3 1/3 1/6

Por tanto, en este método, los puntos intermedios escogidos en el paso n serán: tn,1 = t n tn,2 = t n + 0,5 hn tn,3 == t n + 0,5 hn tn,4 = t n + hn

(= t n,2 )

y como vemos se repite el punto medio. Con los puntos as´ı tomados, se considera y n,1 = y n . El valor de y n,2 se eval´ ua aproximando la integral correspondiente mediante la fórmula del rectángulo con soporte en el extremo izquierdo del intervalo de integración, es decir: tn,2

y(tn,2 ) = y(tn ) +



tn

f (t, y(t))dt

⇒y

n,2 = y n +

h n f (tn , yn ) 2

79


El valor de y n,3 se calcula aproximando la integral correspondiente mediante la fórmula del rectángulo pero esta vez con el punto de soporte en el extremo derecho del intervalo y tomando precisamente yn,2 como valor aproximado de la función y (t) en dicho extremo derecho: tn,3

y(tn,3 ) = y(tn ) +



f (t, y(t))dt

tn

⇒y

n,3 = y n +

hn h n f (tn + , yn,2 ) 2 2

En cuanto al valor de y n,4 lo evaluaremos aproximando la integral correspondiente mediante una fórmula de punto medio, considerando como valor de y (t) en el punto medio del intervalo en el que se plantea la f´ ormula el valor y n,3 que se acaba de calcular, es decir tn +hn

y(tn,4 ) = y(tn ) +



f (t, y(t))dt

⇒y

n,4 = y n +

tn

hn f (tn +

hn , yn,3 ) 2

Por u ´ ltimo, el valor de yn+1 se eval´ ua, con ayuda de los valores anteriores, aproximando la integral correspondiente mediante una fórmula de Simpson en la cual el peso asignado al punto medio del intervalo se reparte por igual entre las dos aproximaciones del valor de y(t) que hemos obtenido en dicho punto medio (y(tn,2 e y n,3 ): tn +hn

y(tn+1 ) = y(tn ) +



f (t, y(t))dt

tn

⇒

h n (f (tn , yn ) + 2 f (tn,2 , yn,2 ) + 2 f (tn,3 , yn,3 ) + f (tn,4 , yn,4 )) 6 h n hn hn = y n + (f (tn , yn ) + 2 (f (tn + , yn,2 ) + f (tn + , yn,3 )) + f (tn + hn , yn,4 )) 6 2 2 Método de Euler Impl´ıcito Mejorado es el método de Runge-Kutta definido mediante la tabla: yn+1 = y n +

1

1 1

que es el método yn,1 = y n + hn f (tn+1 , yn,1 ) yn+1 = y n + hn f (tn+1 , yn + hn f (tn+1 , yn,1 )) En resumen, el algoritmo del método de Runge-Kutta clásico, está definido, en cada paso, por las expresiones: tn,j = t n + cj hn ( j = 1,..,p) p

yn,j = y n + hn



bi,j f (tn,i , yn,i ),

( j = 1,..,p)

j =1

p

yn+1 = y n + hn



aj f (tn,j , yn,j )

j =1

A la hora de realizar un algoritmo de un m´ etodo como el anterior, computacionalmente es m´ as eficaz denotar por: p

W n,j = f (tn,j , yn + hn



bi,j W n,i ),

( j = 1,..,p)

j =1

con lo que:

p

yn+1 = y n + hn



aj W n,j .

j =1

Se trata de dos sistemas de ecuaciones, en general no lineales, el primero nos proporcionará los valores de W n,j en cada paso. Estos valores introducidos en el segundo nos determinarán la solución aproximada yn+1 . Fácilmente puede comprobarse que ambas expresiones son equivalentes a las antes utilizadas para describir el método.

80


7.4. 7.4.1.

Problemas de Valor Inicial y de Contorno Problemas de transporte estacionario 1-dimensionales

Dado el siguiente problema de transporte

−

µu (x) + ηu  (x) + σu(x) = f (x) x

u(a) = u a ,

∈ (a, b)

u(b) = u b

se resuelve mediante el comando: [x,u] = bvp(a,b,N,mu,eta,sigma,f,ua,ub)

en donde mu = µ, eta = η, sigma= σ, f es una función inline y N es el n´ umero de nodos internos (es decir, b−a N = (No. Intervalos) 1 = h 1 en donde h es el paso de discretización.

−

−

En la salida, x = (x1 , . . . , xN +1 ) es el vector que contiene los puntos del mallado y u = (u1 , . . . , u N +1 ) contiene las aproximaciones obtenidas por el método en los puntos del mallado.

7.4.2.

Problemas de difusi´ on evolutiva 1-dimensionales

El siguiente Problema de Valor Inicial y de Contorno (PVIC)

  

∂u ∂ 2 u = µ 2 + f (x, t) ∂t ∂x

∀ (x, t) ∈ (0, L) × (0, T )

u(x, t) = g(x, t)

x = 0,

u(x, 0) = u 0 (x)

x

x = L

∈ (0, L)

se resuelve mediante el comando: [x,u] = heattheta(xspan,tspan,nstep,theta,mu,u0,g,f)

en donde xspan = tspan = nstep = theta = mu = µ u0,g,f

[a,b] [0,T] [Nx,Nt]

θ

- Intervalo espacial. - Intervalo temporal. - Nx , Nt n´ umero de intervalos espaciales y temporales respect. - θ-método de la parte temporal (θ = 0 EE, θ = 0,5 CN, θ = 1 EI) - constante positiva - funciones inline

En la salida Octave muestra directamente la gráfica de la aproximaci´ on de la solución exacta u(x, t) en el instante t = T y almacena en el vector x = (x1 , x2 , . . . , xN +1 ) los puntos del mallado espacial y en T T on buscada en el instante t = T . u = (uT 1 , u2 , . . . , uN +1 ) los valores de la aproximaci´ etodo es estable Observaci´ on. El m´

⇐⇒

µ(1

∆t 1 . − 2θ) (∆x) ≤ 2 2

Cap´ıtulo 8

Otros comandos y funciones u ´ tiles 8.0.1.

Los comandos

max

y

min

Dado un vector cualquiera u los comandos de Octave max y min devuelven los valores máximo y m´ınimo que contiene el vector u. Por ejemplo, si definimos el vector u = (2, 3, 4, 5, 6, 5, 4, 3, 2, 1), entonces max(u) devuelve 6 y min(u) devuelve 1: u = [2,3,4,5,6,5,4,3,2,1] max(u) min(u)

Si además incluimos dos parámetros de salida de la forma [Max,PosMax] = max(u) [Min,PosMin] = min(u)

en las variables Max y Min se almacenan el valor máximo y el m´ınimo de u respectivamente, y en PosMax y PosMin se almacenan las posiciones de u en donde se alcanzan.

8.0.2.

C´ alculo de posiciones. El comando find.

Dado un vector cualquiera u el comando find devuelve la posición (dentro del vector) en la que se encuantra uno o varios valores de u. Es importante tener en cuenta que en Octave los vectores empiezan en la posici´ on 1. Por ejemplo, si tomamos de nuevo el vector u = (2, 3, 4, 5, 6, 5, 4, 3, 2, 1) y escribimos: find(u==6) find(u==4)

81

82


el primer comando devuelve el valor (posición) 5 y el segundo comando devuelve los valores 3 y 7, que son en efecto las posiciones en las que se encuentra el 4 dentro de u. El comando find permite tambi´ en recuperar las posiciones en donde los valores del vector verifican alguna condición >, <, >=, <=, =. Por ejemplo, el comando

∼

find(u>=4)

devuelve las posiciones 3, 4, 5, 6 y 7.

Forma alternativa de calcular posiciones. En la mayor parte de los problemas que trataremos, las

b´ usquedas de posiciones las realizaremos en vectores que representan los mallados t0 < t1 < t2 < .. . < tm < .. . < tN −1 < tN que consideramos para resolver de manera numérica los problemas que se nos plantean. Es este caso, si tomamos como constante el paso de discretización h (constante), podemos recuperar todos los puntos del mallado dentro del vector t en la posición m mediante la fórmula: t(m) = t 0 + (m

− 1)h

la cual nos sirve para calcular la posición m. Por ejemplo, si queremos calcular la aproximación obtenida a partir del método de Euler expl´ıcito en un intervalo [1, 3] en el instante t = 1,6 habiendo tomado un paso de discretización h = 0,1, tenemos que t(m) = 1 + (m

− 1)0,1 = 1,6 =⇒ m = 1,60,1− 1 + 1 = 7

Cap´ıtulo 9

Transformada de Laplace. Funciones definidas a trozos. La funci´ on Heaviside 9.0.1.

Transformada de Laplace. Funciones definidas a trozos. La funci´ on Heaviside

Se puede definir y representar una función derinida a trozos, como f (x) =

 

1 2x , x < 0 1/2 ,x 1 3 1!‘(x + 1) + 2x , x > 1

−

≤− −

de la siguiente forma >> x = linspace(-2, 3, 1000); >> y = (1 2 x). (x > 1); >> plot(x,y) >> plot(x,y,’.’)

∗ ∗

−

Las funciones definidas a trozos tambi´ en pueden escalón: 0 H (x) = 1/2 1

 

− 2 ∗ x). ∗ (x <= −1) + ((x + 1).3 +

definirse usando la funci´ on on Heaviside o funci´ , x < 0 , x = 0 , x > 0

Para utilizarla en Octave usamos la función heaviside.m, la cual podremos combinar para construir todo tipo de funciones definidas a trozos. Por ejemplo, prueba a dibujar las siguientes funciones en el intervalo [0 , 5] usando Octave: f = inline(’heaviside(x-2)’,’x’) f = inline(’heaviside(x-2)-heaviside(x-4)’,’x’) f = inline(´ exp(x).*(heaviside(x-2)-heaviside(x-4))’,’x’)

83

84


9.1.

Transformada de Laplace

Se puede calcular la transformada de Laplace con Octave Symbolic. Con la orden laplace(f(t),t,s) o bien laplace(f) o laplace(f,s) se puede calcular la transformada de Laplace. La transformada inversa se puede calcular con ilaplace. >> load symbolic >> syms t >> f (t) = t 2 + t3 sin(t)

∗

>> laplace(f) >> g = t heaviside(t) + t2 heaviside(t

∗

∗

− 1)

>> laplace(g) >> ilaplace((s2 + s + 1)/(s3 + 4 s2 + 13 s))

∗

>> G(s) = 1/s >> ilaplace(G)

∗

− exp(−2 ∗ s) ∗ s/(s3 − 4 ∗ s2 + 4 ∗ s)

Cap´ıtulo 10

An´ alisis is is Vector ct oria iall 10.0.1 10.0.1..

Repres Represen entac taci´ i´ on de Campos Vectoriales on

De manera similar a cómo omo se representaron los campo de direcciones de una ecuación on diferencial 7.1, −y x se puede representar un campo de vectores, por ejemplo F = x2 +y2 , x2 +y2



>> close >> close all; clear all; clc; >> x =

−1 : 0,1 : 1;

>> y = x = x;; >> [ >> [X, X, Y ] Y ] = meshgrid( meshgrid(x, y ); >> F X = X =

−Y./sqrt( Y./sqrt(X.2 + Y. 2 );

>> F Y = X./sqrt( X./sqrt(X.2 + Y. 2 ); >> quiver( quiver(X , Y , F X , F Y ); ); >> gr >> grid id on;

85



86


Del mismo modo podr´ıamos ıamos representar un campo en R3 , >> x = 0 : 0, 0,2 : 1; y = x = x;; z = 0 : 0, 0,2 : 2; >> [ >> [X X , Y , Z ] = meshgrid( meshgrid(x,y,z); x,y,z); >> U =

Y./sqrt(X.2 + Y. 2 ); −Y./sqrt(

>> V = X./sqrt( X./sqrt(X.2 + Y. 2 ); >> W = Z. 2 /2 + 1; >> h = h = quiver quiver3( 3(X,Y,Z,U,V,W X,Y,Z,U,V,W ); ); >> grid >> grid on; >> set( set(h, linewidth  , 2);



−y x2 y 2

+

x , x2 + ,1+ y2

z2

2



87




Del mismo mi smo modo mo do podr´ p odr´ıamos ıamos representar 0, 0, 9

10.0.2 10.0.2..

− x2 − y2



Repres Represen entac taci´ i´ on de curvas y superficies on

Vemos diversas representaciones que nos pueden ayudar en los ejercicios de Análisis Vectorial. Para representar una curva en polares usamos la orden polar(x,y) >> x = linspace = linspace(( pi, pi, 100);

−

>> y = cos = cos(2 (2.. x);

∗

88


>> polar(x, y)

Una curva en paramétricas [tcos(t),tsin(t)] se podr´ıa representar directamente con plot >> t = [ 4 pi : pi/50 : 4 pi];

− ∗ ∗ >> plot(t. ∗ cos(t), t. ∗ sin(t))

O en tres dimensiones [tcos(t),tsin(t), 2 t] >> plot3(t. cos(t), t. sin(t), 2 t)

∗

∗

∗

>>grid(’ on’);

O bien, en dos dimensiones en polares, r = 3 >> theta = linspace(0, 10 pi, 300);

∗

>> ro = 3

− 5. ∗ theta; >> x = ro. ∗ cos(theta); y = ro. ∗ sin(theta);

− 5 ∗ t, pasándolas a paramétricas

89


>> plot(x, y) >> axis equal

Con la orden mesh se puede representar una superficie, sus curvas de nivel con contour, o ambas cosas en una sóla gráfica con meshc, para z = 9 + x2 + y2 ,



Creamos una malla >> x = [ 2 : 0,1 : 2]; y = x;

−

>> [X, Y ] = meshgrid(x, y) Definimos la componente Z >> Z 3 = sqrt(9 + X.2 Y representamos >> mesh(X , Y , Z 3 ) >> contour(X , Y , Z 3 ) >> meshc(X , Y , Z 3 )

− Y. 2)

90


Con la orden surf se pueden representar superficies en paramétricas. Para representar un toro definimos >> R = 2; r = 1; >> t1 = 0 : pi/25 : 2 pi; p = t1;

∗

>> [T, P ] = meshgrid(t1, p); >> X = R cos(P ) + r cos(T ). cos(P );

∗ ∗ ∗ >> Y = R ∗ sin(P ) + r ∗ cos(T ). ∗ sin(P ); >> Z = r ∗ sin(T ); >> surf (X , Y , Z ) ;

Hay órdenes que representan directamente una superficie conocida, se pueden ver en https://www.gnu.org/software/octave/doc/interpreter/Three 002ddimensional-Geometric-Shapes.html por ejemplo para una esfera basta: >> sphere,axis square

91


10.0.3.

C´ alculos de gradiente, divergencia, rotaciona. Representaciones.

Octave Symbolic nos permite diversos cálculos de análisis vectorial como el gradiente, la divergencia o el rotacional con las órdenes gradient, divergence o curl respectivamente. >> pkg load symbolic >> syms x y z Para la divergencia y rotacional >> F = [x2

− x, y2 − z, x2 + y2 + z2]

>> divergence(F ) ans = (sym) 2*x + 2*y + 2*z - 1 >> curl(F ) [2 y + 1]

∗ [ −2 ∗ x [

0

] ]

Por ejemplo para el gradiente de la función f = x 2 >> f = x 2

− y3 + 3 ∗ x ∗ z

>> g = gradient(f ) g = (sym 3x1 matrix) [2*x + 3*z] [ [

−3 ∗ y2 ] 3∗x ]

Y para su valor en un punto concreto >> g 1 = subs(g,x,y,z, 1, 2, 3)

−

− y3 + 3 x z,


92

g1 = (sym 3x1 matrix) [-7 ] [-12] [3] También lo podr´ıamos haber calculado directamente sin usar la orden, con >> fx = diff (f, x) >> fy = diff (f, y) >> fz = diff (f, z) >> gx = subs(fx,x,y,z, 1, 2, 3)

− >> gy = subs(fy,x,y,z, 1, 2, −3) >> gz = subs(fz,x,y,z, 1, 2, −3) Calculemos el gradiente de f = −3 x2 − 4 y 2 , y representamos la superficie y las curvas de nivel, la

particular que pasa por el punto (1, 2) y los vectores gradiente y tangente >> f =

−3 ∗ x2 − 4 ∗ y2;

>> g = gradient(f ); >> gs = subs(g,x,y, 1, 2) >> mg = sqrt(gs(1)2 + gs(2)2 ) >> c1 = gs(1)/mg >> c2 = gs(2)/mg >> vpa(c1) ans = (sym) -0.3511 >> vpa(c2) ans = (sym) -0.9363 >> t1 =

−3 : 0,2 : 3;

>> [X, Y ] = meshgrid(t1, t1); >> Z =

−3 ∗ X.2 − 4 ∗ Y. 2;

>> surf (X , Y , Z ) y para las curvas de nivel >> cn = contour(X,Y,Z, 10); >> clabel(cn) >> axisequal

Octave Practicas

Recommend Documents