Física Computacional 23/24

Tema 5: Soluciones de ecuaciones

Índice
1. Soluciones de ecuaciones lineales 2
1.1. Eliminación gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1. Ejemplo 5.1: la eliminación gaussiana . . . . . . . . . . . . . . . . . . . . . . . 4
1.2. Método del pivote . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.1. Ejemplo 5.2: incluyendo el pivote . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3. Factorización LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1. Forma matricial de las operaciones elementales sobre matrices . . . . . . . . . 9
1.3.2. Matrices L y U . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3.3. La matriz L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.4. Ejemplo 5.5: aplicando la descomposición LU . . . . . . . . . . . . . . . . . . 12
1.3.5. Ejercicio 5.1: la descomposición LU usando el pivote . . . . . . . . . . . . . . 15
1.3.6. Ejemplo 5.6: Contando el número de operaciones de la descomposición LU . 15
1.4. La inversa de una matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.4.1. Ejemplo 5.7: la inversa de una matriz . . . . . . . . . . . . . . . . . . . . . . . 16
1.5. Matrices tridiagonales y matrices banda . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.5.1. Ejemplo 5.8: la vibración en un sistema unidimensional . . . . . . . . . . . . . 17
1.6. Métodos iterativos para ecuaciones lineales . . . . . . . . . . . . . . . . . . . . . . . . 24
1.6.1. El método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.6.2. Método de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.6.3. Método de Sobrerelajación sucesiva . . . . . . . . . . . . . . . . . . . . . . . . 27

2. Autovalores y autovectores 29
2.1. Descomposición QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1.1. Ejemplo 5.12: estudiando la descomposición QR . . . . . . . . . . . . . . . . . 30
2.2. Descomposiciones QR sucesivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3. Algoritmo QR para el cálculo de autovectores y autovalores . . . . . . . . . . . . . . 33
2.3.1. Ejemplo 5.13: calculando autovalores y autovectores . . . . . . . . . . . . . . 33
2.3.2. Autovalores y autovectores en python . . . . . . . . . . . . . . . . . . . . . . . 35
2.3.3. Ejercicio 5.2: Pozo de potencial asimétrico en Mecánica Cuántica . . . . . . . 36

3. Ecuaciones no lineales 38
3.1. Método del punto fijo o relajación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.1. Ejemplo 5.14: aplicando el método del punto fijo . . . . . . . . . . . . . . . . . 39
3.1.2. No convergencia del método del punto fijo . . . . . . . . . . . . . . . . . . . . 40
3.1.3. Limitaciones del método del punto fijo . . . . . . . . . . . . . . . . . . . . . . 42
3.1.4. Convergencia y estimación del error del método del punto fijo/relajación . . 43
3.1.5. Ejercicio 5.3: propagación de epidemias . . . . . . . . . . . . . . . . . . . . . . 45
3.1.6. El método del punto fijo sobrerelajado . . . . . . . . . . . . . . . . . . . . . . . 45
3.1.7. Método del punto fijo o relajación con más de una variable . . . . . . . . . . . 46
3.2. Método de la bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

1
 3.2.1. Algoritmo para el método de la bisección . . . . . . . . . . . . . . . . . . . . . 48
3.2.2. Convergencia del método de la bisección . . . . . . . . . . . . . . . . . . . . . 49
3.2.3. Limitaciones del método de la bisección . . . . . . . . . . . . . . . . . . . . . . 49
3.2.4. Ejemplo 5.17: visualizando el método de la bisección . . . . . . . . . . . . . . 51
3.2.5. Ejercicio 5.6: la constante de desplazamiento de Wien . . . . . . . . . . . . . . 54
3.3. Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.3.1. Error y convergencia del método de Newton-Raphson . . . . . . . . . . . . . 56
3.3.2. Desventajas del método de Newton-Raphson. . . . . . . . . . . . . . . . . . . 58
3.3.3. Ejemplo 5.18: visualizando el método de Newton-Raphson . . . . . . . . . . . 58
3.3.4. Ejercicio 5.7: función arcotangente hiperbólica . . . . . . . . . . . . . . . . . . 60
3.3.5. Ejercicio 5.8: las raíces de un polinomio . . . . . . . . . . . . . . . . . . . . . . 60
3.3.6. Método de la secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.3.7. Método de Newton-Raphson para más de una variable . . . . . . . . . . . . . 61

4. Máximos y Mínimos de Funciones: minimización 68

4.1. Estrategia general . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.2. Método de la razón áurea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2.1. Algoritmo del método de la razón áurea . . . . . . . . . . . . . . . . . . . . . . 72
4.2.2. Limitaciones del método de la razón áurea . . . . . . . . . . . . . . . . . . . . 72
4.2.3. Ejemplo 5.20: el potencial de Buckingham . . . . . . . . . . . . . . . . . . . . . 72
4.2.4. Ejercicio 5.10: la temperatura de una bombilla . . . . . . . . . . . . . . . . . . 75
4.3. El método de Gauss-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.4. El método del gradiente descendente . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.5. El método del gradiente descendente en varias dimensiones . . . . . . . . . . . . . . 77
4.5.1. Método de máxima pendiente . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.5.2. Método del gradiente conjugado . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5. Ejercicio adicionales 81
5.1. Ejercicio 5.11: minimizando el potencial de ionización (Examen de Enero de 2023) . . 81
5.2. Ejercicio 5.12: funciones complejas con polos (Examen de Junio de 2023) . . . . . . . 81
5.3. Ejercicio 5.13: cuadrando posiciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.4. Ejercicio 5.14: teoría del campo medio para el ferromagnetismo . . . . . . . . . . . . 84
5.5. Ejercicio 5.15: potencial cuadrado cuántico . . . . . . . . . . . . . . . . . . . . . . . . 84
5.6. Ejercicio 5.16: el punto de Lagrange Tierra-Luna . . . . . . . . . . . . . . . . . . . . . 86
5.7. Ejercicio 5.17: ajustando el fondo cósmico de microondas . . . . . . . . . . . . . . . . 87

1. Soluciones de ecuaciones lineales

Queremos resolver numéricamente un sistema de n ecuaciones lineales con n incógnitas de
la forma:

2
 a11 x1 + a12 x2 + · · · + a1n xn = b1 ,
a21 x1 + a22 x2 + · · · + a2n xn = b2 ,
.. .. .. ..
. . . .
an1 x1 + an2 x2 + · · · + ann xn = bn ,
(1)

o en formulación matricial:

A · X = B,

con

     
a11 a12 ··· a1n x1 b1
 a21 a22 ··· a2n   x2   b2 
A= , X= y B= .
     
.. .. .. .. ..  ..
 . . . .   .   . 
an1 an2 · · · ann xn bn

¿ Cuál es el método más efectivo para resolver este sistema numéricamente?

A primera vista puede parecer que el método más sencillo es calcular la inversa de A, es
decir:

X = A−1 · B,

sin embargo no es el caso.

Calcular la inversa de una matriz es en general más complicado que resolver el sistema
inicial de ecuaciones.
Además, la inversa involucra cocientes de diferencias de números que aumenta el error de
redondeo.
Hay métodos más fáciles y rápidos.

1.1. Eliminación gaussiana

La idea del método es utilizar operaciones elementales de matrices que permitan transfor-
mar nuestro sistema de ecuaciones inicial en uno donde la matriz A sea triangular.
Para ello construimos en primer lugar la matriz ampliada:

 
a11 a12 ··· a1n | b1
 a21 a22 ··· a2n | b2 
( A| B) = 
 
.. .. .. .. .. 
 . . . . | . 
an1 an2 · · · ann | bn

3
 Operaciones elementales sobre la matriz extendida:
1. Multiplicar una fila por un escalar (no nulo).
2. Sumar a una fila el múltiplo de otra.
3. Intercambiar las posición de dos filas.
El método de Gauss combina las operaciones 1 y 2 para transformar:

   
a11 a12 a31 ··· a1n | b1 1 ā12 ā31 ··· ā1n | b̄1

 a21 a22 a32 ··· a2n | b2 

 0 1 ā32 ··· ā2n | b̄2 

( A| B) = 
 a31 a32 a33 ··· a3n | b3  ⇒

 0 0 1 ··· ā3n | b̄3 

 .. .. .. .. .. ..   .. .. .. .. .. .. 
 . . . . . | .   . . . . . | . 
an1 an2 an3 · · · ann | bn 0 0 0 ··· 1 | b̄n

donde el segundo sistema de ecuaciones se puede resolver por sustitución inversa.

¿Cómo podemos obtener la matriz triangular?
Dividimos la primera fila por el primer elemento a11 .
Le restamos a la fila j > 1 la primera fila multiplicada por su primer elemento a j1 .
Dividimos la segunda fila por su segundo elemento a22 .
Le restamos a la fila j > 2 la segunda fila multiplicada por su primer elemento a j2 .
Continuamos de la misma fórmula hasta llegar al último elemento.
¿Cómo realizar la sustitución hacía atrás?
El valor de xn = b̄n
El valor de xn−1 = b̄n−1 − ān−1,n xn
Continuamos restándole a las filas anteriores el valor de todas las incógnitas que acabamos
de resolver multiplicadas por su correspondiente coeficiente en la fila.
Algoritmo para obtener la eliminación gaussiana:
1. Dividimos la fila m por su elemento m-esimo.
2. A cada fila i > m le restamos la fila m multiplicada por su elemento ami .
3. El valor de la variable xm = bm − ∑iN=m+1 am,i xi

1.1.1. Ejemplo 5.1: la eliminación gaussiana

Usar el método de la eliminación de Gauss para transformar el sistema de ecuaciones:

2 w + x + 4 y + z = − 4,
3 w + 4 x − y − z =3,
w − 4 x + y + 5 z =9,
2 w − 2 x + y + 3 z =7.
(2)

4
[7]: from numpy import array,zeros,column_stack,dot

def gaussian_elim(A,B):
# Calculamos el número de filas
N=len(B)

# Creamos nuestra matriz ampliada

AB=column_stack((A,B))

# Eliminación gaussiana

for m in range(N):

# 1. dividimos por el elemento

div=AB[m,m]
AB[m,:]/=div

# 2. Sustraemos las filas siguientes

for i in range(m+1,N):
mult=AB[i,m]
AB[i,:]-=mult*AB[m,:]

# 3. Sustitución inversa
x=zeros(N,float)

for m in range(N-1,-1,-1):
x[m]=AB[m,N]
for i in range(m+1,N):
x[m]-=AB[m,i]*x[i]
return x,AB

[8]: # Definimos nuestra matriz ampliada

A=array([[2,1,4,1],[3,4,-1,-1],[1,-4,1,5],[2,-2,1,3]],float)
B=array([-4,3,9,7],float)

AB=column_stack((A,B))
print("Nuestra matriz extendida inicial")
print(AB)
print("")

x,AB2=gaussian_elim(A,B)

print("Nuestra matriz extendida triangular")

print(AB2)
print("")

5
print("La solución del sistema de ecuaciones es:")
print(x)
print("")

print("Comprobamos que es solución")

print(dot(A,x)-B)

Nuestra matriz extendida inicial

[[ 2. 1. 4. 1. -4.]
[ 3. 4. -1. -1. 3.]
[ 1. -4. 1. 5. 9.]
[ 2. -2. 1. 3. 7.]]

Nuestra matriz extendida triangular

[[ 1. 0.5 2. 0.5 -2. ]
[ 0. 1. -2.8 -1. 3.6]
[-0. -0. 1. -0. -2. ]
[-0. -0. -0. 1. 1. ]]

La solución del sistema de ecuaciones es:

[ 2. -1. -2. 1.]

Comprobamos que es solución

[ 4.44089210e-16 1.77635684e-15 -1.77635684e-15 2.66453526e-15]

1.2. Método del pivote

La eliminación gaussiana requiere dividir en cada paso por el elemento aii , lo cual sólo es
posible si el elemento es diferente de cero.
El método del pivote combina la eliminación gaussiana con la operación elemental 3 de las
matrices:
el intercambio de una fila por otro.
Esto garantiza que en ningún paso un elemento de la diagonal sea 0.
No hay una única forma de hacerlo.
La regla general es intercambiar en cada paso la fila m por la fila j > m tal que | a jm | tenga el
valor máximo de todos los akm con m ≤ k ≤ n.

1.2.1. Ejemplo 5.2: incluyendo el pivote

Modificar el algoritmo de la eliminación de gauss para incluir el método del pivote.
1. Comprobar que el resultado del ejercicio 5.1 no cambia.
2. Aplicarlo al sistema:

6
 x + 4 y + z = − 4,
3 w + 4 x − y − z =3,
w − 4 x + y + 5 z =9,
2 w − 2 x + y + 3 z =7.
(3)

y comprobar que el sistema funciona.

[97]: from numpy import array,zeros

def gaussian_pivot(A,B):

# Calculamos el número de filas

N=len(B)

# Creamos nuestra matriz ampliada

AB=column_stack((A,B))

# Eliminación gaussiana con método del pivote.

for m in range(N):

#1. Escogemos la fila i>=m cuyo elemento a_{im} este más alejado de 0.
pivot= m
largest= abs(AB[m,m])

for i in range(m+1,N):
if abs(AB[i,m]>largest):
largest=AB[i,m]
pivot=i

# 2. Intercambiamos la fila pivot por la fila m.

if pivot!=m:
for i in range(N+1):
AB[m,i],AB[pivot,i]=AB[pivot,i],AB[m,i]

# 3. dividimos por el elemento

div=AB[m,m]
AB[m,:]/=div

# 4. Sustraemos las filas siguientes

for i in range(m+1,N):
mult=AB[i,m]
AB[i,:]-=mult*AB[m,:]

7
 # 5. Sustitución inversa
x=zeros(N,float)

for m in range(N-1,-1,-1):
x[m]=AB[m,N]
for i in range(m+1,N):
x[m]-=AB[m,i]*x[i]

return x,AB

[98]: # Definimos nuestra matriz ampliada

A=array([[0,1,4,1],[3,4,-1,-1],[1,-4,1,5],[2,-2,1,3]],float)
B=array([-4,3,9,7],float)

AB=column_stack((A,B))

print("Nuestra matriz extendida inicial")

print(AB)
print("")

x,AB2 = gaussian_pivot(A,B)

print("Nuestra matriz extendida triangular")

print(AB2)
print("")

print("La solución del sistema de ecuaciones es:")

print(x)
print("")

print("Comprobamos que es solución")

print(dot(A,x)-B)

Nuestra matriz extendida inicial

[[ 0. 1. 4. 1. -4.]
[ 3. 4. -1. -1. 3.]
[ 1. -4. 1. 5. 9.]
[ 2. -2. 1. 3. 7.]]

Nuestra matriz extendida triangular

[[ 1. 1.33333333 -0.33333333 -0.33333333 1. ]
[ 0. 1. 4. 1. -4. ]
[ 0. 0. 1. 0.47058824 -0.58823529]
[-0. -0. -0. 1. 1.38095238]]

La solución del sistema de ecuaciones es:

8
[ 1.61904762 -0.42857143 -1.23809524 1.38095238]

Comprobamos que es solución

[ 0.00000000e+00 0.00000000e+00 -3.55271368e-15 0.00000000e+00]

1.3. Factorización LU
La eliminación gaussiana con pivote es un método efectivo y rápido.
Sin embargo, en muchas ocasiones en física queremos resolver el sistema:

A · X = V,

para distintos valores de V.

Como el efecto de la eliminación gaussiana sobre A siempre es el mismo, la pregunta es si
uno puede almacenar de alguna forma las operaciones necesarias para convertir A en una
matriz triangular.

1.3.1. Forma matricial de las operaciones elementales sobre matrices

El punto de partida es la matriz identidad I (n).
Toda operación elemental se puede expresar como una operación sobre I (n) que luego ac-
tuará en A.
1. Multiplicar una fila m por un escalar α.
Por ejemplo, multiplicar la fila segunda por un α.

 
1 0 ··· 0
 0 α ··· 0 
I1 (n, 2) =  .
 
.. .. . .
 . . . 0 
0 0 ··· 1
2. Sumar a una fila k la fila m multiplicada por un escalar.
Por ejemplo, restar a la segunda fila la primera multiplicada por α.

 
1 0 ··· 0
 −α 1 · · · 0 
I2 (n, 2, 1) =  .
 
.. .. . .
 . . . 0 
0 0 ··· 1
3. Intercambiar la fila k por la fila m.
Por ejemplo, intercambiar la fila 2 por la fila 1.

 
0 1 ··· 0
 1 0 ··· 0 
I3 (n, 1, 2) =  .
 
.. .. . .
 . . . 0 
0 0 ··· 1

9
 Definidas esas operaciones. Se puede construir para cada uno de los pasos de la eliminación
gaussiana la matriz que nos permite transforma A en una matriz diagonal.

Ejemplo 5.3: transformaciones de matrices Para una matriz A de dimensión 4x4

 
a11 a12 a13 a14
 a21 a22 a23 a24 
A=
 a31

a32 a33 a34 
a41 a42 a43 a44

definir las matrices asociadas a las 4 operaciones necesarias para convertirla en una matriz trian-
gular superior.
1. Dividir la primera fila por su primer elemento y restarsela al resto de filas al multipicarla
por α j1 con j>1.

 
1 0 0 0
L1 =
1  − a21 a11 0 0 ,
a11  − a31 0 a11 0 
− a41 0 0 a11

que permite

 
1 b12 b13 b14
 0 b22 b23 b24 
L1 · A =  =B
 0 b32 b33 b34 
0 b42 b43 b44
2. Dividir la segunda fila por su segundo elemento y restarsela al resto de filas al multipicarla
por β j2 con j>2.

 
b22 0 0 0
1  0 1 0 0 
L2 = ,
b22  0 −b32 b22 0 
0 −b42 0 b22

que permite

 
1 c12 c13 c14
 0 1 c23 c24 
L2 · B =  =C
 0 0 c33 c34 
0 0 c43 c44

3. Dividir la tercera fila por su tercer elemento y restarsela al la cuarta fila multiplicada por c43 .

10
  
c33 0 0 0
1  0 c33 0 0 
L3 = ,
c33  0 0 1 0 
0 0 −c43 c33

que permite

 
1 d12 d13 d14
 0 1 d23 d24 
L3 · C =  =D
 0 0 1 d34 
0 0 0 d44
4. Dividir la cuarta fila por su cuarto elemento.

 
d44 0 0 0
1  0 d 44 0 0 
L4 =  ,
d44  0 0 d44 0 
0 0 0 1

que nos permite convertir la matriz D en una matriz triangular superior.

1.3.2. Matrices L y U
Por tanto, para una matriz con n filas, necesitamos definir n matrices Li con i = 1, · · · , n.
Una vez obtenidas el problema reside en operar con ellas sobre la matriz extendida y realizar
la sustitución hacía atrás:

Ln · · · L2 · L1 · A · X = Ln · · · L2 · L1 · V,

Puesto que el producto Ln · · · L2 · L1 · A siempre es el mismo, no hace falta volverlo a calcular

si se cambia el valor de V.
En la práctica sin embargo se intenta minimizar las operaciones sobre V,
de forma que sus valores sólo se utilizan para la sustitución hacía atrás.
Ello se consigue definiendo la matrices:

L = L1−1 L2−1 · · · L−
n
1
y U = Ln · · · L2 · L1 · A,

que satisfacen que:

L · U = A, ⇒ L·U·X = V

U, tal y como hemos visto, es una matriz triangular superior.

En principio el cálculo de L requiere calcular la inversa de cada una de las matrices elemen-
tales y operar con ellas.

11
 Ejemplo 5.4: cálculo de L Calcular el valor de L para la matriz del ejercicio 5.3
1. Calculamos la inversa de cada una de las matrices elementales.

       
a11 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0
 a21 1 0 0   0 b22 0 0   0 1 0 0   0 1 0 0 
L1−1 =
 a31
, L2−1 =
 0 b32
, L3−1 = , L4−1 = 
0 1 0  1 0   0 0 c33 0   0 0 1 0 
a41 0 0 1 0 b42 0 1 0 0 c43 1 0 0 0 d44
2. Operamos con todas ellas:

 
a11 0 0 0
 a21 b22 0 0 
L = L1−1 L2−1 L3−1 L4−1 =
 a31 b32 c33 0  .


a41 b42 c43 d44

1.3.3. La matriz L
Pero no es necesario:
• la matriz L es la matriz triangular cuyas columnas se generan en cada paso al convertir
la matriz A en una matriz triangular superior, es decir, cuando se calcula U.
La matriz L es una matriz diagonal inferior mientras que la matriz U es una matriz diagonal
superior.
Por tanto, cada una de ellas puede resolverse por medio de la sustitución, hacía atrás para U
y hacía adelante para L.
Específicamente, queremos resolver:

L · U · X = V,

que se puede descomponer en:

U·X =Y

y el problema:

L · Y = V.

1.3.4. Ejemplo 5.5: aplicando la descomposición LU

Aplicar la descomposición LU a la matriz del ejercicio 5.1.

[13]: def LU(A):

from numpy import zeros,copy

12
 N=len(A)

# Inicialiazmos nuestras dos matrices.

L=zeros([N,N],float) # L tiene que ser una matriz triangular inferior.
U=copy(A) # U será la matriz A convertida en triangular␣
,→superior.

# Creamos la matrix L
for m in range(N):
L[m:N,m] = U[m:N,m] # para cada iteracción (para cada fila),
# la columna m de L queda fijada por el valor que␣
,→tiene U.

# Convetimos ahora U en una matrix triangular superior.

# Para ello usamos la eliminación guassiana.

# 1. Dividimos la fila m por el elemento m,m

div=U[m,m]
U[m,:]/=div

# 2. Sustraemos la fila m a las filas i>m multiplicadas por el elemento␣

,→ i,m
for i in range(m+1,N):
mult=U[i,m]
U[i,:]-= mult*U[m,:]

return L,U

[14]: # Definimos nuestra matriz

from numpy import dot

A=array([[2,1,4,1],[3,4,-1,-1],[1,-4,1,5],[2,-2,1,3]],float)

# Generamos las matrices L y U y comprobamos que son como deberían ser.

L,U=LU(A)
print(L)
print()
print(U)
print()

# Comprabamos que L U es efectivamente A.

print(dot(L,U))
print()
print(A)

[[ 2. 0. 0. 0. ]
[ 3. 2.5 0. 0. ]

13
 [ 1. -4.5 -13.6 0. ]
[ 2. -3. -11.4 -1. ]]

[[ 1. 0.5 2. 0.5]
[ 0. 1. -2.8 -1. ]
[-0. -0. 1. -0. ]
[-0. -0. -0. 1. ]]

[[ 2. 1. 4. 1.]
[ 3. 4. -1. -1.]
[ 1. -4. 1. 5.]
[ 2. -2. 1. 3.]]

[[ 2. 1. 4. 1.]
[ 3. 4. -1. -1.]
[ 1. -4. 1. 5.]
[ 2. -2. 1. 3.]]

[15]: # programa que usa las matrices LU

# para resolver un sistema de ecuaciones

def LU_sol(L,U,V):

from numpy import empty

N=len(V)

# 1. Sustitución hacía adelante. Operando con L.

y=empty(N,float)

for m in range(N):
y[m]=V[m]
for i in range(m):
y[m]-= L[m,i]*y[i]
y[m]/=L[m,m]

# 2. Sustitución hacía atrás. Operando con U.

x=empty(N,float)

for m in range(N-1,-1,-1):
x[m]=y[m]
for i in range(m+1,N):
x[m]-=U[m,i]*x[i]

return(x)

[16]: V=array([-4,3,9,7],float)
print(LU_sol(L,U,V))

14
 [ 2. -1. -2. 1.]

1.3.5. Ejercicio 5.1: la descomposición LU usando el pivote

Aplicar la descomposición LU incluyendo el pivote a la matriz del ejemplo 5.1.
(Ayuda: al no estar usando en este caso la matriz extendida es necesario almacenar la información
de que columnas se han ido intercambiando en cada paso).

[ ]:

Como era de esperar python ya tiene incorporada una rutina para resolver un sistema de
ecuaciones.
La descomposición LU se usa en la función solve del módulo linalg del paquete numpy.

[20]: from [Link] import solve

A=array([[0,1,4,1],[3,4,-1,-1],[1,-4,1,5],[2,-2,1,3]],float)
V=array([-4,3,9,7],float)
print(solve(A,V))
print()

[ 1.61904762 -0.42857143 -1.23809524 1.38095238]

1.3.6. Ejemplo 5.6: Contando el número de operaciones de la descomposición LU

Calcular el número de operaciones necesarias para realizar la descomposición LU.
1. Para la primera fila, dividimos cada columna por el primer elemento, y se la restamos a las
demás filas multiplicadas por el primer elemento de cada una de ellas:

N + ( N − 1) · 2N = N (2N − 1).
2. Para la segunda fila, repetimos el procedimiento pero empezando en la segunda columna:

N − 1 + ( N − 2) · 2( N − 1) = ( N − 1) · (2N − 3).
3. Para la fila i-esima:

N − i + 1 + ( N − i ) · 2( N − i + 1) = ( N − i + 1) · (2 ( N − i ) + 1)) .
4. Sumando todos los términos:

N N
2 N 2
#operaciones LU = ∑ ( N − i + 1) (2 ( N − i) + 1)) = ∑ i (2i − 1) = 3 N 3 + 6
(3N − 1) ∼ N 3 .
3
i =1 i =1

La descomposición LU es un cálculo computacionalmente pesado.

15
 1.4. La inversa de una matriz
La matriz de inversa se calcula como:

1
A −1 = [C ] T
det( A)

con Cij el cofactor de Aij .

Sin embargo su aplicación directa es computacionalmente pesada.
Además, conlleva errores de redondeo que pueden ser no despreciables.
En aquellos casos donde sea necesario calcularla el procedimiento es resolver el sistema:

A · A −1 = I

aplicando la eliminación gaussiana o la descomposición LU a cada una de las columnas de

A−1 y de I.

1.4.1. Ejemplo 5.7: la inversa de una matriz

Encontrar la inversa de la matriz del ejercicio 5.1 usando el método que se prefiera.

[21]: A=array([[2,1,4,1],[3,4,-1,-1],[1,-4,1,5],[2,-2,1,3]],float)
L,U=LU(A)
N=len(A)
Ainv=zeros([N,N],float)

for i in range(N):
I=zeros(N,float)
I[i]=1
Ainv[:,i]=LU_sol(L,U,I)

[22]: from numpy import set_printoptions

set_printoptions(precision=2)
set_printoptions(suppress=True)

print(dot(A,Ainv))

[[ 1. 0. 0. 0.]
[ 0. 1. 0. 0.]
[-0. -0. 1. 0.]
[-0. 0. 0. 1.]]
Final de la clase del 27/10/21

16
1.5. Matrices tridiagonales y matrices banda
En física es muy común encontrar sistemas de ecuaciones:

A · X = B,

con A siendo una matriz tridiagonal, es decir:

 
a11 a12
 a21 a22 a23 
 
 a32 a33 a34 
A=
 
.. .. .. 

 . . . 

 an−2n−1 an−1n−1 an−1n 
ann−1 ann

donde todos los elementos que no aparecen son cero.

Por ejemplo en el algoritmo para definir los splines cúbicos o en general en todos aquellos
sistemas donde sólo hay interacciones entre vecinos.
Los problemas tridiagonales son perfectos para aplicar la eliminación gaussiana, y su algo-
ritmo pueden simplificarse con respecto al general.
La idea es que para un paso dado de la eliminación gaussiana (para una fila dada),
no es necesario recorrer todas las demás, sólo la siguiente.
De la misma forma en la sustitución hacía atrás, la ecuación para la variable xn sólo involu-
crará a xn+1 .
Las matrices banda son aquellas en las que:

aij = 0 si j < i − k1 || j > i + k2 con k1 , k2 > 0.

El algoritmo para resolverlas será igual al de las matrices tridiagonales pero recorriendo en
cada paso las k1 filas siguientes.

1.5.1. Ejemplo 5.8: la vibración en un sistema unidimensional

Sea un sistema de N partículas idénticas en una dimensión unidas por muelles. Una vez des-
plazado del equilibrio, cada una de las partículas vibrará con respecto a las demás debido a los
muelles.
(Este modelo describe de forma bastante precisa la vibración de los átomos en algunos sólidos.)
La dinámica de este sistema viene dada por la ecuaciones de Newton. Llamando xi al desplaza-
miento de la partícula i con respecto a su posición de reposo, la segunda ley impone:

17
 m ẍ1 = k ( x2 − x1 ) + F1 ,
m ẍi = k ( xi+1 − xi ) + k ( xi−1 − xi ) + Fi , con 1 < i < N,
m ẍ N = k ( x N −1 − x N ) + FN ,
(4)

siendo m la masa de las partículas, k la constante del muelle y Fi cualquier fuerza externa aplicada
sobre la partícula i.

[71]: from [Link] import display,Image

display(Image(filename='[Link]'))

Asumamos ahora que accionamos el sistema al aplicar una fuerza armónica (sinusoidal) sobre la
primera partícula de la forma:

F1 = Cei ω t ,

siendo C una constante y ω la frecuencia de oscilación (empleamos la forma compleja para sim-
plificar los cálculos entendiendo que al final tomaremos la parte real).
La consecuencia de aplicar esta fuerza será que el sistema oscilará con la frecuencia w, de forma
que la posición de cada partícula vendrá dada por:

x i = a i ei ω t ,

siendo ai la amplitud del movimiento de la partícula i.

Introduciendo la solución en las ecuaciones de Newton obtenemos:

−m ω 2 a1 = k ( a2 − a1 ) + C,
− m ω 2 a i = k ( a i +1 − a i ) + k ( a i −1 − a i ), con 1 < i < N,
2
− m ω a N = k ( a N −1 − a N ),
(5)

que, simplificando, pueden escribirse como:

18
 (k − α) a1 − k a2 =C,
−k ai−1 + α ai − k xi+1 =0, con 1 < i < N,
−k a N −1 + (α − k) a N =0,
(6)

con α = 2k − m ω 2 que en notación matricial toma la forma:

(α − k) −k
     
a1 C
 −k α −k   a2   0 
     
 −k α −k   a3   0 
· = ,
     
 . .. ... ... .. ..


 
  .  
  . 

 −k α −k   a N −1   0 
−k (α − k) aN 0
1. Para un sistema con N=30 partículas de masa m=1 Kg, ligadas por muelles con constantes
de elasticidad k = 5 N/m, y que sufren una fuerza armónica de módulo 1 N con frecuencia
angular ω = 2 rad/s calcular la posición de cada una de las partículas con respecto a su
posición de equilibrio.
Para ello escribir un programa que optimice la eliminación gaussiana para el caso de tener
una matriz tridiagonal.
2. Generalizar el programa para el caso de tener una matriz banda con up elementos no nulos
por encima de la diagonal y down por debajo.
3. Representar la evolución del sistema por medio de una animación.
[Ayuda: para mayor claridad colocar las partículas en su posición de reposo separadas dos
unidades en una recta].
[38]: # programa que optimiza la eliminación gussiana para matrices tridiagonales.

def tridiagonal_gaussian(A,v):
N=len(v)

for i in range(N-1):
# 1. Dividimos la fila por su primer elemento
# sólo tenemos que hacerlo con el elemento a la derecha
# no hace falta que mandemos el elmento A[i,i] a 1-
# Es suficiente con tenerlo en cuenta al hacer la eliminación hacía␣
,→atrás.

A[i,i+1]/=A[i,i]
v[i]/=A[i,i]

# 2. Le sustraemos a la fila i+1 la fila i multiplicada por A[i+1,i]

19
 A[i+1,i+1]-=A[i+1,i]*A[i,i+1] # no hace falta mandar el elmento␣
,→ A[i+1,i] a cero
v[i+1]-=A[i+1,i]*v[i] # es suficiente asumirlo
# en la sustitución hacía atrás

# 3. Nuestro for va hasta N-2, sólo nos falta operar sobre el último␣
,→elemento de V

v[N-1]/=A[N-1,N-1]

# 4. Sustitución hacía atrás.

x=zeros(N,float)
x[N-1]=v[N-1]
for i in range(N-2,-1,-1):
x[i]=v[i]-A[i,i+1]*x[i+1]
return x

[39]: from numpy import zeros

# Definimos las constantes del sistema

N=30
C=1
m=1
k=5
omega=2

alpha=2*k-m*omega**2

# Inicializamos nuestra matriz y la matriz V

A=zeros([N,N],float)
v=zeros(N,float)

for i in range(N-1):
A[i,i]=alpha # Dado el valor de un elemento de la diagonal.
A[i+1,i],A[i,i+1]=-k,-k # le damos valores al elemento arriba y a la␣
,→izquierda.

A[0,0]+=-k # El primer y el último elemento son␣

,→diferentes.

A[N-1,N-1]=alpha-k
v[0]=C # El único elemento no nulo de V es el primero.

x=tridiagonal_gaussian(A,v)

20
[40]: from [Link] import plot,show

plot(x)
plot(x,"ro")
show()

[41]: # Para optimizar el código, una matriz a banda con up=2 y down =2 estará escrita␣
,→como:

# ( - - A02 A13 A24 ...

# ( - A01 A12 A23 A34 ...
# ( A00 A11 A22 A33 A44 ...
# ( A10 A21 A32 A43 A54 ...
# ( A20 A31 A42 A53 A64 ...

# donde los elementos con guiones no son relevantes para el cálculo.

# Los elementos de la diagonal están en la fila up.
# El número de filas será 1+up+down.
# El número de columnas será N.

def band_gaussian(A,v,up,down):

N = len(v)

# Aplicamos la eliminación gaussiana

21
 for m in range(N):

# dividimos por el elemento m de la diagonal

# que en nuestra interpretación es el elemento A[up,m]
div = A[up,m]

# Actualizamos el vector v dividiendo por div,

# y sustrayendo los elementos por debajo.
v[m] /= div
for k in range(1,down+1): # para todos los elemementos por␣
,→debajo down veces

if m+k<N: # del elemento m de la diagonal

v[m+k] -= A[up+k,m]*v[m] # procedemos con la eliminación␣
,→gaussiana

# Ahora aplicamos la sustracción a los elementos de A.

for i in range(up):
j = m + up - i
if j<N:
A[i,j] /= div
for k in range(1,down+1):
A[i+k,j] -= A[up+k,m]*A[i,j]

# Finalmente hacemos la sustitucioń hacía

for m in range(N-2,-1,-1):
for i in range(up):
j = m + up - i
if j<N:
v[m] -= A[i,j]*v[j]

return v

[42]: # Creamos nuestra matriz banda siguiendo la notacioń

A_band=zeros([3,N],float)
v=zeros(N,float)

A_band[0,:]=-k
A_band[1,:]=alpha
A_band[2,:]=-k

A_band[1,0]=alpha-k
A_band[1,N-1]=alpha-k

v[0]=C

x2= band_gaussian(A_band,v,1,1)

22
[43]: plot(x2)
plot(x,"ro")
show()

[44]: from vpython import canvas,sphere,vector,rate

from numpy import empty,cos

scene=canvas(width=300, height=300,
center=vector(0,0,0))

s=empty(N,sphere)
mu=empty(N,float)

for i in range(N):
mu[i]=2*i-N
s[i]=sphere(pos=vector(mu[i],0,0),radius=0.5)

t=0.0
tmax=20
framerate=30

while t<tmax:
rate(framerate)
t+=1/framerate
for i in range(N):
s[i].pos=vector(mu[i]+x[i]*cos(omega*t),0,0)

23
<[Link] object>

<[Link] object>

<[Link] object>

<[Link] object>

<[Link] object>

<[Link] object>

<[Link] object>

<[Link] object>

1.6. Métodos iterativos para ecuaciones lineales

La factorización LU conlleva 32 N 3 operaciones, lo que para sistemas de ecuaciones grandes
puede ser muy costoso computacionalmente.
Como alternativa, se suelen usar los métodos iterativos, que en general sólo son prácticos en
ciertas situaciones, por ejemplo con matrices con muchos elementos no diagonales nulos.

1.6.1. El método de Jacobi

El método de Jacobi se basa en descomponer una matriz cuyos elementos diagonales no sean
nulos (hay que aplicar el pivote si no es el caso) en:

A = D+L+U

donde D es diagonal y L/U sólo tiene elementos no nulos por debajo/encima de la diagonal.
Por tanto el sistema de ecuaciones A · x = v se puede escribir como

A · x = (D + L + U) · x = v ⇒ D · x = v − ( L + U ) · x,

Como D no tiene determinante nulo, es invertibles y por tanto:

!
1
x=D −1
(v − ( L + U ) · x ) ⇒ xi = vi − ∑ aij x j , con i = 1 · · · N.
aii j 6 =i

En vez de buscar la solución exacta, la ecuación permite desarrollar un método iterativo.

24
 Definiendo un punto de partida x (0) , por ejemplo x (0) = 0 si no hay más información, obte-
nemos para la estimación k-esima es:
!
1
xik+1 = vi − ∑ aij x kj ,
aii j 6 =i

Sin embargo el sistema sólo convergerá si se cumple la condición:

D −1 · ( L + U ) < 1.

Una condición suficiente pero no necesaria para su convergencia viene dada por la condi-
ción:

aii > ∑ |aij |,

j 6 =i

es decir, que sea una matriz diagonal dominante.

Ejemplo 5.9: aplicando el método de Jacobi Aplicar el método de Jacobi al sistema de ecuacio-
nes:

8 w + x + 2 y + z = − 1,
w + 5 x − y − z =3,
w − 4 x + 6y + z =2,
w − 2 x + y + 5 z =1.
(7)

Calcular el número de iteraciones necesarias para obtener una precisión mayor de 10−6 .

[5]: from numpy import zeros,diag,diagflat,dot,copy,array

def jacobi(A,v,eps):

# eps es el error de convergencia

N=len(v)
# Creamos nuestra descompsición D,L,U
D=diag(A) # diag crea un vector con los elementos de la␣
,→diagonal

LU=A-diagflat(D) # diagflat crea una matriz diagonal cuyo argumento␣

,→es un vector.

x0=zeros(N,float) # creamos la primera estimación de X

err=1e6
it=0
while err>eps:

25
 x=(v-dot(LU,x0))/D # iteramos con X
err=max(abs(x-x0))
x0=copy(x)
it+=1
return x,it

[6]: A=array([[8,1,2,1],[1,5,-1,-1],[1,-4,6,1],[1,-2,1,5]],float)
V=array([-1,3,2,1],float)
x1,it=jacobi(A,V,1e-6)
x2,AB=gaussian_elim(A,V)
print(x1,x2,sep=", ")
print("El número de itracciones del método de jacobi es: ",it)

[-0.57409676 1.02366086 1.02490595 0.51930229], [-0.57409714 1.02366127

1.0249066 0.51930262]
El número de itracciones del método de jacobi es: 25

1.6.2. Método de Gauss-Seidel

Corresponde a una variación del método de Jacobi, consistente en notar que la matriz D+L,
es triangular inferior y por tanto se puede resolver por sustitución hacía adelante:

( L + D ) · x = v − U · x,

Por tanto, partiendo de una estimación inicial para x, podemos resolver el sistema exacta-
mente, lo que nos permite acelerar la convergencia del método de jacobi.
!
i −1 N
1
xik+1 = vi − ∑ aij x kj +1 − ∑ aij x kj ,
aii j =1 j = i +1

que converge si

( D + L)−1 · U < 1,

condiciones que se cumplen si A es positiva definida o dominada por los elementos de la

diagonal.

Ejemplo 5.10: aplicando el método de Gauss-Seidel Aplicar el método de Gauss-Seidel al siste-

ma del ejemplo 5.9.
Calcular cuantas iteraciones son necesarias para la misma precisión.

[7]: from numpy import tril

def gauss_seidel(A,v,eps):

# eps es el error de convergencia

26
 N=len(v)
# Creamos nuestra descompsición D,L,U
DL=tril(A) # tril crea un matriz triangular inferior
U=A-DL # creamos la matrix U
x0=zeros(N,float) # creamos la primera estimación de x
# inicializamos x
err=1e6
it=0

while err>eps:
x=(v-dot(U,x0))
# Sustitución hacía adelante
for m in range(N):
for i in range(m):
x[m]-= DL[m,i]*x[i]
x[m]/=DL[m,m]

err=max(abs(x-x0))
x0=copy(x)
it+=1

return x,it

[8]: A=array([[8,1,2,1],[1,5,-1,-1],[1,-4,6,1],[1,-2,1,5]],float)
V=array([-1,3,2,1],float)
x1,it=gauss_seidel(A,V,1e-6)
x2,AB=gaussian_elim(A,V)
print(x1,x2,sep=", ")
print("El número de itracciones del método de Gauss-Seidel es: ",it)

[-0.57409702 1.02366117 1.02490654 0.51930256], [-0.57409714 1.02366127

1.0249066 0.51930262]
El número de itracciones del método de Gauss-Seidel es: 13

1.6.3. Método de Sobrerelajación sucesiva

El método de Gauss-Seidel se puede escribir como:
!
i −1 N
1
xik+1 = xik + vi − ∑ aij x kj +1 −∑ aij x kj ,
aii j =1 j =i

donde simplemente hemos sumado y restado la estimación x k en la solución.

Los métodos de relajación usan que la combinación D + ωL, siendo ω un parámetro arbi-
trario, es una matriz triangular inferior para generalizar el método de Gauss-Seidel para el
caso ω 6= 0:

27
 !
i −1 N
ω
xik+1 = (1 − ω ) xik + vi − ∑ aij x kj +1 − ∑ aij x kj ,
aii j =1 j = i +1

que permite encontrar el valor de ω que optimiza la convergencia.

Además, para valores de 0 < ω < 2 el método de sobrerelajación sucesiva converge para
cualquier matriz simétrica definida positiva.

Ejemplo 5.11: aplicando el método de la sobrerelajación sucesiva Aplicar el método de la so-

brerelajación sucesiva al sistema de ecuaciones del ejemplo 5.9.
Calcular una vez más el número de itercciones necesarias para alcanzar la misma precisión.

[9]: from numpy import tril

def gauss_seidel_rel(A,v,w,eps):

# eps es el error de convergencia

N=len(v)
# Creamos nuestra descompsición D,L,U
DL=tril(A) # tril crea un matriz triangular inferior
U=A-DL # creamos la matrix U
x0=zeros(N,float) # creamos la primera estimación de x
# inicializamos x
err=1e6
it=0

while err>eps:
x=(v-dot(U,x0))
# Sustitución hacía adelante
for m in range(N):
for i in range(m):
x[m]-= DL[m,i]*x[i]
x[m]=(1-w)*x0[m]+w/DL[m,m]*x[m]
err=max(abs(x-x0))
x0=copy(x)
it+=1

return x,it

[10]: A=array([[8,1,2,1],[1,5,-1,-1],[1,-4,6,1],[1,-2,1,5]],float)
V=array([-1,3,2,1],float)
w=1.1
x1,it=gauss_seidel_rel(A,V,w,1e-6)
x2,AB=gaussian_elim(A,V)
print(x1,x2,sep=", ")
print("El número de itracciones del método de la sobrerelajación sucesiva es:␣
,→",it)

28
[-0.57409704 1.02366121 1.02490661 0.51930259], [-0.57409714 1.02366127
1.0249066 0.51930262]
El número de itracciones del método de la sobrerelajación sucesiva es: 9

2. Autovalores y autovectores
Los problemas de autovalores son fundamentales en muchos campos de la física.
En la mayoría de los casos físicos, involucran matrices reales simétricas o complejas hermí-
ticas.
En esos casos la matriz es diagonalizable, y el problema es equivalente a encontrar el sistema
de coordenadas para el que la matriz es diagonal.
En particular, para una matriz simétrica de dimensión N × N el problema de autovalores:

A · vi = λi vi , i = 1, · · · , N

con vi el autovector i y λi ∈ R su correspondiente autovalor, se puede escribir como:

A · V = D · V,

con D una matriz diagonal cuyos valores son los autovalores λi y V una matriz ortogonal,
es decir que satisface,

V T · V = V · V T = I,

cuya columna i-esima corresponde al autovector vi .

2.1. Descomposición QR
El método más común usado para calcular los autovalores de una matriz es la descomposi-
ción QR.
Es parecida a la factorización LU, pero en este caso la matriz A se descompone como

A = Q · R,

siendo Q una matriz ortogonal y R una matriz triangular superior.

La factorización QR corresponde a aplicar el método de Gramm-Schmidt para la obtención
de una base ortonormal.
Para ello entendamos la matriz A como un conjunto de N vectores columna a1 , · · · , a N

 
| | ··· |
A =  a1 a2 · · · an 
| | ··· |

29
 que podemos usar para construir los vectores ui y qi con i = 1, · · · , N definidos por

u1
u1 = a1 , q1 = ,
| u1 |
u2
u2 = a2 − ( q1 · a2 ) q1 , q2 = ,
| u2 |
u3
u3 = a3 − ( q1 · a3 ) q1 − ( q2 · a3 ) q2 , q3 = ,
| u3 |
(8)

y de forma general

i −1
ui
ui = ai − ∑ a j · q j q j ,


qi = ,
j =1
| ui |

Los vectores qi forma una base de vectores ortonormales, que permiten escribir:

  

| | ··· |
 
| | ··· | | u1 | ( q1 · a2 ) · · · ( q1 · a n )
A=  a1 a2 · · · an  =  q1 q2 · · · qn   0
· | u2 | ··· ( q2 · a n ) 
 = Q·R
| | ··· | | | ··· | .
.. .
.. .. ..
. .

2.1.1. Ejemplo 5.12: estudiando la descomposición QR

Dada la matriz

 
1 4 8 4
 4 2 3 7 
A=
 8
,
3 6 9 
4 7 9 2

Escribir un programa que la factorice en las matrices Q y R.

[2]: from numpy import array,zeros,dot,sqrt,identity,copy,absolute,max

from [Link] import norm

# Función para calcular la factorización QR

def QR_fact(A):

N=len(A)

# Inicializamos los tres vectores que necesitamos para Gram-Schmidt

U = zeros([N,N],float)
Q = zeros([N,N],float)

30
 R = zeros([N,N],float)

for m in range(N):
U[:,m] = A[:,m] # Creamos los vectores de U
for i in range(m):
R[i,m] = dot(Q[:,i],A[:,m]) # los vamos calculando en cada␣
,→iteracción

U[:,m] -= R[i,m]*Q[:,i]
R[m,m] = norm(U[:,m])
Q[:,m] = U[:,m]/R[m,m]
return Q,R

[3]: A=array([[1,4,8,4],
[4,2,3,7],
[8,3,6,9],
[4,7,9,2]],float)

Q,R=QR_fact(A)
print(Q)
print()
print(R)
print()
print(dot(Q,R)-A)

[[ 0.10153462 0.558463 0.80981107 0.1483773 ]

[ 0.40613847 -0.10686638 -0.14147555 0.8964462 ]
[ 0.81227693 -0.38092692 0.22995024 -0.37712564]
[ 0.40613847 0.72910447 -0.5208777 -0.17928924]]

[[ 9.8488578 6.49821546 10.55960012 11.37187705]

[ 0. 5.98106979 8.4234836 -0.484346 ]
[ 0. 0. 2.74586406 3.27671222]
[ 0. 0. 0. 3.11592335]]

[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]

[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 8.88178420e-16 -1.77635684e-15]
[ 0.00000000e+00 -8.88178420e-16 0.00000000e+00 0.00000000e+00]]

2.2. Descomposiciones QR sucesivas

Una vez hemos factorizado nuestra matriz A en una matriz ortogonal Q1 y una matriz trian-
gular R1

A = Q1 · R1 ,

multiplicando Q1 por la derecha y por Q1T por la izquierda tenemos

31
 Q1T · A · Q1 = Q1T · Q1 · R1 · Q1 = R1 · Q1 ,

donde hemos usado que Q1 es una matriz ortogonal.

Llamando A1 = Q1T · A · Q1 , lo relevante es notar que A y A1 tienen los mismos autovalores
(el mismo espectro) y son por tanto matrices semejantes.
Llamando v1 y λ1 a uno de los autovectores y autovalores de A1 tenemos que

A1 · v1 = Q1T · A · Q1 · v1 = λ1 v1 ⇒ A · ( Q1 · v1 ) = λ1 ( Q1 · v1 ) ,

y por tanto λ1 también es un autovalor de A con autovector Q1 · v1 .

Esto implica que el problema de autovalores sigue siendo el mismo, y por tanto nos podemos
centrar en la nueva matriz A1 , la cual podemos descomponer una vez más

A1 = Q1T · A · Q1 = Q2 · R2 ,

y volviendo a multiplicar por Q2 y Q2T por la derecha e izquierda respectivamente, encon-

tramos una nueva matriz A2

A2 = Q2T · A1 · Q2 = Q2T · Q1T · A · Q1 · Q2 = R2 · Q2 ,

Podemos repetir este proceso n veces obteniendo en cada iteración:

A1 = Q1T · A · Q1 ,
A2 = Q2T · Q1T · A · Q1 · Q2 ,
A3 = Q3T · Q2T · Q1T · A · Q1 · Q2 · Q3 ,
.. ..
. .,
 
An = QnT · · · Q1T · A · ( Q1 · · · Qn ) ,
(9)

Una matriz ortogonal no es más que un cambio de coordenadas, es decir, una rotación.
Por ejemplo, en dos dimensiones, N=2:

 
cos α − sin α
Q= ,
sin α cos α

Al rotar el sistema de coordenadas de nuestra matriz k veces, estamos transformando nues-

tra matriz en una matriz diagonal.

[6]: from [Link] import HTML

HTML('<img src="[Link]">')

32
[6]: <[Link] object>

La demostración de que en el limite cuando k → ∞ la matriz resultante es diagonal, es

complicada.
Sin embargo, se puede entender notando que con cada rotación, mientras que fuera de la
diagonal estamos en cada iteración multiplicando por un producto de senos y cosenos, que
acaban cancelándose, en la diagonal siempre sumamos senos y cosenos al cuadrado.
Por tanto, la idea es continuar el proceso hasta que la matriz Ak se vuelva diagonal.
Los términos fuera de la diagonal se van haciendo cada vez más pequeños hasta que se
hacen cero o tan cercanos a cero que los podemos despreciar.
Por tanto, dada una precisión podemos diagonalizar la matriz mediante transformaciones
de semejanza hasta obtener una matriz diagonal.
Por tanto, definiendo la matriz ortogonal:

k
V = ( Q1 · · · Q n ) = ∏ Qi ,
i =1

obtenemos que:

V T · A · V = D ⇒ A · V = D · V.

2.3. Algoritmo QR para el cálculo de autovectores y autovalores

Elegida una precisión e para los términos de fuera de la diagonal.
1. Crear una matriz V N × N para contener los autovectores. Inicializar V igual a la matriz
identidad.
2. Calcular la factorización QR de A: A = Q · R.
3. Actualizar A a un nuevo valor A = R · Q.
4. Multiplicar V por la derecha por Q.
5. Comprobar los términos no-diagonales de la nueva A. Si son en valor absoluto menores que
e, hemos acabado.
En caso contrario, volver al paso 2.

2.3.1. Ejemplo 5.13: calculando autovalores y autovectores

Dada la matriz del ejemplo 5.12

 
1 4 8 4
 4 2 3 7 
A=
 8
,
3 6 9 
4 7 9 2

33
 Escribir un programa que utilice la factorización QR para obtener los autovalores de A asegurando
que los elementos de fuera de la diagonal son menores de 10−6 .

[4]: A=array([[1,4,8,4],
[4,2,3,7],
[8,3,6,9],
[4,7,9,2]],float)

Q,R=QR_fact(A)
print(Q)
print()
print(R)
print()
print(dot(Q,R)-A)

[[ 0.10153462 0.558463 0.80981107 0.1483773 ]

[ 0.40613847 -0.10686638 -0.14147555 0.8964462 ]
[ 0.81227693 -0.38092692 0.22995024 -0.37712564]
[ 0.40613847 0.72910447 -0.5208777 -0.17928924]]

[[ 9.8488578 6.49821546 10.55960012 11.37187705]

[ 0. 5.98106979 8.4234836 -0.484346 ]
[ 0. 0. 2.74586406 3.27671222]
[ 0. 0. 0. 3.11592335]]

[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]

[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 8.88178420e-16 -1.77635684e-15]
[ 0.00000000e+00 -8.88178420e-16 0.00000000e+00 0.00000000e+00]]

[5]: # Algoritmo QR para la diagonalizcion

def QR_diag(A,eps):

N=len(A)
V = identity(N,float) # Empezamos con la identidad
delta = 1.0
while delta > eps:
# Hacemos un paso del algoritmo QR

Q,R = QR_fact(A)
A = dot(R,Q)
V = dot(V,Q)

# Encontramos el mayor valor absoluto de los elementos fuera de la␣

,→diagonal

Ac = copy(A)
for i in range(N):

34
 Ac[i,i]=0.0
delta = max(absolute(Ac))
return A,V

[6]: from numpy import diag

D,V=QR_diag(A,1e-12)

print(D)
print()
print(V)
print()
print("Nuestros autovalores son: ",diag(D))

[[ 2.10000000e+01 9.83584264e-13 2.41644555e-15 -4.44737642e-14]

[ 9.97947346e-13 -8.00000000e+00 5.70355716e-14 -5.87329266e-15]
[ 2.60052149e-26 6.51467609e-14 -3.00000000e+00 -7.21588694e-15]
[ 1.75424525e-41 2.17973835e-28 -4.31727622e-15 1.00000000e+00]]

[[ 0.43151697 -0.38357064 -0.77459667 -0.25819889]

[ 0.38357064 0.43151697 -0.25819889 0.77459667]
[ 0.62330229 0.52740963 0.25819889 -0.51639778]
[ 0.52740963 -0.62330229 0.51639778 0.25819889]]

Nuestros autovalores son: [21. -8. -3. 1.]

[7]: # Comprobemos que son autovalores y autovectores:

for i in range(len(A)):
print(dot(A,V[:,i])- diag(D)[i]*V[:,i])

[-3.80140364e-13 4.36983782e-13 5.18696197e-13 -6.11066753e-13]

[3.72146758e-13 3.56159546e-13 6.27942143e-13 5.48894263e-13]
[-2.08721929e-14 2.13162821e-14 3.24185123e-14 -3.77475828e-14]
[-1.18238752e-14 -1.74305015e-14 -3.24185123e-14 -2.32591724e-14]

2.3.2. Autovalores y autovectores en python

En Python en el módulo [Link] hay funciones para calcular las autofunciones y auto-
valores de una matriz A:
• eigh calcula los autovalores y autovectores de un matriz simétrica o hermítica optimi-
zando el algoritmo QR para este tipo de matrices.
• eigvalsh calcula solo los autovalores. Devuelve un array con los autovalores.
• Para matrices no simétricas o no hermíticas uno puede usar las funciones: eig y eigvals.

35
 2.3.3. Ejercicio 5.2: Pozo de potencial asimétrico en Mecánica Cuántica
La Mecánica Cuántica puede formularse como un problema matricial que puede resolverse en un
ordenador usando los métodos del álgebra lineal.
Estudiemos por ejemplo una partícula de masa m en un potencial unidimensional de de longitud
L descrito por una función V ( x ).

[1]: from [Link] import Image,display

display(Image(filename="[Link]", width=600, height=600))

En general estos problemas no pueden resolverse analíticamente y tienen que ser analizados en
un ordenador. Un autoestado de energía E satisface la ecuación de Schrödinger independiente del
tiempo:

Ĥψ( x ) = E ψ( x ),

donde el hamiltoniano satisface:

36
 h̄2 d2
Ĥ = − + V ( x ).
2M dx2

Asumamos que las paredes del potencial tienen una altura infinita, que nos permite definir que la
función de onda se anula para x = 0 y L.
En ese caso la función de onda puede escribirse como una series de Fourier:

∞
πnx
ψ( x ) = ∑ ψn sin L
,
n =1

donde ψi con i = 1, · · · , ∞ son los coeficientes de Fourier.

1. Teniendo en cuenta que para dos enteros positivos m y n

ZL 
πmx πnx L/2 si m = n
sin sin dx =
L L 0 de otra forma.
0

mostrar que la ecuación de Schrödinger Ĥψ = Eψ implica que

∞ ZL
πmx πnx 1
∑ ψn sin
L
Ĥ sin
L
dx = LE ψm .
2
n =1 0

2. Definiendo la matriz H como aquella con elementos:

ZL
2 πmx πnx
Hmn = sin Ĥ sin
L L L
0
ZL
" #
2 πmx h̄2 d2 πnx
= sin − 2
+ V ( x ) sin ,
L L 2M dx L
0

muestra que la ecuación de Schrödinger se puede escribir de forma matricial

H · ψ = Eψ

siendo ψ = (ψ1 , ψ2 , · · · ), que por tanto es un autovector de la matriz del Hamiltoniano H

con autovalor E.
Conociendo los autovalores de la matriz estamos encontrando los diferentes niveles de ener-
gía posibles de nuestra partícula.

37
 3. Asumiendo ahora que V ( x ) = ax/L, resolver la integral analíticamente,que permite encon-
trar una expresión para el elemento de matriz Hmn . Mostrar que la matriz es real y simétrica.
Tener en cuenta que:


ZL  0 si m 6= n y ambos son pares o impares
2L 2
πmx πnx 
mn


x sin sin dx = − 2 si m 6= n y uno es par y otro impar
L L  π ( m2 − n2 )
0
 2
L /4 si m = n

4. Escribir un programa que calcule Hmn para cualquier valor de m y n para un electrón en un
pozo de anchura L = 5 Å y a = 10 eV.
(Ayuda me = 9,1094 × 10−31 Kg y e = 1,6022 × 10−19 C).
5. La matriz H es infinitamente grande por lo que no es posible calcular sus autovalores. Sin
embargo, se puede obtener una solución aproximada muy buena si sólo consideramos los
primeros autovalores.
Crear una matriz 10 × 10 de los elementos de matriz hasta m, n = 10 y calcular sus autova-
lores y auto vectores.
6. Hacerlo ahora una matriz 100 × 100 y comparar los resultados.
7. Representar la densidad de probabilidad |ψ( x )|2 para función de onda ψ( x ) normalizada,
para el estado fundamental y los dos primeros estados excitados.

[ ]:

3. Ecuaciones no lineales
Por ahora sólo hemos aprendido a resolver sistemas de ecuaciones lineales.
Pero muchas ecuaciones que aparecen en Física son no-lineales.
Las ecuaciones no lineales son mucho más complicadas de resolver.
En particular, resolver una única ecuación no lineal puede ser ya complicado, por lo que
resolver varias a la vez aún más dífícil.
En este apartado veremos algunas técnicas para enfrentarnos a estos problemas.
Empezaremos por el caso de una única ecuación no-lineal.

3.1. Método del punto fijo o relajación

Se aplica a ecuaciones del tipo:

x = f ( x ),

y consiste en busca una solución iterativamente. Es decir, partiendo de una estimación de la

solución x0 la siguiente es:

38
 x1 = f ( x0 ),

de forma que tras k iteraciones, la solución exacta se ha encontrado si:

xk+1 − xk = f ( xk ) − xk = 0.

Lo relevante es por tanto estudiar bajo que condiciones ese proceso iterativo converge a la
solución exacta.
Para ello supongamos que tenemos una estimación de la solución xi y definamos la solución
exacta como x̄.
Expandiendo en serie de Taylor en torno a la solución, se tiene:

xi+1 = f ( xi ) = f ( x̄ ) + ( xi − x̄ ) f 0 ( x̄ ) + · · · = x̄ + ( xi − x̄ ) f 0 ( x̄ ) + · · · ,

por tanto, la diferencia entre nuestra siguiente estimación y la solución exacta es:

xi+1 − x̄ = ( xi − x̄ ) f 0 ( x̄ ) + · · · ,

En general, para funciones suficientemente suaves, la distancia de las diferentes estimacio-

nes a la solución exacta irá disminuyendo si y sólo si

| f 0 ( x̄ ) | < 1.

3.1.1. Ejemplo 5.14: aplicando el método del punto fijo

Aplicar el método del punto fijo para encontrar la solución a la ecuación

x = 2 − e− x ,

que no tiene solución analítica.

¿Por qué converge el método?

[15]: from numpy import exp

def f(x):
return 2-exp(-x)

eps=1e-6 # definimos la precisión objetivo

it=0
err=1.0
x0=1.0
while err>eps:
x1=f(x0)
err=abs(x1-x0)

39
 it+=1
x0=x1
print(it,x1)

1 1.6321205588285577
2 1.8044854658474119
3 1.8354408939220457
4 1.8404568553435368
5 1.841255113911434
6 1.8413817828128696
7 1.8414018735357267
8 1.8414050598547234
9 1.8414055651879888

[16]: # Veamos cual es el valor de la derivada en la solución

h=1e-5
derx0=(f(x0+h/2)-f(x0-h/2))/h

print(derx0,exp(-x0),sep=", ")

0.1585943546711377, 0.15859435466898786

3.1.2. No convergencia del método del punto fijo

Si la derivada de la función en la solución cumple

| f 0 ( x̄ )| > 1,

el método del punto fijo no convergerá hacía la solución.

Sin embargo, si la inversa de f existe,

x = f ( x ) ⇒ x = f −1 ( x ),

y desarrollando en serie de Taylor en torno a la solución se obtiene:

0 1
xi+1 = f −1 ( xi ) = f −1 ( x̄ ) + ( xi − x̄ ) f −1 ( x̄ ) + · · · = x̄ + ( xi − x̄ ) +··· ,
f 0 ( x̄ )

que por tanto convergerá a la solución si:

| f 0 ( x̄ ) | > 1.

40
 Ejemplo 5.15: aplicando el método del punto fijo por segunda vez aplicar el método del punto
fijo para encontrar la solución a la ecuación

2
x = e1− x .

La ecuación no tiene solución analítica. Sin embargo se ve trivialmente que x = 1 es una solución.
La derivada de la función en ese punto es:

2
f 0 (x) x =1
= − 2 x e1 − x = −2,
x =1

y por tanto, la aplicación directa del método del punto fijo no va a funcionar. Comprobémoslo con
las 10 primeras iteraciones:

[27]: from math import exp

x = 0
eps = 1e-6
for k in range(10):
x = exp(1-x**2)
err= abs(x-exp(1-x**2))
if err < eps:
break
print(k+1,x,sep=", ")

1, 2.718281828459045
2, 0.0016798410570681988
3, 2.71827415784286
4, 0.001679911110814363
5, 2.718274157203078
6, 0.00167991111665744
7, 2.7182741572030253
8, 0.0016799111166579221
9, 2.7182741572030253
10, 0.0016799111166579221
La estimación va oscilando continuamente sin obtener ningún patrón de convergencia.
Haciendo uso de que la inversa de la función existe, podemos reescribir el problema como:

p
x= 1 − log x,

cuya solución si podemos obtener por el método del punto fijo.

[29]: from numpy import sqrt,log

def f(x):
return sqrt(1-log(x))

41
eps=1e-6 # definimos la
it=0
err=1.0
x0=0.5
while err>eps:
x1=f(x0)
err=abs(x1-x0)
it+=1
x0=x1
print(it,x1)

1 1.3012098910475378
2 0.8583154914892762
3 1.0736775779454883
4 0.9637999044091371
5 1.0182689104343374
6 0.990906635925747
7 1.004557096969838
8 0.997724037576543
9 1.0011386299421705
10 0.9994308469350205
11 1.000284617043603
12 0.9998577016016549
13 1.000071151730577
14 0.9999644247674951
15 1.0000177877744567
16 0.9999911061523217
17 1.0000044469337268
18 0.9999977765356085
19 1.0000011117328136
20 0.9999994441337476
21 1.000000277933165

3.1.3. Limitaciones del método del punto fijo

Puede haber más de una solución:
• En ese caso el valor de la solución encontrada por el método de relajación dependerá
del valor inicial.
• Así, si tenemos una idea aproximada de la solución que estábamos buscando debere-
mos usarla para el valor inicial.
Hay funciones cuya inversa no existe o no sé puede determinar analíticamente.
• En ese caso, si | f 0 ( x̄ )| > 1 no será posible encontrar una solución.
Sin embargo, hay ocasiones en la que, incluso en ese caso, aún es posible encontrar una
solución reordenando parcialmente la función.

42
 Ejemplo 5.16: aplicando el método del punto fijo por tercera vez dada la ecuación:

x = x2 + sin 2x,

cuya solución es x = 0, explicar si es resoluble por medio del método de relajación/punto fijo.
En caso negativo, intentar reescribir la ecuación de forma que el método converja.
La derivada de la función en x=0 es | f 0 ( x )| = 2, por lo que el método de relajación no funcionará.
Además, la función f ( x ) = x2 + sin 2x no es invertible.
Sin embargo, reescribiendola como:

1
sin−1 x − x2 ,


x = g( x ) =
2

se puede comprobar que en ese caso: | g0 ( x )| = 1/2 por lo el método del punto fijo en este caso si
va a converger hacía la solución.

[4]: from numpy import arcsin

def f(x):
return arcsin(x-x*x)/2

eps=1e-6 # definimos la
it=0
err=1.0
x0=1e-4
while err>eps:
x1=f(x0)
err=abs(x1-x0)
it+=1
x0=x1
print(it,x1)

1 4.999500008330834e-05
2 2.4996250302049482e-05
3 1.249781274606156e-05
4 6.248828275531731e-06
5 3.12439461385879e-06
6 1.5621924260110852e-06
7 7.810949927832724e-07

3.1.4. Convergencia y estimación del error del método del punto fijo/relajación
Una vez hemos visto las condiciones por las cuales el método del punto fijo converge, el
paso siguiente es saber como de rápido lo hace.
Dada una iteración i la diferencia con la solución exacta es

43
| xi+1 − x̄ | = | xi − x̄ | | f 0 ( x̄ ) | + · · · = | xi−1 − x̄ | | f 0 ( x̄ ) |2 + · · · = · · · = | x0 − x̄ | | f 0 ( x̄ ) |i+1 + · · · ,
0 ( x̄ )| 0
= | x0 − x̄ | e(i+1) log | f + · · · = | x0 − x̄ | e−(i+1) log(1/| f (x̄)|) + · · ·
(10)

que implica que el método del punto fijo converge exponencialmente.

Aunque sabemos que el método converge rápidamente, también es interesante saber cuantas
iteraciones necesitamos hasta llegar a una determinada precisión. O lo que es lo mismo, cual
es el error de una estimación.
Dada la estimación i-esisma, sabemos que la solución exacta es:

x̄ = xi + ei = xi+1 + ei+1 ,

sin embargo ya hemos visto que:

ei + 1
ei = ,
f 0 ( x̄ )

por lo que se obtiene que

ei + 1 x i − x i +1
x i + 1 + ei + 1 = x i + 0
⇒ ei + 1 = ,
f ( x̄ ) 1 − 1/ f 0 ( x̄ )

Sin embargo, en general no conocemos la solución exacta, que es la que estamos buscando.
Aproximando que nuestra estimación k-esima no esta lejos de la solución, podemos desarro-
llar f 0 ( x̄ ) en serie de Taylor

f 0 ( x̄ ) = f 0 ( xi ) − ( x̄ − xi ) f 00 ( xi ) + · · · ' f 0 ( xi ),

que nos permite aproximar que:

x i − x i +1
ei + 1 ' .
1 − 1/ f 0 ( xi )

Por tanto sólo tenemos que repetir la iteración hasta que ei+1 < e siendo e nuestra tolerancia,
sin necesidad de calcular la siguiente iteración i + 2.
Sin embargo, en el caso de no conocer la forma funcional de f ( x ) tendremos que estimar su
derivada numéricamente:

f ( x i ) − f ( x i −1 ) x − xi
f 0 ( xi ) ' = i +1 ,
x i − x i −1 x i − x i −1

que nos permite aproximar el error de estimación i + 1

44
 ( x i − x i +1 )2
ei + 1 ' ,
2xi − xi+1 − xi−1

a partir de las dos estimaciones anteriores.

3.1.5. Ejercicio 5.3: propagación de epidemias

La ecuación x = 1 − e−cx aparece en modelos de procesos de contacto, de epidemias y en la teoría
de gráficos aleatorios.
1. Escribir un programa que resuelva la ecuación utilizando el método del punto fijo para c =
2. Calcular la solución con una precisión de al menos 10−6 .
2. Modificar el programa para calcular la solución para valores de c de 0 a 3 en pasos de 0.01 y
hacer un plot de x como función de c.
3. ¿Existe una transición de Fase?. (Es la llamada transición de percolación o el umbral epidé-
mico).
[ ]:

3.1.6. El método del punto fijo sobrerelajado

Al igual que con el método de Gauss-seidel para los sistemas de ecuaciones lineales, el mé-
todo de la sobrerelajación puede ayudar a hacer el método del punto fijo más rápido.
El método de sobrerelajación acelera la convergencia del sistema iterativo utilizando en cada
paso la dirección en que la iteración esta yendo. Dada la estimación

x i +1 = f ( x i )

la diferencia entre dos estimaciones es

∆xi+1 = xi+1 − xi = f ( xi ) − xi , ⇒ xi+1 = xi + ∆xi+1 .

El método de sobrerelajación modifica esta ecuación introduciendo un parámetro ω de for-

ma que:

xi+1 = xi + (1 + ω ) ∆xi+1 = xi + (1 + ω ) ( f ( xi ) − xi ) = (1 + ω ) f ( xi ) − ω xi .

Si el parámetro ω > 0 la iteración cambia por un factor ω∆xi+1 más de lo que hace el método
de relajación normal.
Para que el método funcione ω tiene que escogerse correctamente.
Lamentáblemebte no hay una teoría general que permite obtener el valor óptimo de ω, que
normalmente se obtiene por prueba y error.

45
 Ejercicio 5.4: aplicando el método del punto fijo sobrerelajado
1. Obtener la ecuación que da el error de una estimación dada usando el método de sobrerela-
jación.
2. Modificar el programa del Ejercicio 5.3 para que calcule el número de iteracciones necesarias
para llegar a la solución. Calcularlas para c=2.
3. Modificar de nuevo el programa para buscar la solución de la ecuación del Ejercicio 5.3 con
c=2 usando el método de la sobrerelajación. Empezar con un valor de ω = 0,5 y variarlo
hasta encontrar el valor que minimiza el número de iteracciones.
4. ¿Hay alguna situación en la que el método de la sobrerelajación pueda funcionar para ω <
0?.
[ ]:

3.1.7. Método del punto fijo o relajación con más de una variable
El método del punto fijo o relajación se puede extender fácilmente para el estudio de sistemas
de ecuaciones no lineales en varias variables.
Dadas N ecuaciones con N incógnitas:

x1 = f 1 ( x1 , · · · , x N ),
x2 = f 2 ( x1 , · · · , x N ),
.. ..
. .
x N = f N ( x1 , · · · , x N ),
(11)

(0)
partiendo de una estimación inicial para cada una de las variables, xi con i = 1, · · · , N,
aplicamos las ecuaciones iterativamente hasta encontrar una solución.
Como en una dimensión, la convergencia no esta garantizada.
En este caso se puede demostrar que el sistema convergerá si todos los autovalores λi de la
matriz jacobiana J, aquella cuyos elementos de matriz son:

d f k ( x1 , · · · , x N )
Jkl = ,
dxl

evaluados en la solución cumplen que |λi | < 1.

Si el sistema no converge, uno puede modificar el sistema de ecuaciones para intentar con-
seguir que lo haga.

46
 Ejercicio 5.5: El proceso de la glucólisis La glucólisis es el proceso bioquímico que permite oxi-
dar la glucosa para extraer de ella energía y puede ser modelado por las ecuaciones:

dx dy
= − x + ay + x2 y, = b − ay − x2 y,
dt dt

donde x y y representan la concentración de los compuestos adenosín difosfato (ADP) y fructosa-

6-fosfato (F6P o éster de Neuberg).
Este tipo de ecuaciones vienen caracterizadas por los puntos estacionarios, aquellos en los que las
derivadas con respecto a todas las variables son cero, y a partir de los cuales las ecuaciones no
cambian.
Los puntos estacionarios de nuestro modelo para glucólisis vienen dados por:

0 = − x + ay + x2 y, 0 = b − ay − x2 y.
1. Obtener los valores de los puntos estacionarios.
2. Probar que en los puntos estacionarios, las ecuaciones de la glucólisis se pueden escribir
como:

b
x = y ( a + x 2 ), y= .
a + x2

3. Escribir un programa para obtener la solución de la glucólisis en los puntos estacionarios

usando el método del punto fijo.
Aplicarlo para a = 1 y b = 2. ¿Converge?.
4. Reordenar la ecuaciones de una forma alternativa de forma que el método sea ahora conver-
gente.

[ ]:

3.2. Método de la bisección

El método de relajación o de punto fijo converge rápidamente y es fácil de programar.
Sin embargo, sólo funciona en ciertas ocasiones.
El método de la bisección, es más robusto para resolver ecuaciones no lineales de una varia-
ble.
Requiere especificar el intervalo en el que se quiere buscar la solución de una ecuación y
garantiza que la encontrará siempre y cuando esta sea única.
Consiste en encontrar una solución de la ecuación:

f ( x ) = 0,

el método se basa por tanto en encontrar las raíces de una ecuación.

47
 3.2.1. Algoritmo para el método de la bisección
1. Definimos el intervalo en el que queremos encontrar una solución (x1 , x2 ).
Comprobamos que los valores de la función en esos puntos, f ( x1 ) y f ( x2 ) tienen signos
diferentes.
En ese caso, sí la función es continua, debe haber como mínimo un punto en el que función
se anula (teorema de Bolzano).
Se dice por tanto que el intervalo [ x1 , x2 ] contiene la raíz.
2. Calculamos el punto intermedio del intervalo, x 0 = 1
2 ( x1 + x2 ), y evaluamos la función
f ( x 0 ).
3. Si f ( x 0 ) tiene el mismo signo de f ( x1 ), hacemos la sustitución x1 = x 0 . En caso contrario
hacemos x2 = x 0 .
4. Si | x1 − x2 | > e, con e la precisión objetivo de nuestra solución, volvemos al paso 2.
En caso contrario nuestra estimación de la solución es el nuevo punto medio: x 0 =
1
2 ( x1 + x2 ).

[16]: from [Link] import display,Image

display(Image(filename="[Link]", width=600, height=600))

48
 3.2.2. Convergencia del método de la bisección
Al igual que con el método del punto fijo la precisión mejora exponencialmente.
La distancia entre dos estimaciones se divide por dos en cada iteración.
Por tanto, llamando e0 = | x1 − x2 | al error asociado a la estimación inicial, al cabo de n
iteraciones:

e0
en = = e0 e−n log 2 .
2n

Esto nos permite calcular el número de iteraciones necesarias N para llegar a una precisión
objetivo e:

e0
N = log2 .
e

Como el logaritmo es una función que crece muy lentamente, esto nos asegura que con un
número de pasos razonable podremos encontrar la solución con una determinada precisión.
Por ejemplo, sí e0 = | x1 − x2 | = 1010 y nuestra precisión objetivo e = 10−10 el número de
iteraciones necesarias será:
[20]: from math import log

log(10**20,2)

[20]: 66.43856189774725

3.2.3. Limitaciones del método de la bisección

Si dado nuestro intervalo inicial f ( x1 ) y f ( x2 ) tienen el mismo signo el método no funciona.
Además esto se puede deber a que no hay solución, o que hay un número par de soluciones.
Lamentablemente no es fácil distinguir entre ambas posibilidades.

[22]: display(Image(filename="[Link]", width=600, height=600))

49
 Igualmente, es muy posible que no se tenga una estimación inicial del intervalo en el que
debe encontrarse la raíz.
Cualquier par de puntos para los que la función tenga diferente signo valdrán, pero encon-
trarlos puede no ser fácil.
Reglas generales que facilitan encontrar el intervalo inicial son:
• Si se tiene información sobre la función que permite encontrar donde puede estar la
raíz, tenerla en cuenta.
Por ejemplo, si se conoce el signo de la función por consideraciones físicas en un deter-
minado intervalo.
• Si se busca una solución en torno a un punto eso ayuda.
Tomamos entonces un entorno alrededor del punto y vamos aumentando el tamaño
hasta que contenga una raíz.
De la misma forma, el método no permite encontrar soluciones de raíces múltiples pares.
Estos problemas son equivalentes a encontrar el mínimo de una función.

[23]: display(Image(filename="[Link]", width=600, height=600))

50
 Finalmente, es un método que no se puede generalizar a funciones de más de una variable.

3.2.4. Ejemplo 5.17: visualizando el método de la bisección

Encontrar la raíz de la función:

f ( x ) = ex − 2,

usando el método de la bisección con una precisión mayor de 10−6 .

[11]: from numpy import exp,linspace

from [Link] import plot,show,axhline,subplots

def f(x):
return exp(x)-2

x=linspace(-2,2,100)
fig,ax=subplots(1,1,figsize=(12, 4)) # subplots(#filas,#columas)

51
 # nos permite crear varias gráficas en una␣
,→ sólo llamada.
# y cambiar facilmente sus propiedades
[Link](x,f(x),lw=1.5)
[Link](0,ls=':',color="k")
ax.set_xticks([-2, -1, 0, 1, 2])
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$f(x) = e^x - 2$', fontsize=18)
show()

[12]: # En el intervalo [-2,2] tenemos seguro una raíz.

# Definimos nuestro intervalo inicial

a,b=-2,2
fa,fb=f(a),f(b)

fig2,ax2=subplots(1,1,figsize=(12, 4))

[Link](x,f(x),lw=1.5)

# Lo añadimos a la figura.

[Link](a,fa, 'ko')
[Link](b,fb, 'ko')
[Link](a,fa + 0.5, r"$a$", ha='center', fontsize=18)
[Link](b,fb + 0.5, r"$b$", ha='center', fontsize=18)
show()

52
[13]: from numpy import sign

# Definimos nuestra tolerancia

eps=1e-6

# Encontramos la raíz usando el método de la bisección y visualizamos

# los pasos del método en el gráfico

fig3,ax3=subplots(1,1,figsize=(12, 4))
[Link](x,f(x),lw=1.5)
[Link](a,fa, 'ko')
[Link](b,fb, 'ko')
[Link](a,fa + 0.5, r"$a$", ha='center', fontsize=18)
[Link](b,fb + 0.5, r"$b$", ha='center', fontsize=18)

it=0

while b-a>eps:
m=a+(b-a)/2
fm=f(m)
[Link](m,fm,'ko')
[Link](m,fm-0.5, r"$m_%d$" %it,ha='center')

it += 1
if sign(fa)==sign(fm):
a,fa=m,fm
else:
b,fb= m,fm
print(it)

[Link](m,fm,'r*',markersize=10)

53
[Link]("Raíz aproximada en %.3f" % m,fontsize=14, family="serif",
xy=(m, fm), xycoords='data',
xytext=(-150, +50), textcoords='offset points',
arrowprops=dict(arrowstyle="->", connectionstyle="arc3, rad=-.5"))

ax3.set_title("Método de la Bisección")

show()

22

3.2.5. Ejercicio 5.6: la constante de desplazamiento de Wien

La ley de radiación de Planck nos dice que la intensidad de radiación por unidad de área y por
unidad de longitud de onda λ para un cuerpo negro a temperatura T es:

2π h c2 λ−5
I (λ) = ,
ehc/λk B T − 1

donde h es la constante de Planck, c es la velocidad de la luz y k B es la constante de Boltzmann.

1. Probar que la longitud de onda de la radiación emitida máxima es solución de la ecuación

hc
5e−hc/λk B T + − 5 = 0.
λk B T

Haciendo el cambio de variable, x = hc/λk B T la ecuación se puede reescribir como

5e− x + x − 5 = 0,

y la longitud de onda de máxima radiación satisface la ley de desplazamiento de Wien:

54
 b
λ= ,
T

siendo

hc
b= ,
kB x

la constante de desplazamiento de Wien y x la solución de la ecuación no lineal anterior.

2. Escribir un programa para obtener el valor de x con una precisión mayor de 10−6 usando el
método de la bisección.
Calcular el valor de la constante de desplazamiento de Wien.
3. La ley del desplazamiento de Wien es la base de la Pirometría, que permite medir de la tem-
peratura de objetos mediante la observación de la radiación térmica que emiten. El método
se utiliza para estimar la temperatura de la superficie de objetos astronómicos, tales como el
Sol.
La longitud de onda del máximo de radiación del Sol esta en λ = 502 nm.
Utilizando los resultados anteriores estime la temperatura de la superficie del Sol.
[Ayuda: h = 6,62607 × 10−34 J·s, c = 2,99792 × 108 m· s−1 , k B = 1,38064 × 10−23 J · K−1 ]

[ ]:

3.3. Método de Newton-Raphson

Tanto el método del punto fijo como el de la bisección son rápidos y simples, pero tienen
inconvenientes.
El método del punto fijo o relajación no siempre converge a la solución.
El método de la bisección requiere definir un intervalo inicial que contenga la raíz de la
función.
El método de Newton-Raphson es incluso más rápido que los dos anteriores.
Consiste al igual que el método de la bisección en buscar una solución a la ecuación

f ( x ) = 0,

es decir, en buscar las raíces de la función, las cuales se estiman extrapolando la pendiente
de la función en un punto dado.

[34]: display(Image("[Link]",width=600,height=600))

55
 Dada una estimación inicial x la nueva estimación viene dada por:

f (x)
x 0 = x − ∆x = x − ,
f 0 (x)

que por tanto requiere conocer la derivada de la función.

Repitiendo iterativamente el procedimiento se obtiene (de forma general) una estimación
mejor y mejor de la raíz.

3.3.1. Error y convergencia del método de Newton-Raphson

Para estimar el error del método de Newton-Raphson, sólo tenemos que expandir una vez
más en serie de Taylor.
Llamando x̄ a la solución exacta a la ecuación y x a una estimación dada, tenemos:

1
f ( x̄ ) = f ( x ) + ( x̄ − x ) f 0 ( x ) + ( x̄ − x )2 f 00 ( x ) + · · · = 0.
2

56
Despejando, se obtiene:

f (x) 00 00
1 2 f (x) 0 1 2 f (x)
x̄ = x − − ( x̄ − x ) + · · · = x − ( x̄ − x ) +···
f 0 (x) 2 f 0 (x) 2 f 0 (x)

siendo x 0 la nueva estimación que se obtiene a partir de x usando el método de Newton-

Raphson.
Esto implica que nuestra raíz se diferenciará de la nueva estimación con un error:

−1 f 00 ( x )
e = x̄ − x 0 = ( x̄ − x )2 0 .
2 f (x)

De la misma forma el error de la nueva estimación x 0 será:

−1 f 00 ( x 0 ) 1 f 00 ( x 0 ) 2
e0 = x̄ − x 00 = ( x̄ − x 0 )2 0 0 = e ,
2 f (x ) 2 f 0 (x0 )

el error tiene una convergencia cuadrática.

Esto implica que converge mucho más rápidamente que el método del punto fijo y de la
bisección.
Asumiendo que el cociente f 00 ( x )/ f 0 ( x ) = k es aproximadamente constante cuando uno esta
cerca de la raíz, se obtiene que después de N iteraciones

 2 N
2 k
eN ' e0 ,
k 2

es decir, el error varia con N como la exponencial de una exponencial.

Por tanto, la diferencia entre dos estimaciones sucesivas será:
 
0 k 2 k
x − x = e − e = e 1 − e ' e,
2 2

siempre y cuando e sea pequeño.

Así, la estimación del error e para un valor x, viene dado por la diferencia con lo estimación
siguiente x 0 − x.
Por tanto, para implementar el método de Newton-Raphson iteramos la fórmula de la apro-
ximación hasta que la diferencia entre dos valores sucesivos sea comparable o más pequeña
que la precisión deseada.

57
 3.3.2. Desventajas del método de Newton-Raphson.
Necesitamos conocer la derivada de f ( x ).
El método no siempre converge:
por ejemplo si f 0 ( x ) es muy pequeña, el método iterativo hará que el error de cada iteración
sea mayor que en la anterior.
Puede ocurrir que la pendiente de la función apunte en el sentido contrario de la raíz.

[35]: display(Image("[Link]",width=600,height=600))

Nota: cada método tiene sus ventajas e inconvenientes.

Lo importante es disponer de métodos diferentes para utilizar en cada caso el más adecuado,
dependiendo del problema concreto.

3.3.3. Ejemplo 5.18: visualizando el método de Newton-Raphson

Encontrar la raíz de la función:

f ( x ) = ex − 2,

58
 usando el método de Newton-Raphson con una precisión mayor de 10−6 .

[37]: from numpy import exp,linspace

from [Link] import subplots,plot,show

# Definimos nuestra función

def f(x):
return exp(x)-2

# Definimos su derivada

def fp(x):
return exp(x)

x=linspace(-2,2,100)
fig,ax=subplots(1,1,figsize=(12, 4)) # subplots(#filas,#columas)
# nos permite crear varias gráficas en una␣
,→sólo llamada.

# y cambiar facilmente sus propiedades

[Link](x,f(x),lw=1.5)
[Link](0,ls=':',color="k")
ax.set_xticks([-2, -1, 0, 1, 2])
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$f(x) = e^x - 2$', fontsize=18)

# Aplicamos el método de newton hasta que nuestro error sea menor que la␣
,→precisión

eps=1e-6 # Precisión deseada

it = 0 # Número de iteracciones
x1=2 # Estimación inicial. Sabemos que la función es creciente.
# Empezamos en el límite superior del intervalo

while f(x1)>eps:
x2 =x1-f(x1)/fp(x1)
[Link](x1,f(x1),'ko')
[Link](x2,0,'ko')
[Link]([x1,x1],[0,f(x1)],color='k',ls=':')
[Link](x1,-.5, r'$x_%d$' %it,ha='center')
[Link]([x1,x2], [f(x1), 0], 'k-')
x1=x2
it+=1

print("Hemos necesitado",it,"iteracciones.")
[Link](x1, f(x1), 'r*', markersize=15)

59
 [Link]("Raíz aproximada en %.3f" %x1,fontsize=14, family="serif",␣
,→xy=(x1,f(x1)), xycoords='data',

xytext=(-150, +50), textcoords='offset points',arrowprops=dict(arrowstyle="->",␣

,→connectionstyle="arc3, rad=-.5"))

ax.set_title("Método de Newton-Raphson")
fig.tight_layout()
show()

Hemos necesitado 5 iteracciones.

3.3.4. Ejercicio 5.7: función arcotangente hiperbólica

Utilizar el método de Newton para calcular la arcotangente hiperbólica de un número u:

x = arctanh u =⇒ u = tanh x.

Representar la arctanh en el intervalo (-1,1).

[ ]:

3.3.5. Ejercicio 5.8: las raíces de un polinomio

Dado el polinomio de grado 6

P( x ) = 924 x6 − 2772 x5 + 3150x4 − 1680 x3 + 420 x2 − 42 x + 1,

1. Representarlo en el intervalo (0,1) y encontrar aproximadamente la posición de las 6 raíces.
2. Escribir un programa que las calcule con una precisión mayor de 10−12 usando el método de
Newton.
[ ]:

60
3.3.6. Método de la secante
Una de las limitaciones del método de Newton-Raphson es que debemos conocer la derivada
de f ( x ).
Si no es el caso, aún la podemos estimar numéricamente usando las técnicas estudiadas en
el Tema 4,
que es lo que llamamos el método de la secante.
Para calcularla, necesitamos como mínimo puntos x1 y x2 que a diferencia del método de la
bisección, no tienen que contener a la raíz.
Con ellos calculamos la aproximación de la derivada:

f ( x2 ) − f ( x1 )
f 0 ( x2 ) = .
x2 − x1

lo cual nos permite obtener una nueva estimación de la raíz:

( x2 − x1 )
x3 = x2 − f ( x2 ) ,
f ( x2 ) − f ( x1 )

que como se puede ver depende de las dos estimaciones anteriores.

Al igual que el método de Newton-Raphson, converge muy rápidamente a la solución, y
tiene sus misma limitaciones.

3.3.7. Método de Newton-Raphson para más de una variable

Sea un conjunto de N ecuaciones no-lineales para N variables x1 · · · , x N

f 1 ( x1 , · · · , x N ) = 0 ,
f 2 ( x1 , · · · , x N ) = 0 ,
.. ..
. .
f N ( x1 , · · · , x N ) = 0 ,
(12)

cuyas raíces con x̄1 , · · · , x̄ N .

Expandiendo cada una de las ecuaciones en serie de Taylor:

∂ fi
f i ( x̄1 , · · · , x̄ N ) = f i ( x1 , · · · , x N ) + ∑ x̄ j − x j +··· ,
j
∂x j

que en notación vectorial se puede expresar como

F ( x̄ ) = F ( x ) + J · ( X̄ − X ) + · · · ,

61
 siendo J la matriz jacobiana del sistema de ecuaciones, es decir la matriz N × N con elemen-
∂f
tos Jij = ∂xij .

Puesto que el vector X̄ es solución del sistema, se tiene que F ( x̄ ) = 0, y

J · ( X − X̄ ) = J · ∆X = F ( x ),

es exactamente la generalización del método de Newton-Raphson para más de una variable,

que se puede resolver con los métodos que hemos visto para sistemas de ecuaciones lineales.
Una vez que hemos obtenido ∆X, podemos obtener la nueva estimación a partir de:

X 0 = X − ∆X.

En el caso de que la matriz Jacobiana no sea conocida, se puede estimar por medio del
método de las derivadas finitas, que generaliza por tanto el método de la secante para más
de una variable.

Ejemplo 5.19: representación gráfica del método de Newton-Raphson en varias variables. Re-
solver el sistema de ecuaciones:

y − x3 − 2 x2 + 1 = 0,
y + x2 − 1 = 0,
(13)

usando el método de Newton-Raphson en varias variables.

[1]: from numpy import linspace

from [Link] import subplots,plot,show,axhline

# definimos nuestro sistema de ecuaciones

def f(x):
return [x[1] - x[0]**3 - 2 * x[0]**2 + 1, x[1] + x[0]**2 - 1]

# definimos el jacobiano del sistema

def fJ(x):
return [[- 3*x[0]**2 - 4*x[0],1],[ 2*x[0],1]]

# representamos nuestra dos curvas

x=linspace(-3, 2, 5000)
y1=x**3+2*x**2-1
y2=-x**2+1

62
 fig,ax=subplots(figsize=(8,4))
[Link](x,y1,'b',lw=1.5, label=r'$y = x^3 + 2x^2 - 1$')
[Link](x,y2,'g',lw=1.5, label=r'$y = -x^2 + 1$')
[Link](loc=0)
axhline(0,ls=':',color="k")
show()

[2]: from numpy import array

# Obtenemos de ahí las estimaciones inciales:

x0=array([[-2.5, -7],[-1, 2],[1, 2]],float)

fig,ax=subplots(figsize=(8,4))
[Link](x,y1,'b',lw=1.5, label=r'$y = x^3 + 2x^2 - 1$')
[Link](x,y2,'g',lw=1.5, label=r'$y = -x^2 + 1$')
axhline(0,ls=':',color="k")

for i in x0:
[Link](i[0],i[1],'ko',markersize=10)
show()

63
[3]: # Aplicamos el método de Newton Rapshon
from numpy import copy
from [Link] import solve

eps=1e-6
x1=copy(x0)
it=0
for i in x1:
err=1
while err>eps:
dx=solve(fJ(i),f(i))
i-=dx
err=abs(max(dx))
x1[it]=i
it+=1

[4]: # Lo representamos

fig,ax=subplots(figsize=(8,4))
[Link](x,y1,'b',lw=1.5, label=r'$y = x^3 + 2x^2 - 1$')
[Link](x,y2,'g',lw=1.5, label=r'$y = -x^2 + 1$')
axhline(0,ls=':',color="k")

for k in range(3):

[Link](x1[k,0],x1[k,1],'r*', markersize=15)
[Link](x0[k,0],x0[k,1], 'ko')
[Link]("", xy=(x1[k,0],x1[k,1]), xytext=(x0[k,0],x0[k,1]),

64
 arrowprops=dict(arrowstyle="->", linewidth=2.5))

[Link](loc=0)
ax.set_xlabel(r'$x$',fontsize=18)
fig.tight_layout()
show()

Obviamente el método de Newton-Raphson para sistemas de ecuaciones no lineales ya esta

implementado.
Consiste en la función optimize del módulo scipy

[5]: from scipy import optimize

sol1=[Link](f,x0[0])
sol2=[Link](f,x0[1])
sol3=[Link](f,x0[2])

print(sol1,x1[0])
print(sol2,x1[1])
print(sol3,x1[2])

[-2.73205081 -6.46410162] [-2.73205081 -6.46410162]

[-1.00000000e+00 -8.07793567e-28] [-1. 0.]
[0.73205081 0.46410162] [0.73205081 0.46410162]
Vamos a aplicarlo para ver que solución obtenemos dependiendo de las condiciones inicia-
les.

65
[9]: fig,ax=subplots(figsize=(8, 4))
[Link](x, y1, 'k', lw=1.5, label=r'$y = x^3 + 2x^2 - 1$')
[Link](x, y2, 'k', lw=1.5, label=r'$y = -x^2 + 1$')

colors = ['r', 'b', 'g']

for m in linspace(-4, 3, 80):

for n in linspace(-15, 15, 40):

x00=[m, n]
sol=[Link](f,x00)

for idx,s in enumerate([sol1, sol2, sol3]):

if abs(s-sol).max() < 1e-5:
[Link](sol[0], sol[1], colors[idx]+'*', markersize=20)
[Link](x00[0], x00[1], colors[idx]+'.')

ax.set_xlabel(r'$x$',fontsize=18)
fig.tight_layout()
show()

Diferentes puntos iniciales en general llevan a diferentes soluciones

Ejercicio 5.9: sistema de circuitos no lineales Mientras que las resistencias introducen un efecto
lineal en la intensidad de corriente respecto a la diferencia de potencial (Ley de Ohm),

I = ∆V/R,

los diodos obedecen una ecuación no lineal:

66
  
I = I0 e∆V/VT − 1 ,

siendo V la diferencia de voltaje entre los extremos del diodo e I0 , V0 parámetros dependientes de
su composición.
Dado el siguiente circuito:

[4]: from [Link] import Image,display

display(Image("[Link]",width=600,height=600))

1. Aplicando la ley de conservación de la corriente de Kirchhoff (dado cualquier nodo, la co-

rriente que entra debe ser igual a la que sale), obtener las ecuaciones que han de cumplir los
potenciales V1 y V2 .
2. Usando el método de Newton-Raphson obtener el valor de V1 e V2 usando los valores:

67
 V+ = 5 V, R1 = 1 kΩ, R2 = 4 kΩ, R3 = 3 kΩ, R4 = 2 kΩ, I0 = 3 nA, VT = 0,005 V.

3. Comprobar que la diferencia de voltaje entre el diodo es aproximadamente 0.6 V, que es la

regla general aplicada por los ingenieros electrónicos.

[ ]:

4. Máximos y Mínimos de Funciones: minimización

La búsqueda de mínimos aparece en Física continuamente:
• soluciones de equilibrio,
• mínimos del potencial,
• problemas variacionales,
• estado fundamental de la función de ondas de un sistema físico en Mecánica Cuántica,
• sistemas de ecuaciones con más datos que parámetros.
• cálculo de cotas inferiores o superiores.
Minimización de funciones de una variable f ( x ) o de varias variables f ( x1 , x2 , · · · ).
Mínimos globales: el mínimo de la función en todo su dominio.
Mínimos locales: un punto en el que la función toma un valor más pequeño que cualquiera
de su entorno (posiciones de equilibrio inestable).
Mientras que en general sólo suele haber un mínimo global, puede haber varios mínimos
locales.
Los métodos numéricos no nos dicen si es local o global.

[2]: from [Link] import Image,display

display(Image("[Link]",width=600,height=600))

68
4.1. Estrategia general
Si es posible calcular la derivada:
• Buscar los puntos estacionarios en los que se anulen la derivada o las derivadas parcia-
les.
Si estamos buscando el mínimo de una función: f ( x1 , x2 , · · · ) resolvemos el sistema:

∂ f ( x1 , x2 · · · )
= 0, j = 1, 2, · · ·
∂x j

• Si el sistema de ecuaciones es lineal aplicamos la factorización LU o los métodos

iterativos que hemos visto en la sección 5.1.

• Si es un sistema no-lineal podemos utilizar el método del punto fijo/relajación o el mé-

todo de Newton-Raphson.
La aplicación del método de Newton-Raphson a la búsqueda de mínimos se llama mé-
todo de Gauss-Newton.
Si no es posible calcular la derivada, es necesario buscar el mínimo con métodos numéricos.

69
 4.2. Método de la razón áurea
Es análoga al método de la bisección pero para la búsqueda de mínimos.
Sólo sirve para funciones de una variable: buscamos los mínimos de f ( x ).
Al igual que el método de la bisección, requiere partir de dos valores de la variable x que
contengan el mínimo.
La situación que consideramos es la siguiente:

[3]: display(Image("[Link]",width=600,height=600))

1. Partiendo de las estimaciones iniciales x1 y x4 , asumimos que el mínimo de la función f ( x )

esta contenido en el intervalo [ x1 , x4 ].
2. Seleccionamos dos puntos x2 y x3 interiores en el intervalo [ x1 , x4 ].
3. Si la función en al menos uno de ellos f ( x2 ) o f ( x3 ) es menor que en f ( x1 ) y f ( x4 ), entonces
podemos asumir que como mínimo hay un mínimo de la función (local o global).
Si f ( x2 ) < f ( x3 ) entonces el mínimo estará en el intervalo [ x1 , x3 ].
Si f ( x2 ) > f ( x3 ) entonces el mínimo estará en el intervalo [ x2 , x4 ].
Esto nos permite reducir el intervalo en el que se encuentra el mínimo y eligiendo un punto
adicional, podemos repetir el proceso iterativamente.

70
 4. El problema reside en optimizar la posición de los puntos x2 y x3 .
Optimizando la posición de los puntos interiores
1. Como no sabemos en que lado del intervalo [ x1 , x4 ] se encuentra el mínimo, x2 y x3 deben
estar distribuidos simétricamente:

x2 − x1 = x4 − x3 .

2. Si x2 y x3 se encuentran muy cerca del centro del intervalo [ x1 , x4 ], uno de los intervalos de
la siguiente iteración, para el segmento [ x1 , x3 ] o [ x2 , x4 ] será mucho más grande que el otro.
Definimos el cociente entre los intervalos en dos iteraciones consecutivas:

x4 − x1 x2 − x1 + x3 − x1 x2 − x1
z= = = + 1, si f ( x2 ) < f ( x3 ),
x3 − x1 x3 − x1 x3 − x1
x4 − x1 x − x2 + x4 − x3 x − x3
z= = 4 = 4 + 1, si f ( x2 ) > f ( x3 ),
x4 − x2 x4 − x2 x4 − x2
(14)

3. Volviendo a definir el mismo cociente en la siguiente iteración (asumiendo que la función en

el nuevo punto a introducir es más pequeño):

x3 − x1
z0 = , si f ( x2 ) < f ( x3 ),
x2 − x1
x4 − x2
z0 = , si f ( x2 ) > f ( x3 ),
x4 − x3
(15)

Si queremos que el algoritmo sea igual de eficiente en dos iteraciones sucesivas entonces,
debemos imponer que z = z0 que implica que

z2 − z − 1 = 0,

como z>1 sólo podemos coger la solución positiva que implica que:

√
1+ 5
z= = 1,618...
2

que es la razón aúrea.

4. Una vez que tenemos el valor de z podemos obtener unívocamente el valor de los puntos
interiores:

x4 − x1 x4 − x1
x2 = x4 − , x3 = x1 +
z z

71
4.2.1. Algoritmo del método de la razón áurea
1. Escogemos una precisión para la obtención del mínimo.
2. Escogemos dos puntos en los extremos del intervalo en el que vamos a buscar el mínimo x1
y x4 .
Calculamos los puntos del interior x2 y x3 con la regla da la razón aúrea.
Evaluamos f ( x ) en los cuatro puntos y comprobamos que al menos uno de los valores
de la función en x2 o x3 es menor que en x1 y x4 .
3. Si f ( x2 ) < f ( x3 ) entonces el mínimo esta entre x1 y x3 , en este caso x3 se transforma en el
nuevo x4 , x2 se transforma en el nuevo x3 y habrá un nuevo x2 , obtenido por la regla de la
razón aúrea. Evaluamos f ( x ) en este nuevo punto.
4. En caso contrario, el mínimo esta entre x2 y x4 , entonces x2 se convierte en el nuevo x1 , x3
se convierte en el nuevo x2 , y habrá un nuevo x3 elegido con la regla de la razón aúrea.
Evaluamos f ( x ) en este nuevo punto.
5. Si x4 − x1 es mayor que la precisión requerida, repetimos desde el paso 3.
En caso contrario, calculamos 21 ( x2 + x3 ) y esta será nuestra estimación final de la posición
del mínimo.

4.2.2. Limitaciones del método de la razón áurea

Sufre problemas similares a los del método de la bisección:
Necesitamos conocer un intervalo que contenga el mínimo.
Es necesario conocer/hacer uso de las propiedades generales de la función permite estimar
este intervalo.
No se puede generalizar a funciones de varias variables.

4.2.3. Ejemplo 5.20: el potencial de Buckingham

El potencial de Buckingham es una representación aproximada de la energía de interacción entre
los átomos de un sólido ó gas como función de la distancia r entre ellos:
  
σ 6 −r/σ
V (r ) = V0 −e ,
r

El potencial contiene dos términos:

1. Uno positivo correspondiente a una fuerza repulsiva de corta distancia.
2. Otro negativo, que describe la fuerza atractiva de larga distancia.
Eso implica que hay un punto intermedio en el que el potencial tendrá un mínimo. El valor de r
asociado a ese mínimo corresponderá por tanto a la distancia en reposo entre los átomos del sólido
o gas. El valor de este mínimo sólo se puede obtener numéricamente.
Usar el método de la razón áurea para encontrar el mínimo del potencial de Buckingham para un
sólido con σ = 1 nm.

72
[1]: # 1- Constantes del problema

from numpy import exp,sqrt,linspace

from [Link] import plot,show,ylim

sigma=1.0 # Valor de sigma en nm

eps=1e-6 # precisión objetivo
z = (1+sqrt(5))/2 # razón áurea

# 2. definimos nuestro potencial normalizado (el mínimo no depende de V0)

def f(r):
return (sigma/r)**6-exp(-r/sigma)

# 3. Lo representamos para estimar el intervalo

x = linspace(0.1,10,1000)
ylim(-0.2,0.1)
plot(x,f(x))
show()

[4]: # 4. Posición incial de nuestros puntos, y valor de la función en ellos

x1=sigma/10
x4=sigma*10

73
x2=x4-(x4-x1)/z
x3=x1+(x4-x1)/z

f1=f(x1)
f2=f(x2)
f3=f(x3)
f4=f(x4)

# 5. Aplicamos el algoritmo de la razón áurea para acortar el intervalo

while x4-x1>eps:
if f2<f3:
x4,f4=x3,f3
x3,f3=x2,f2
x2=x4-(x4-x1)/z
f2=f(x2)
else:
x1,f1=x2,f2
x2,f2=x3,f3
x3=x1+(x4-x1)/z
f3=f(x3)

xmin=(x1+x4)/2
print("El mínimo cae en",0.5*(x1+x4),"nm")

plot(x,f(x))
plot(xmin,f(xmin),"ko")
ylim(-0.2,0.1)
show()

El mínimo cae en 1.630516067174875 nm

74
4.2.4. Ejercicio 5.10: la temperatura de una bombilla
Tradicionalmente, una bombilla esta formada por un filamento de tungsteno, que se calienta cuan-
do la corriente eléctrica pasa por él, emitiendo radiación térmica. Esencialmente toda la energía
consumida es radiada electromagnéticamente, pero una parte de esa radiación no es emitida en el
espectro visible.
Definamos la eficiencia de una bombilla como la energía radiada en el visible entre la total. Asu-
miendo que un filamento de tungsteno cumple la ley de Planck, la potencia radiada por unidad
de longitud de onda viene dada por:

λ −5
I (λ) = 2πAhc2 ,
ehc/λk B T − 1

siendo A la superficie del área del filamento de tungsteno, h la constante de Planck, c la velocidad
de la luz, k B la constante de Boltzmann.
Teniendo en cuenta que radiación visible se encuentra entre las longitudes de onda λ1 = 390 nm
y λ2 = 750 nm.
hc
1. Llamando x = x (λ) = , y xi = x (λi ) con i = 1, 2, probar que la eficiencia de una
λk B T
bombilla de tungsteno η viene dada por la integral:

R x2
x
x3 /(ex − 1) dx 15
Z x2
x3
η = R ∞1 = dx,
0
x3 /(ex − 1) dx π4 x1 (e x − 1)

75
 donde en el último paso hemos hecho uso de que la integral del denominador se conoce
analíticamente.
2. Escribir un programa que tome la temperatura como argumento y que calcule el valor de η
para dicho valor.
Como la integral no se puede resolver analíticamente, usar el método más conveniente para
integrarla numéricamente.
Representar la eficiencia de la bombilla para temperaturas entre 300 K y 10000 K.
Comprobar que hay una temperatura intermedia para la que la eficiencia es máxima.
3. Calcular la temperatura que corresponde al máximo de la eficiencia con una precisión de 1
K usando la razón áurea.
4. ¿Tiene sentido practico calentar el tungsteno a la temperatura del máximo?.

[ ]:

4.3. El método de Gauss-Newton

El método de la razón áurea es rápido y robusto, sin embargo:
• Requiere empezar por un intervalo que contenga al máximo.
• No se puede generalizar a más de una variable.
Sin embargo un método de minimización consiste en encontrar solución de la ecuación:

f 0 ( x ) = 0,

en el caso de una variable, y de forma similar para más de una variable.

Solución que se puede encontrar usando el método de Newton-Raphson, que aplicado a la
derivada se llama método de Gauss-Newton.
Dada una estimación inicial x la siguiente vendrá dada por:

f 0 (x)
x0 = x − ,
f 00 ( x )

Si conocemos la función y sus derivadas, el método nos permite encontrar los máximos y
mínimos de una función de forma rápida y fiable.
El método de Gauss-Newton converge de forma cuadrática y puede generalizarse a varias
variables.
Sin embargo es bastante común que la primera y/o la segunda derivada no sean conocidas.
En ese caso estás se pueden aproximar usando los métodos de la diferencias finitas del Tema
3.
Da lugar a un método de la secante para un problema de minimización.

76
4.4. El método del gradiente descendente
El método del gradiente descendente es una modificación del método de Gauss-Newton
para el caso en que la segunda derivada no se puede calcular ni aproximar correctamente.
En ese caso, se aproxima por una constante γ que debe tener el orden de magnitud correcto
de forma que:

x 0 = x − γ f 0 ( x ).

El valor de γ no tiene que estar determinado muy precisamente, es suficiente con que el
orden de magnitud sea correcto.
Si γ tiene signo positivo el método le sustraerá el valor de la derivada γ veces a nuestra
estimación x, y el método convergerá hacía un mínimo.
Si γ tiene signo negativo, el método le añadirá dicho valor, y convergerá a un máximo.
Si γ es grande podemos converger más rápidamente, pero si lo es demasiado pasaremos por
encima del mínimo y el método fallará.
El valor de γ tiene que ser optimizado para maximizar la convergencia sin saltarnos en
ningún momento el máximo.
Para ello, el método más rápido es probar con varios valores hasta encontrar el valor más
óptimo.
Si la primera derivada tampoco se puede calcular, se estima una vez más numéricamente:

f ( x2 ) − f ( x1 )
f 0 (x) ' ,
x2 − x1

por lo que la siguiente estimación vendrá dada por

f ( x2 ) − f ( x1 )
x3 = x2 − γ
x2 − x1

4.5. El método del gradiente descendente en varias dimensiones

Los algoritmos de búsqueda de mínimos (optimización) en varias dimensiones usan una
misma idea básica:
1. Dada una estimación inicial X0 = x01 , x02 , · ·
· , encontramos una dirección D0 a lo largo

de la cual el valor de f ( X0 ) = f x01 , x02 , · · · decrece.

2. Después, buscamos el menor valor de f ( X ) en esa dirección empezando en X0 .
Esto hace que el problema sea equivalente a uno unidimensional:
Hay que minimizar la función:

f ( X0 + α D0 ) , α > 0,

77
 Esta optimización unidimensional se denomina búsqueda lineal y el valor de α que
minimiza f ( X0 + α D0 ), α0 , se denomina longitud de paso.
3. Una vez se ha encontrado α0 tenemos una nueva estimación:

X1 = X0 + α0 D0

el proceso se repite de nuevo empezando por el punto 1 generando una secuencia de itera-
ciones:
$$X_{k+1}= X_{k} +\alpha_k D_k.$$
Los diferentes métodos para encontrar mínimos se diferencian en como escoger las direccio-
nes Dk .
La virtud de estos métodos es que remplazan el problema n-dimensional inicial por una
secuencia de problemas unidimensionales, más sencillos.

4.5.1. Método de máxima pendiente

El problema por tanto se resume en encontrar la dirección en la que la función decrece.
Si podemos calcular las derivadas parciales de la función, la dirección en que la función varia
más rápidamente viene dada por el gradiente.
De esa forma el método consiste en tomar:

Dk = − ∇ f ( X )| X =Xk ,

4.5.2. Método del gradiente conjugado

Es una mejora del método de máxima pendiente. En esta caso la dirección de cambio se elige
como:

Dk = − ∇ f ( X )| X =Xk + β k Dk−1 ,

donde:

∇ f ( X )T ∇ f (X ) Xk
βk = ,
∇ f (X )T ∇ f (X ) | X k −1

De esta forma las direcciones elegidas en cada iteración son conjugadas con respecto a la
matriz Hessiana de f ( x ):

DkT · H ( X )| XK · Dk−1 = 0.

Que implica que lo que bajamos en una dirección no lo deshacemos en las siguientes.

78
[17]: from [Link] import display,Image
display(Image("Conjugate_gradient.png",width=400))

79
80
 El paquete scipy en su módulo optimize proporciona la función minimize donde se pueden esco-
ger diferentes método más avanzados:
Método de Brendt.
Método de Nelder–Mead.
Método de Powell.
Método BFGS.
[Link]

5. Ejercicio adicionales
5.1. Ejercicio 5.11: minimizando el potencial de ionización (Examen de Enero de 2023)
El potencial de ionización utilizado en un experimento para atrapar átomos viene dado por la
integral:

Z     
1 1 1
Vion ≡ Vion (λ) = α cos2 + λ2 sin2 +y dx dy
D λ2 2xy − 1 xy

siendo α = 10−3 eV, D el conjunto de puntos pertenecientes a R2 tal que

D ≡ ( x, y) x2 + y2 ≤ 1

y λ un parámetro libre adimensional experimental que se puede ajustar.

El parámetro λ se puede ajustar experimentalmente para garantizar que el valor del potencial de
ionización sea lo menor posible.
Encontrar el valor de λ que hace mínimo Vion .

[ ]:

5.2. Ejercicio 5.12: funciones complejas con polos (Examen de Junio de 2023)
Dada la función

g2
t(s) = , M = 0,75 GeV, g = 0,35 GeV,
M 2 − s + g2 J ( s )

definida para valores de s ∈ C y con J(s) dada por

2 σ(s) σ(s) − 1
J (s) = + log ,
π π σ(s) + 1

donde

81
 r
sth
σ(s) = 1− , con sth = 4m2 , m = 0,140 GeV.
s

Encontrar el valor de s ∈ C para el que la funcion t(s) tiene un polo.

[ ]:

5.3. Ejercicio 5.13: cuadrando posiciones

Un barco enemigo está atacando la costa Este de la isla de las pirulétas, sin embargo, dada la
distancia de la nave enemiga no somos capaces de determinar su posición.
No obstante, a partir de las posiciones donde han caído los proyectiles queremos determinar desde
donde nos atacan. El [Link], en la carpeta de datos del tema, guarda en dos columnas
las coordenadas ( x, y) de los cráteres dejados por los proyectiles, medidos respecto de un origen
O tal y como se muestra en el siguiente esquema (los cráteres tienen un diámetro aproximado de
10 m) .
[15]: from [Link] import Image,display
display(Image("[Link]",width=400,height=400))

Por el tipo de proyectil sabemos que se trata de un cañón que dispara a una velocidad de 320 m/s
y que el ángulo con la horizontal al realizar el disparo es tal que se maximiza el alcance. Bajo estas
circunstancias se tiene que los proyectiles deben caer en los puntos que cumplen

82
 ( x − d )2 + ( y − c )2 = R2 ,
donde R es la distancia máxima alcanzada.
1. Considerando que la altura de lanzamiento y la de impacto son la misma, obtener la expre-
sión de R. Usando esta definir una función de python tal que represente el modelo x̃ (y|d, c),
donde d y c son los parámetros del modelo.
2. A partir del modelo, uno puede construir la función χ2 , definida como

[ xi − x̃ (yi |d, c)]2

χ2 (d, c) = ∑ σi
,
i

donde el índice i hace referencia a los datos del archivo [Link]. Definir una función de
python que tome como entrada, por lo menos, los parámetros del modelo, es decir, d y c y
calcule el valor de χ2 (d, c).
3. Leer los datos contenidos en el fichero [Link] y representar en el plano los puntos de
impacto.
4. Construir una matriz que contenga el valor de χ2 (d, c) para distintos valores de los pará-
metros y representar mediante la función imshow. Mediante un proceso manual, ajustar los
límites del cuadrado [d0 , d1 ] × [c0 , c1 ], para localizar el mínimo de la función χ2 (d, c) de forma
aproximada.
AYUDA: limitar el valor máximo del colorbar puede ser útil para determinar la región del
mínimo.
Para tratar de localizar el mínimo vamos a usar el método del gradiente, para lo cual necesitamos
calcular el gradiente de la función χ2 (d, c).
5. Definir una función que determine el valor del gradiente de la función χ2 (d, c) en un punto
concreto y con un paso espacial determinado. Emplear un método de diferencias finitas para
determinar el valor del gradiente.
El método del gradiente para buscar mínimos de una función resulta un proceso iterativo,
donde a cada paso actualizamos el valor de los parámetros, de forma que:

~θi+1 = ~θi − γ∇χ2 (~θi ) ,

siendo γ el ratio de aprendizaje. Dicho proceso iterativo se suele detener cuando la diferencia
en los parámetros entre una iteración y la siguiente es menor que una cierta tolerancia.
6. Tomando γ ∼ 0,01 definir una función que tome un punto aleatorio del cuadrado [d0 , d1 ] ×
[c0 , c1 ] y realice el proceso iterativo para encontrar el mínimo de la función χ2 (d, c).
7. Usando la función antes definida, repetir el proceso de búsqueda del mínimo para varios
puntos iniciales del cuadrado [d0 , d1 ] × [c0 , c1 ] y guardar los valores de los parámetros d y c
obtenidos cada vez.
8. Mediante el estudio estadístico de los resultados obtenidos determinar el valor de d y c, así
como el error.

83
[ ]:

5.4. Ejercicio 5.14: teoría del campo medio para el ferromagnetismo

En la teoría de campo medio del ferromagnetismo, la intensidad de magnetización M de un ma-
terial ferromagnético depende da la temperatura según la ecuación:

JM
M = µ tanh ,
kB T

donde µ es el momento magnético, J es una constante de acoplo y k B es la constante de Boltzmann.

Definiendo m = M/µ y C = µJ/k B la magnetización por unidad de momento magnético viene
dada:

Cm
m = tanh .
T

Cogiendo C = 1, encontrar la solución de la ecuación más allá de la trivial, m=0, para temperaturas
0 < T < 2 K con una tolerancia para la solución e = 10−6 .

[ ]:

5.5. Ejercicio 5.15: potencial cuadrado cuántico

Considerar un potencial cuadrado de anchura L y profundidad V

[11]: from [Link] import Image,display

display(Image(filename="Finite_Square_Potential.png",width=600,height=600))

84
 1. Usando las ecuaciones de Schrödinger probar que las energías permitidas de una partícula
de masa m atrapada en el pozo son las soluciones de

q  p
tan L2 2
m E/2h̄ = p(V − E) / E para un número par de estados,
− E / (V − E) para un número impar de estados,

donde los estados están numerados empezando por 0, el estado fundamental, 1, el primer
estado excitado, etc.
2. Para un electrón, me = 9,1094 × 10−23 kg, en un pozo de profundidad V = 20 eV y anchura
L = 1 nm, representar en una misma gráfica las funciones

q q q
y1 = tan L2 m E/2h̄2 , y2 = (V − E) / E, y3 = − E / (V − E )

como función de E en el intervalo (0, 20) eV. A partir de la gráfica estimar el valor de los 6
primeros niveles de energía.
3. Escribir un programa que calcule los 6 primeros niveles de energía con una precisión de 10−6
usando el método de la bisección.
[ ]:

85
 5.6. Ejercicio 5.16: el punto de Lagrange Tierra-Luna
El punto de Lagrange del sistema Tierra-Luna esta definido como aquel que siempre se encuentra
a la misma distancia entre ambos.
Es decir, aquel punto en el que un satélite orbitaria en torno a la Tierra en perfecta sincronía con la
Luna.
Esto se debe a que la atracción gravitatoria en ese punto debido a la Tierra y la Luna se combinan
para crear la fuerza centrípeta que mantiene al satélite en una órbita constante.

[70]: from [Link] import Image,display

display(Image("[Link]",width=600,height=600))

1. Asumiendo una orbita circular, probar que la distancia r desde el centro de la tierra al punto
de Lagrange L satisface:

GM Gm
− = ω 2 r,
r 2
( R − r )2

siendo R la distancia entre la Tierra y la Luna, M la masa de la Tierra, m la de la Luna, G la

constante de la gravitación,
y ω la velocidad angular del satélite (que tiene que ser igual a la de la Luna).
2. La ecuación anterior es un polinomio de grado 5. Escribir un programa que use el método
de Newton y el de la secante para obtener el valor de r con una precisión de 10−6 .
Los valores de los diferentes parámetros del problema son:
G = 6,674 × 10−11 N m2 Kg−1 ,
M = 5,974 × 1024 Kg,
m = 7,348 × 1022 Kg,
R = 3,844 × 108 m,
v = 1,023 Km/s,
siendo v la velocidad con la que la Luna se desplaza en su orbita.

86
[ ]:

5.7. Ejercicio 5.17: ajustando el fondo cósmico de microondas

El satélite Cosmic Background Explorer (COBE) de la NASA llevaba un instrumento, FIRAS (Far-
Infrared Absolute Spectrophotometer) para medir la radiación del fondo cósmico de microondas
(CMB), que se confirmó que se distribuía según una curva de cuerpo negro de acuerdo con la
teoría del big bang:

2 h ν3 c2
I (ν, T ) =  
exp khBc Tν − 1

donde la frecuencia de radiación ν se expresa en números de onda, cm−1 , y la velocidad de la luz,

c, se toma en cm s−1 .
El archivo de datos “[Link]” en la carpeta de datos del tema 5 contiene mediciones de
I (ν, T ) basadas en las observaciones de FIRAS. Nótese que las unidades de I en este archivo son
ergs−1 cm−1 sr −1 , siendo sr el esteradian, la unidad de ángulo sólido, y teniendo en cuenta que
1J ≡ 107 erg.
Determinar la temperatura del CMB minimizando las mediciones con la distribución del cuerpo
negro.

[ ]:

Final del Tema 5

87