Métodos Numéricos Problemas Resueltos y Prácticas
Métodos Numéricos Problemas Resueltos y Prácticas
Métodos Numéricos Problemas Resueltos y Prácticas
numéricos
Problemas
resueltos
y prácticas
Isaac A. García
Susanna Maza
Subido por:
https://www.facebook.com/pages/Interfase-
IQ/146073555478947?ref=bookmarks
Métodos numéricos
ISBN: 978-84-8409-
Maquetación
Edicions i Publicacions de la UdL
Diseño de portada
cat & cas
La reproducción total o parcial de esta obra por cualquier procedimiento, incluidos la reprografía y el tratamiento
informático, y la distribución de ejemplares mediante alquiler o préstamo público, queda rigurosamente prohibida
sin la autorización de los titulares del “copy-right”, y será sometida a las sanciones establecidas por la Ley.
Índice general
1. Prólogo 11
3. Sistemas Lineales 27
3.1. Resumen de teorı́a . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.1. Algoritmo de eliminación Gaussiana . . . . . . . . . . . . 27
3.2. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . . 31
4. Interpolación Polinómica 43
4.1. Resumen de teorı́a . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . . 49
5. Derivación Numérica 71
5.1. Resumen de teorı́a . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.2. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . . 72
6. Integración Numérica 85
6.1. Resumen de teorı́a . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.2. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . . 90
5
6 ÍNDICE GENERAL
Prólogo
11
12 Prólogo
Uno de los lenguajes utilizados en las prácticas será el Fortran (el nombre
que proviene de las primeras letras de las palabras inglesas FORmula y TRANs-
lation) que es un lenguaje de traducción de fórmulas, ver por ejemplo [3]. El
nombre sugiere la aplicación más importante de dicho lenguaje: calcular expre-
siones matemáticas. Es, por lo tanto, un lenguaje de programación muy usual
en cálculo numérico. Se utilizará también el programa de libre distribución Oc-
tave que es también muy adecuado para cálculos numéricos. Otra opción que se
utiliza consiste en usar también manipuladores algebraicos (en concreto Math-
ematica) que tanto se han aconsejado para otras asignaturas de matemáticas
de carreras técnicas, de manera que en muchas de las prácticas propuestas se
pide realizar programas que se podrı́an sustituir por la utilización de un único
comando de tales manipuladores. Pero, sin embargo, creemos que debe de ser
un objetivo, cuando se trata de asignaturas de carácter numérico, el utilizar
un lenguaje de programación y profundizar más en la propia estructura de los
métodos. Si hacemos esto es porque estamos convencidos, y la experiencia con
nuestros alumnos ası́ nos lo demuestra, de que sólo cuando uno se enfrenta a la
programación efectiva de los métodos es capaz de entenderlos en profundidad.
De cualquier forma, en el estudio de los métodos de aproximación puede ser
sumamente útil contar con un paquete de cómputo simbólico. Por ejemplo, mu-
chos métodos numéricos tienen lı́mites de error que exigen acotar una derivada
ordinaria o parcial de una función. A veces se trata de una tarea tediosa y de
poco valor didáctico una vez dominados los métodos del cálculo. Estas derivadas
pueden obtenerse rápidamente en forma simbólica, y un poco de ingenio per-
mitirá que una computación simbólica facilite también el proceso de acotación.
Hemos escogido Mathematica para este uso debido a que los estudiantes ya
lo han utilizado en prácticas realizadas en las asignaturas de Cálculo y Álgebra
Lineal de primer curso. Una guı́a completa sobre las posibilidades de Math-
ematica se encuentra en el manual del programa contenido en la página web
13
http://documents.wolfram.com/v5/TheMathematicaBook/index.html.
Errores, Estabilidad y
Condicionamiento
Sea x∗i un valor aproximado del valor xi . Denotamos el error absoluto por
Δxi := |xi − x∗i | y al error relativo de la forma xi := |Δxi /xi | siempre que
xi = 0.
Sea f : D ⊂ Rn → R con f ∈ C 1 (D) tal que y = f (x1 , x2 , . . . , xn ). Las
fórmulas de propagación de errores son:
n
Error Absoluto: Δy ≈ j=1 ∂f∂x(x) j
Δxj .
n xj ∂f (x)
Error Relativo: y ≈ j=1 f (x) ∂xj xj .
15
16 Errores, Estabilidad y Condicionamiento
Se sabe que, con el ordenador, no podemos trabajar con los números reales,
sino sólo con un subconjunto de los números racionales. Como existe un número
representable por el ordenador máximo en valor absoluto, cuando el ordenador
intenta utilizar uno de mayor valor absoluto se produce un error llamado over-
flow. También se tiene el caso contrario que consiste en utilizar un número no
nulo pero menor en va-lor absoluto que el mı́nimo número representable en
valor absoluto. Entonces el ordenador devuelve un error llamado underflow. Se
dispone del concepto epsilon de la máquina, que es el mayor valor positivo que
verifique
1 + x = 1 ∀x ∈ (0, ) .
Tenemos pues que la precisión de la aritmética en coma flotante se caracteriza
por y dependerá del ordenador y del tipo de variable.
∂g 4π 2 ∂g 8π 2
= 2 , =− 3 ,
∂ T ∂T T
se obtiene
g ≈ + 2T .
Según los datos del enunciado
0,001 0,01
= = 0,001 , T = = 0,005 ,
1 2
de manera que
g ≈ 0,001 + 2 × 0,005 = 0,011 .
2.2 Problemas resueltos 17
Solución. (i) Puesto que las únicas variables del problema que contienen
errores son μ1 y R, la fórmula de propagación de errores absolutos permite
aproximar el error absoluto Δv con el que obtiene la velocidad v de la forma
∂v
Δv = Δμ + ∂v ΔR ,
∂μ1 1 ∂R
Puesto que
T
v(T ) = 331 1+ ,
273
y su derivada es
dv 331 1
= ,
dT 273 2 1 + T
273
se tiene que
T
κT = .
2(273 + T )
2x + ay = 1,
x + by = 0,
con 2b − a = 0.
(i) Si a = 1 ± 0,1 y b = 2 ± 0,2, aproximar con qué exactitud se puede calcular
x + y.
(ii) Obtener el rango de valores de a y de b para los cuales el problema de hallar
x + y está mal condicionado.
b −1
x= , y= ,
2b − a 2b − a
de manera que, definimos su suma como la función
b−1
f (a, b) = x + y = .
2b − a
Aplicando la fórmula de propagación de errores absolutos, se tiene que el error
absoluto Δf al evaluar la suma x + y está realcionado con los errores absolutos
de a y de b mediante
∂f ∂f
Δf ≈ Δa + Δb .
∂a ∂b
Puesto que las derivadas parciales de f de primer orden son
∂f b−1 ∂f 2−a
= , = ,
∂a (a − 2b)2 ∂b (a − 2b)2
obtenemos
b−1
Δf ≈ Δa + 2 − a Δb .
(a − 2b)2 (a − 2b)
2
por lo tanto
a b(a − 2)
κa =
, κb = .
2b − a (a − 2b)(b − 1)
De aquı́ se deduce que, si a − 2b ≈ 0 o bien b − 1 ≈ 0 entonces el problema de
hallar x + y está mal condicionado.
n
|x |
= i i ,
n
i=1 j=1 xj
Por lo que se refiere a la pregunta ¿Es bueno desde el punto de vista numérico
restar dos números reales muy próximos? la respuesta es NO a la vista del
resultado anterior puesto que x1 + x2 ≈ 0 lo cual implica κxi >> 1.
despejamos D y obtenemos
1
D = f (k, T, θ) = ln cosh[T kg sin θ] .
k
El valor aproximado de D será
En definitiva se obtiene
D = D∗ ± ΔD ≈ 6,060 ± 1,212 .
n 2 6 10 25 550
xn 2 3.1365 3.1415 3.1425 0
Solución. Puesto que las únicas magnitudes que contienen errores son x y R,
2.2 Problemas resueltos 25
tomaremos la función
σ x
E(x, R) = 1− √ .
20 x + R2
2
Sistemas Lineales
Métodos directos, que son algoritmos finitos para el cálculo de las solu-
ciones de un sistema lineal. El método de eliminación Gaussiana es un
ejemplo de este tipo.
27
28 Sistemas Lineales
(3) Sumar a una ecuación una combinación lineal de las ecuaciones restante.
Es conocido que las transformaciones elementales no modifican la solución del
sistema lineal.
n
|aii | > |aij | para i = 1, . . . , n.
j=1
j=i
Esta última ecuación nos permite asegurar que la matriz original de coeficientes
del sistema puede descomponer en el producto de dos matrices triangulares de
la forma
⎛ ⎞
a11 a12 a13 · · · a1n
⎜ a21 a22 a23 · · · a2n ⎟
⎜ ⎟
⎜ a31 a32 a33 · · · a3n ⎟
⎜ ⎟= (3.8)
⎜ .. .. .. .. .. ⎟
⎝ . . . . . ⎠
an1 an2 an3 · · · ann
⎛ ⎞ ⎛ a a12 a13 · · · a1n
⎞
1 0 0 ··· 0 11
⎜ m21 ⎜ 0 a(1) a(1) · · · a2n ⎟
(1)
⎜ 1 0 ··· 0 ⎟ ⎟ ⎜ 22 23 ⎟
⎜ m31 m32 ⎜ ⎟
⎜ 1 ··· 0 ⎟ ⎟˙⎜ 0 0
(2)
a33 · · ·
(2)
a3n ⎟
⎜ .. . ⎟ ⎜ ⎟
. .. ⎠ ⎜ ⎟
.. .. .. . .. .. .. ..
⎝ . . . ⎝ .. . . . . ⎠
mn1 mn2 mn3 · · · 1 0 0 0
(n−1)
· · · ann
con > 0. Utilizando la norma del infinito y suponiendo que ≤ 0,01, demostrar
que una pequen̄a perturbación relativa del término independiente b puede pro-
ducir una perturbación relativa 40401 veces mayor en la solución del sistema
lineal.
3.2 Problemas resueltos 33
−1 2 3 −1
−1 2 3 −1 0 3 3 2
⎛ ⎞
1 1 0 3
⎜ 0 −1 −1 −5 ⎟
A(2) = ⎜
⎝ 0 0
⎟ = A(3) .
3 13 ⎠
0 0 0 −13
Solución. Recordemos que, dada una matriz A ∈ Mn (R), si todos sus de-
teminantes menores principales Δk con k = 1, 2, . . . , n, son diferentes de cero,
entonces existen matrices triangulares únicas L = (ij ) y U = (uij ) triangulares
inferiores y superiores respectı́vamente con ii = 1 para i = 1, 2, . . . , n tales que
A = LU .
1 2
Δ1 = 1 = 0 , Δ2 = = 0 , Δ3 = det A = 175 = 0 .
4 8
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
1 2 6 1 0 0 u11 u12 u13
⎝ 4 8 −1 ⎠ = ⎝ 21 1 0 ⎠ ⎝ 0 u22 u23 ⎠ .
−2 3 5 31 32 1 0 0 u33
1 =
u11 ,
4 = m21 u11 = m21 ,
−2 = m31 u11 = m31 .
2 = u12 ,
8 = m21 u12 + u22 = 4 × 2 + u22 , u22 = 0 ,
3 = m31 u12 + m32 u22 = −2 × 2 + m32 × 0 = −4 ,
(i) Averiguar para qué valores de los parámetros a y b es posible hallar la des-
composición LU de la matriz de coeficientes del sistema a partir del algo-
ritmo de eliminación Gaussiana sin pivotación. Para esos valores de los
parámetros obtener las matrices L y U . ¿Es única dicha descomposición
si a = 2 y b = 7?
(ii) Utilizando la norma L1 , calcular el número de condición de la matriz de
coeficientes del sistema cuando a = 1 y b = 3.
(iii) Si a = 2 ± 0,5 y b = 5 ± 0,1, aproximar la precisión con la que se puede
determinar x + y.
Δ1 = 7 = 0 , Δ2 = 12 = 0 ,
(ii) Si a = 1 y b = 3 la matriz A es
3 2
A= .
1 1
3.2 Problemas resueltos 37
La matriz inversa de A es
−1 1 2+ −2
A = .
−1 1
38 Sistemas Lineales
A1 = máx{2, 4 + } = 4 + ,
1 3+
A−1 1 = máx{3 + , 3} = .
Se concluye que, el número de condición κ(A) asociado a A en la norma L1 vale
3+
κ(A) = A1 A−1 1 = (4 + ) .
Para el caso = 10−5 se obtiene
(i) Supongamos que se han de hallar las posiciones de equilibrio x∗i con i =
1, 2, 3 de los osciladores en 1000 experimentos diferentes, donde la única
variación entre los experimentos es la longitud d entre las paredes. Ra-
zonar qué método se programarı́a en un ordenador para resolver los 1000
sistemas lineales asociados a cada experimento.
(ii) Si todos los muelles son iguales de constante elástica k y longitud natural ,
utilizar el método elegido en el apartado anterior para hallar las posiciones
de equilibrio de los osciladores.
(ii) Si todos los muelles son iguales de constante elástica k y longitud natural
, el sistema lineal a resolver es
⎛ ⎞⎛ ∗ ⎞ ⎛ ⎞
−2k k 0 x1 0
⎝ k −2k k ⎠ ⎝ x∗2 ⎠ = ⎝ 0 ⎠ .
0 k −2k x∗3 −kd
se tiene
⎛ ⎞ ⎛ ⎞
−2k k 0 −2k k 0
A = ⎝ k −2k k ⎠∼⎝ 0 −3/2k k ⎠
0 k −2k 0 k −2k
⎛ ⎞
−2k k 0
∼ ⎝ 0 −3/2k k ⎠=U .
0 0 −4/3k
Las transformaciones elementales por filas realizadas han sido en primer lugar
fi → fi − mi1 f1 para i = 2, 3, con m21 = −1/2 y m31 = 0 y en segundo lugar
f3 → f3 − m32 f2 con m32 = −2/3. De este modo se tiene que
⎛ ⎞
1 0 0
L = ⎝ −1/2 1 0 ⎠ .
0 −2/3 1
0 0 4 −23 −21
Permutamos ahora la fila 2 y la 3 con lo cual
⎛ ⎞
12 −8 6 −10 34
⎜ 0 −11 7,5 5,5 18,5 ⎟
A(2) ∼ A(3) = ⎜⎝ 0
⎟ .
2 −1 9 −5 ⎠
0 0 4 −23 −21
Interpolación Polinómica
Son muchos los problemas en ingenierı́a en los que hay que estudiar funciones
que no se conocen explı́citamente. En su lugar sólo se conoce una tabulación
de los valores que toma la función en ciertos puntos (por ejemplo como resulta-
do de medidas experimentales). Dichos problemas son buenos candidatos para
la utilización de las técnicas interpolatorias que explicaremos a continuación.
Nosotros sólo vamos a centrarnos en la interpolación polinómica que es la más
sencilla.
Teorema 4.1 Supongamos conocido el valor de una función f (x) en un con-
junto de puntos distintos dos a dos {x0 , x1 , . . . , xn }. Entonces, existe un único
polinomio Pn (x) de grado menor o igual que n que interpola a la función en
estos puntos, es decir Pn (xi ) = f (xi ) para i = 0, . . . , n.
43
44 Interpolación Polinómica
f [xi ] = fi , (4.3)
f [xi+1 , xi+2 , . . . , xi+j ] − f [xi , xi+1 , . . . , xi+j−1 ]
f [xi , xi+1 , . . . , xi+j ] = .
xi+j − xi
x0 f [x0 ]
f [x0 , x1 ]
x1 f [x1 ] f [x0 , x1 , x2 ]
f [x1 , x2 ] f [x0 , x1 , x2 , x3 ]
x2 f [x2 ] f [x1 , x2 , x3 ]
f [x2 , x3 ]
x3 f [x3 ]
Teorema 4.2 Sea f una función de clase C n+1 [a, b], y sea Pn un polinomio de
grado menor o igual que n que interpola a la función f en los siguientes n + 1
puntos {x0 , x1 , . . . , xn } ⊂ [a, b] distintos dos a dos. Entonces, para cualquier
x ∈ [a, b], existe un punto ξx ∈ (a, b) tal que
f (n+1) (ξx )
n
f (x) − Pn (x) = (x − xi ) . (4.4)
(n + 1)! i=0
Teorema 4.3 Sea f : [a, b] ⊂ R → R una función derivable en [a, b]. Con-
sideremos n + 1 puntos {x0 , x1 , . . . , xn } ⊂ [a, b] distintos dos a dos. Existe
un único polinomio H2n+1 (x) de grado menor o igual que 2n + 1 tal que
H2n+1 (xk ) = f (xk ) y H2n+1 (xk ) = f (xk ), k = 0, 1, . . . , n. La expresión ex-
plı́cita de H2n+1 (x) es la siguiente:
n
n
H2n+1 (x) = fk hk (x) + fk h̃k (x) ,
k=0 k=0
siendo hk (x) = [1 − 2k (xk )(x − xk )]2k (x), h̃k (x) = (x − xk )2k (x), donde k (x)
denota al polinomio k-ésimo básico de Lagrange. Además, si f ∈ C 2n+2 [a, b],
entonces
f (2n+2) (ξx )
n
f (x) − H2n+1 (x) = (x − xk )2 ,
(2n + 2)!
k=0
n
πn (x)∞ := máx |x − xi |
a≤x≤b
i=0
sea la más pequen̄a de entre cualquier otra posible elección de nodos. La respues-
ta a esta pregunta es afirmativa y su solución se basa en una teorı́a desarrollada
por el matemático ruso P.L. Chebychev que exponemos a continuación.
Tomando θ = arc cos x y sumando las dos anteriores ecuaciones se tiene que
es decir, las funciones gn (x) verifican la misma recurrencia que los Tn (x). Además,
satisfacen las dos condiciones iniciales de la recurrencia puesto que
Demostración. Por reducción al absurdo, supongamos que |p(x)| < 21−n para
todo x ∈ [−1, 1] y lleguemos a una contradicción. Sea q(x) = 21−n Tn (x) siendo
Tn (x) el polinomio de Chebychev de grado n. Entonces, q(x) es mónico y de
4.1 Resumen de teorı́a 47
Teorema 4.7 Sea f una función de clase C n+1 [−1, 1], y sea Pn un polinomio
de grado menor o igual que n que interpola a la función f en los siguientes n + 1
nodos {x0 , x1 , . . . , xn } ⊂ [−1, 1] distintos dos a dos. Entonces, si se eligen los
nodos xk como las raı́ces del polinomio de Chebychev Tn+1 (x), es decir, xk =
cos((2k + 1)π/(2n + 2)), con k = 0, 1, . . . , n, se minimiza la norma del infinito
del error de interpolación máx−1≤x≤1 {|f (x)−Pn (x)|} respecto de cualquier otra
elección de nodos en [−1, 1].
Demostración. Se sabe a partir del Teorema 4.2 que, para todo x ∈ [−1, 1],
existe ξx ∈ (−1, 1) tal que
f (n+1) (ξx )
f (x) − Pn (x) = πn (x) ,
(n + 1)!
n
siendo πn (x) = i=0 (x − xi ). Entonces,
M
máx {|f (x) − Pn (x)|} ≤ máx {|πn (x)|} ,
−1≤x≤1 (n + 1)! −1≤x≤1
A partir de (4.5) y (4.6) se tiene que el valor mı́nimo de máx−1≤x≤1 {|πn (x)|}
se obtiene cuando
n
πn (x) := (x − xi ) = 2−n Tn+1 (x) ,
i=0
es decir, los nodos xk a elegir son las raı́ces del polinomio de Chebychev Tn+1 (x).
Es además obvio que dichas raı́ces son xk = cos((2k + 1)π/(2n + 2)) a partir de
la expresión Tn+1 (x) = cos((n + 1) arc cos x).
Notemos que el Teorema 4.7 está formulado para realizar una interpolación
con nodos únicamente en el intervalo [−1, 1]. Sin embargo, no existe pérdida
de generalidad, puesto que el intervalo x ∈ [−1, 1] se transforma en el intervalo
y ∈ [a, b] mediante el cambio lineal de variables
a+b b−a
y= + x.
2 2
De este modo es fácil mostrar la siguiente generalización del Teorema 4.7
Teorema 4.8 Dada una función f ∈ C n+1 [a, b], el menor error en la norma
del máximo en la interpolación de f se comete cuando se toman las siguientes
abscisas de interpolación {y0 , y1 , . . . , yn } ⊂ [a, b]
En la práctica los splines más utilizados son los cúbicos (k = 3). Se puede ver
que, en la construcción de splines cúbicos existen 2 grados de libertad. Si se
toma S (x0 ) = S (xn ) = 0 los splines se llaman splines cúbicos naturales.
(ii) Sabemos que existe ξ ∈ (0,4, 0,8) tal que el error de interpolación viene
dado por
3
f (iv) (ξ)
ln 0,6 − P3 (0,6) = (0,6 − xi ) .
4! i=0
−6
Como f (iv) (ξ) = ξ4 , obtenemos un error
−0,0004 1
ln 0,6 − P3 (0,6) = .
4 ξ4
1 1 1 1 1 1
0,4 < ξ < 0,8 =⇒ < < =⇒ < 4 < ,
0,8 ξ 0,4 0,84 ξ 0,44
Problema 4.2 Considerar una función f : [a, b] → R de clase C 2 [a, b]. Sean
tres puntos distintos x0 , x1 , x2 ∈ [a, b]. Demostrar que existe un único polinomio
P (x) de grado menor o igual que tres verificando
Solución. Sea P (x) ∈ R3 [x] un polinomio de grado menor o igual que tres,
de manera que P (x) = a0 + a1 x + a2 x2 + a3 x3 , P (x) = a1 + 2a2 x + 3a3 x2 y
4.2 Problemas resueltos 51
Como los tres puntos x0 , x1 , x2 son distintos dos a dos, tenemos que det A = 0
y el sistema lineal es compatible y determinado con lo que demostramos la
existencia y unicidad del polinomio P (x).
Problema 4.3 Consideremos una función f ∈ C n+1 [x0 , xn ]. Sea P (x) el poli-
nomio que interpola a la función f (x) en los nudos x0 < x1 < · · · < xn .
Consideremos un valor x̄ ∈ (x0 , xn ) que no coincida con ninguno de los nudos
anteriores, es decir, x̄ = xi para i = 0, 1, . . . , n.
(i) Utilizando el formalismo de diferencias divididas de Newton, demostrar que
n
f (x̄) − P (x̄) = f [x0 , x1 , . . . , xn , x̄] w(x̄) , siendo w(x) = (x − xi ) .
i=0
de lo que se deduce
f (n+1) (ξ)
f (x̄) − P (x̄) = w(x̄) ,
(n + 1)!
f (n) (ξ)
f [x0 , x1 , . . . , xn ] = .
n!
Problema 4.4 Sea Pn (x) el polinomio de grado menor o igual que n que in-
terpola a la función f (x) = sinh x en cualquier conjunto de n + 1 abscisas en el
intervalo [−1, 1].
(i) Suponiendo que n es impar, hallar una cota del error absoluto |f (x)−Pn (x)|
para x ∈ [−1, 1] que sólo dependa de n.
(ii) Calcular, por el método de las diferencias divididas de Newton, el polinomio
que interpola a la siguiente tabulación.
xi -1 0 1
sinh(xi ) -1.1752 0.0000 1.1752
Como n es impar, n+1 es par y por lo tanto f (n+1) (ξx ) = sinh ξx . Además, como
ξx ∈ (−1, 1) y la función sinh x es impar y creciente se tiene que | sinh ξx | ≤ sinh 1
de manera que
sinh 1
n
|f (x) − Pn (x)| ≤ (x − xi ) .
(n + 1)! i=0
de manera que
h2 a2
Δ≤ <,
2
√
de lo que se deduce h < 2/a. Finalmente, puesto que en la tabulación {0 =
x0 , x1 , . . . , xn = a}
√ se tiene h = a/n, el mı́nimo número de abscisas n viene
dado por n > a2 / 2. Es decir
2
a
n=E √ +1 ,
2
donde E[ ] denota la parte entera.
54 Interpolación Polinómica
1=b+c , 2=b+c ,
3
1=b+c+d , 2=b+c+ d ,
4
Solución. Sea P2 (x) el polinomio de grado menor o igual que 2 que interpola
a todos los puntos de la tabulación del enunciado, es decir, f (xi ) = P2 (xi ) para
i = 0, 1, 2. Como f (x) ∈ C 3 [1, 3], si se realiza la aproximación f (x) ≈ P2 (x)
para cualquier x ∈ [1, 3], el error Δ cometido viene dado por
(3)
f (ξx )
Δ = |f (x) − P2 (x)| = (x − 1)(x − 2)(x − 3) ,
3!
siendo ξx ∈ (x0 , x2 ) = (1, 3). Puesto que f (x) = ln x, se tiene que f (3) (ξx ) =
2/ξx3 . Por otra parte, la aproximación que se quiere efectuar es ln 1,5 = f (1,5) ≈
P2 (1,5), de manera que el error cometido es
2/ξx3 0,375
Δ=
(1,5 − 1)(1,5 − 2)(1,5 − 3) = .
3! 3ξx3
Como 1/ξx3 es una función decreciente en el intervalo ξ ∈ (1, 3), se concluye que
0,375
Δ< = 0,125 .
3
Problema 4.8 Se ha de realizar una interpolación polinomial cuadrática de
la función f (x) = ex sin x para valores de x ∈ [0, 1]. Hallar una cota superior
del error absoluto que se cometerı́a al aproximar la función f por el polinomio
interpolador.
En definitiva
|f (ξx )|
Δ≤ ,
3!
donde la tercera derivada de la función f viene dada por
Es fácil ver que la función f (ξx ) es decreciente en el intervalo (0, 1) puesto que
f iv (ξx ) = −4eξx sin ξx < 0 para todo ξx ∈ (0, 1). De este modo se tiene que los
extremos absolutos de f (ξx ) en el compacto [0, 1] ocurren en los extremos de
dicho intervalo. Puesto que f (0) = 2 y f (1) = −1,63 se tiene la siguiente
acotación del error
|f (0)| 2 1
Δ< = = .
3! 3! 3
Problema 4.9 Se han realizado las siguientes mediciones de la capacidad
calorı́fica molal a presión constante c del nitrógeneo para diferentes temperaturas
T obteniéndose
|f (ξ)| 3 4
Δ≤ 2 = |f (ξ)| ,
3! 3
siendo ξ ∈ (0, 2). Calculando la tercera derivada de la función f (x) = x4 (5 −
x)/60 se obtiene f (ξ) = −(ξ − 1)2 + 1. La gráfica de f (ξ) es una parábola
con un máximo relativo en ξ = 1, de manera que podemos realizar una última
acotación de la forma
4 4
Δ≤ f (1) = = 1,3333 .
3 3
Nota: En la acotación efectuada se ha utilizado que |π(t)| := |t(t − 1)(t −
2)| < 23 para todo t ∈ [0, 2]. Existe, por supuesto, otra forma mejor de acotar
|π(t)| para todo t ∈ [0, 2] que consiste en obtener su máximo absoluto en el
2
intervalo compacto√ [0, 2]. Ası́, ya que las raı́ces del polinomio π (t) = 2 − 6t + 3t
son t± = 1 ± 1/ 3 ∈ [0, 2], es claro que π(t) tiene extremos relativos en t± .
Además, para obtener los extremos absolutos calculamos los valores de π(t) en
los extremos del intervalo cerrado π(0) = π(2) = 0 y en los extremos relativos
|π(t± )| = 0,3849. Se tiene en definitiva que |π(t)| ≤ 0,3849 para todo t ∈ [0, 2]
y por lo tanto una mejor cota del error
Problema 4.11 Sean f (x) y g(x) dos funciones reales definidas sobre el inter-
valo [a, b]. Definamos Pn (x) y Qn (x) como los polinomios de grado menor o igual
que n que interpolan a f y g respectivamente en los nodos {x0 , x1 , . . . , xn } ⊂
[a, b] distintos entre sı́.
(i) ¿Es αPn + βQn el polinomio de grado menor o igual que n que interpola a
la función αf + βg en los nodos {x0 , x1 , . . . , xn }, siendo α, β ∈ R?
(ii) ¿Es Pn Qn el polinomio de grado menor o igual que n que interpola a la
función f g en los nodos {x0 , x1 , . . . , xn }?
Solución. Por definición, si Pn (x) y Qn (x) son los polinomios que interpolan a
f y g respectivamente en los nodos {x0 , x1 , . . . , xn }, entonces, por ser los nodos
4.2 Problemas resueltos 59
distintos entre sı́, Pn (x) y Qn (x) son polinomios de grado menor o igual que
n, es decir, Pn , Qn ∈ Rn [x]. Además Pn (xi ) = f (xi ) y Qn (xi ) = g(xi ) para
i = 0, 1, . . . , n.
es el spline cúbico natural que interpola a los puntos (0, 1), (1, 1), (2, 0) y (3, 10).
es el spline cúbico natural que interpola a los puntos (x0 , f0 ), (x1 , f1 ), (x2 , f2 )
y (x3 , f3 ) si se verifica lo siguiente: (i) si (x) ∈ R3 [x] para i = 0, 1, 2; (ii) S ∈
C 2 [x0 , x2 ]; (iii) S(xi ) = fi para i = 0, 1, 2; (iv) S (x0 ) = S (x2 ) = 0.
Los puntos de interpolación son (x0 , f0 ) = (0, 1), (x1 , f1 ) = (1, 1), (x2 , f2 ) =
(2, 0) y (x3 , f3 ) = (3, 10), de modo que también se verifica la condición
(iii) ya que
Como lı́mx→2− S (x) = lı́mx→2+ S (x) se tiene que S(x) es derivable dos veces
con continuidad en x = 2.
Se concluye que la función S verifica la condición (ii). Finalmente, puesto que
S (x0 ) = S (0) = s0 (0) = 0 y S (x3 ) = S (3) = s2 (3) = 0 también se cumple
la condición (iv). En resumen, S(x) es el spline cúbico natural que interpola a
los puntos (0, 1), (1, 1), (2, 0) y (3, 10).
Problema 4.13 Sea Pn (x) el polinomio de grado menor o igual que n que
interpola a la función f (x) = sinh x en n + 1 puntos del intervalo [−2, 0]. Si
n es impar, hallar una cota que sólo dependa de n del error absoluto Δn (x)
cometido en la aproximación f (x) ≈ Pn (x) para cualquier x ∈ [−2, 0]. ¿Se
puede afirmar a partir de la cota de error anterior que lı́mn→∞ Pn (x) = f (x)
para todo x ∈ [−2, 0]?
4.2 Problemas resueltos 61
de modo que
lı́m Pn (x) = f (x) para todo x ∈ [−2, 0] .
n→∞
Calculemos ahora f (n+1) (x) suponiendo cierta la ecuación (4.8). Puesto que
d (n) d n! (n + 1)!
f (n+1) (x) = f (x) = (−1)n = (−1)n+1
dx dx (1 + x)n+1 (1 + x)n+2
continua obedeciendo la ley dada en (4.8), el método de inducción concluye que
(4.8) es correcta.
Δ= |x − xi | ,
3! i=0
donde ξx ∈ (0, 3) y xi son los nodos. Existe una acotación inmediata del término
con el productorio: puesto que x, xi ∈ [0, 3] para i = 0, 1, 2, entonces |x − xi | ≤ 3
2
y por lo tanto i=0 |x − xi | ≤ 33 . En definitiva se obtiene
|f (ξx )| 3 9
Δ≤ 3 = |f (ξx )| .
3! 2
Puesto que
1
f (x) = 1 + x(6 + 3x + 5x2 ) + sinh(1 − x) ,
6
se obtiene derivando tres veces que f (ξx ) = 5 − cosh(1 − ξx ). Calculemos a
continuación el máximo absoluto de la función continua |f (x)| en el intervalo
compacto [0, 3]. Para ello calculamos primero los extremos relativos de f en el
intervalo [0, 3], es decir, resolvemos la ecuación f (iv) (x) = sinh(1 − x) = 0 para
x ∈ [0, 3]. La solución de dicha ecuación es x = 1. Entonces, el máximo absoluto
M de |f (x)| en el intervalo [0, 3] viene dado por
M = máx{|f (0)|, |f (1)|, |f (3)|} = máx{3,45692, 4, 1,2378} = 4 .
Concluimos que
9
Δ≤ × 4 = 18 .
2
(ii) Los nodos xi ∈ [0, 3] (i = 0, 1, 2) de modo que se minimize f (x) −
Q2 (x)∞ = máx0≤x≤3 {|f (x) − P2 (x)|} son, por definición, los nodos de Cheby-
chev en el intervalo [a, b] = [0, 3], es decir,
a+b b−a 2i + 1
xi = + cos π , i = 0, 1, . . . , n .
2 2 2n + 2
En nuestro caso n = 2, de modo que
3 1
x0 = 1 + cos π = 2,79904 ,
2 6
3 1
x1 = 1 + cos π = 1,5 ,
2 2
3 5
x2 = 1 + cos π = 0,200962 .
2 6
Puesto que f (x0 ) = 23,0517, f (x1 ) = 5,9164 y f (x2 ) = 2,11474, se obtiene el
siguiente esquema triangular de diferencias divididas de Newton
64 Interpolación Polinómica
2.79904 23.0517
13.1907
1.5 5.9164 3.95068
2.92652
0.200962 2.11474
Problema 4.16 Sea Pn (x) el polinomio de grado menor o igual que n que
interpola a la función f (x) en los nodos {x0 , x1 , . . . , xn } ⊂ [a, b] distintos dos a
dos.
(i) Si f (x) = xn , calcular Pn (x).
(ii) Si f (x) = xn+1 , calcular Pn (x) utilizando la fórmula del error en la inter-
polación.
(iii) Si f (x) = xn+2 y [a, b] = [−2, −1,5], acotar en función de n el error
absoluto cometido en la aproximación f (x) ≈ Pn (x) para todo x ∈ [a, b].
¿Qué se puede decir de lı́mn→∞ Pn (x)?
f (n+1) (ξx )
f (x) − Pn (x) = πn (x) ,
(n + 1)!
n
siendo ξx ∈ (a, b) y πn (x) = i=0 (x − xi ). Como f (x) = x
n+1
, se tiene
(n+1)
f (x) = (n + 1)!, de modo que
n
Pn (x) = f (x) − πn (x) = xn+1 − (x − xi ) .
i=0
(iii) Ahora f (x) = xn+2 , de modo que f (n+1) (x) = (n + 2)!x. Entonces,
utilizando de nuevo la fórmula del error en la interpolación, para todo x ∈ [a, b],
el error absoluto Δn (x) cometido en la aproximación f (x) ≈ Pn (x) es
(n+1)
f (ξx )
Δn (x) = |f (x) − Pn (x)| = πn (x) = (n + 2)|ξx ||πn (x)| ,
(n + 1)!
4.2 Problemas resueltos 65
con ξx ∈ (a, b). Como [a, b] = [−2, −1,5] y x, xi ∈ [a, b] es claro que
n
n+1
1
|πn (x)| = |x − xi | < .
i=0
2
Por otra parte, como ξx ∈ (a, b) = (−2, −1,5) se tiene |ξx | < |−2| = 2. Utilizando
las dos acotaciones halladas se concluye que
n+1
1 n+2
Δn (x) < 2(n + 2) = .
2 2n
Finalmente, puesto que 0 ≤ Δn (x) < n+2
2n y lı́mn→∞ n+2
2n = 0, se tiene
lı́mn→∞ Δn (x) = 0. Entonces
lı́m |f (x) − Pn (x)| = lı́m Δn (x) = 0 ,
n→∞ n→∞
de modo que
lı́m Pn (x) = f (x) .
n→∞
√
Problema 4.17 Considerar la función f (x) = (1 − x)/x2 .
(i) Calcular el polinomio P2 (x) de grado menor o igual que 2 que interpola a
f (x) en los nodos {2, 3, 4}.
(ii) Hallar una cota del error absoluto cometido al aproximar f (x) ≈ P2 (x)
para cualquier x ∈ [2, 4].
Solución. (i) Sea√ P2 (x) el polinomio de grado menor o igual que 2 que interpola
a f (x) = (1 − x)/x2 en los nodos {2, 3, 4}. Entonces, la gráfica de P2 (x)
pasa por los puntos (2, f (2)) = (2, −0,103553), (3, f (3)) = (3, −0,081339) y
(4, f (4)) = (4, −0,0625). Es fácil obtener (mediante un esquema de diferencias
divididas o bien con el método de Lagrange) la expresión
P2 (x) = −0,1581 − 0,0306x − 0,0016x2 .
(ii) Sea Δ = |f (x) − P2 (x)|. Como f ∈ C ∞ para todo x > 0, en parti-
cular, f ∈ C 3 ([2, 4]). Entonces, podemos utilizar la expresión del error en la
interpolación
|f (ξx )|
Δ= |π2 (x)| ,
3!
siendo π2 (x) = (x − 2)(x − 3)(x − 4) y ξ ∈ (2, 4). Es evidente que, para todo
x ∈ [2, 4] se tiene que |π2 (x)| < 23 . Entonces, se tiene una primera acotación
|f (ξx )| 3 M 3
Δ< 2 ≤ 2 ,
3! 3!
siendo M = máx2≤ξx ≤4 {|f (ξx )|}. Derivando tres veces la función f se tiene
√
3 35 ξx − 64
f (ξx ) = .
8 ξx5
66 Interpolación Polinómica
En definitiva,
0,169951 3
Δ< 2 = 0,226601 .
3!
Problema 4.18 Se desea construir una tabulación de la función f (x) =
x
0
exp(−t2 )dt en abscisas equiespaciadas de modo que posteriormente se uti-
lizará para aproximar los valores de f (x) para cualquier x ∈ R mediante in-
terpolación lineal entre dos abscisas consecutivas. ¿Cuál debe ser la longitud de
paso en dichas abscisas si se desea que el error absoluto cometido en cualquiera
de las anteriores aproximaciones sea menor que 10−3 ?
Solución. Sea h > 0 la longitud de paso en las abscisas consecutivas que defin-
imos como {x0 , x1 , . . . , xn }, Se tiene pues que xi = x0 + ih para i = 0, 1, . . . , n.
Tomemos un punto arbitrario x entre dos abscisas consecutivas, es decir, x ∈
[i]
[xi , xi+1 ]. Sea P1 (x) el polinomio de grado menor o igual que 1 que interpola
a f (x) en los nodos {xi , xi+1 }. Puesto que la función exp(−x2 ) ∈ C ∞ (R), se
tiene que f (x) ∈ C ∞ (R) y, en particular, f (x) ∈ C 2 ([xi , xi+1 ]). Entonces, si de-
[i]
notamos por Δ(x) el error absoluto cometido en la aproximación f (x) ≈ P1 (x)
para cualquier x ∈ [xi , xi+1 ], se tiene que existe un número ξx ∈ (xi , xi+1 ) tal
que
|f (ξx )|
Δ(x) = |w1 (x)| , siendo w1 (x) = (x − xi )(x − xi+1 ) .
2!
Es claro que |w1 (x)| ≤ (xi+1 − xi )2 = h2 para todo x ∈ [xi , xi+1 ]. Entonces, se
tiene una primera acotación
h2
Δ(x) ≤ |f (ξx )| .
2
Utilizando el Teorema Fundamental del Cálculo se tiene que f (x) = exp(−x2 ).
Entonces, derivando de nuevo, se llega a que
Hemos de hallar (si es que existe) una cota superior de |f (x)| para todo x ∈ R
puesto que, aunque x ∈ [xi , xi+1 ], se tiene que el intervalo [xi , xi+1 ] ⊂ R es
arbitrario. Recordemos que R no es un intervalo compacto (cerrado y acotado),
de modo que en un principio no tiene ni siquiera por qué existir la cota superior
que buscamos. Sin embargo, para la función que tenemos sı́ que existe como
4.2 Problemas resueltos 67
0.5
3 2 1 1 2 3
0.5
Problema 4.19 Consideremos la función f (x) = (4 + x)5 (328 − 19x + x2 ).
(i) Hallar, mediante el método de las diferencias divididas de Newton, el poli-
nomio P2 (x) de grado menor o igual que 2 que interpola a la función f (x)
en los tres nodos equiespaciados de la partición del intervalo [0, 2].
(ii) Acotar el error absoluto Δ(x) cometido en la aproximación f (x) ≈ P2 (x)
para todo x ∈ [0, 2].
(iii) Hallar el polinomio Q2 (x) de grado menor o igual que 2 tal que la cota del
error absoluto Δ̄(x) cometido en la aproximación f (x) ≈ Q2 (x) para todo
x ∈ [0, 2] sea la mı́nima posible.
Δ(x) = |x − i| ,
3! i=0
donde ξx ∈ (0, 2). Puesto que x ∈ [0, 2], entonces para i = 0, 1, 2, se tiene
2
|x − i| ≤ 2 y por lo tanto i=0 |x − i| ≤ 23 . En definitiva se obtiene
|f (ξx )| 3 4
Δ(x) ≤ 2 = |f (ξx )| ,
3! 3
con ξx ∈ (0, 2). Derivando tres veces la función f (x) se llega a
315 x(x − 1)
f (x) = √ .
8 4+x
4.2 Problemas resueltos 69
Calculemos a continuación el máximo absoluto de la función continua |f (x)| en
el intervalo compacto [0, 2]. Para ello calculamos primero los extremos relativos
de f en el intervalo [0, 2], es decir, resolvemos la ecuación
315 −8 + 15x + 3x2
f (iv) (x) = =0
16 (4 + x)3
para x ∈ [0, 2]. Las soluciones de dicha ecuación son las soluciones de la ecuación
−8 + 15x + 3x2 = 0 que vienen dadas por
1 √
x = (−15 ± 321) ,
6
es decir, x = 0,486079 ∈ [0, 2] puesto que la otra solución −5,48608 ∈ [0, 2].
Entonces, el máximo absoluto M de |f (x)| en el intervalo [0, 2] viene dado por
M = máx{|f (0)|, |f (2)|, |f (0,486079)|} = máx{0, 32,1496, 4,64398}
= 32,1496 .
Concluimos que
4
Δ(x) ≤ × 32,1496 = 42,8661 .
3
(iii) El polinomio Q2 (x) deberá interpolar a la función f (x) en los 3 nodos
xi (i = 0, 1, 2) de Chebychev en el intervalo [a, b] = [0, 2]. Entonces,
a+b b−a 2i + 1
xi = + cos π , i = 0, 1, 2 .
2 2 2×2+2
En nuestro caso
1
x0 = 1 + cos π = 1,86603 ,
6
1
x1 = 1 + cos π =1,
2
5
x2 = 1 + cos π = 0,133975 ,
6
Derivación Numérica
71
72 Derivación Numérica
xi 5 -7 -6 0
f (xi ) 1 -23 -54 -954
(i) Calcular el polinomio de tercer grado que interpola a todos los puntos de la
tabla mediante el método de las diferencias divididas de Newton.
(ii) Aproximar el valor de la derivada f (−6,5) mediante alguna de las siguien-
tes fórmulas: diferencias hacia adelante, diferencias hacia atrás o diferen-
cias centradas. Explicar vuestra elección.
5.2 Problemas resueltos 73
Problema 5.3 Se está analizando un circuito eléctrico que tiene una bobina,
una resistencia y un generador de corriente. Aplicando las leyes de Kirchhoff se
sabe que dicho circuito está gobernado por la siguiente ecuación
dI
(t) = L + RI ,
dt
siendo (t) la fuerza electromotriz en función del tiempo del generador, L el
coeficiente de autoinducción de la bobina, R la resistencia e I la intensidad de
corriente eléctrica. Se realiza con un cronómetro y un amperı́metro las siguientes
mediciones
ti (segundos) 0 1 2
I(ti ) (amperios) 2 4 10
Δ(0) ≈ LΔI(0)
˙ + |I(0)|ΔR = 0,5ΔI(0)
˙ + 2 × 0,3 . (5.7)
donde ξ(t) ∈ (0, 2). Derivando esta ecuación respecto del tiempo t y particula-
rizándola para t = 0 se obtiene
(iii)
I (ξ(0)) 1 (iii)
ΔI(0) := | ˙
I(0) − Ṗ (0)| = (−1)(−2) = I (ξ(0)) ,
˙ 2 3! 3
f (x + h) − f (x − h)
f (x) = + O(h2 ) .
2h
(i) Aplicar la extrapolación de Richardson a la fórmula de diferencias centradas
para obtener una nueva fórmula de aproximación para f (x), obteniendo
también el orden del error.
(ii) Sea f ∈ C 4 [x − 2h, x + 2h], con h > 0. Utilizando desarrollos de Taylor
hasta cuarto orden para f (x ± h) y f (x ± 2h), obtener la misma fórmula
que en el apartado anterior.
f (x + h) − f (x − h)
F (h) = .
2h
Realizando la extrapolación de Richardson con duplicación del paso, se obtiene
el primer término del error, es decir,
F (h) − F (2h)
a1 h2 = + O(h4 )
22 − 1
−f (x + 2h) + 2f (x + h) − 2f (x − h) + f (x − 2h)
= + O(h4 ) .
12h
Se tiene pues que
−f (x + 2h) + 8f (x + h) − 8f (x − h) + f (x − 2h)
f (x) = + O(h4 ) .
12h
(ii) Una forma de obtener expresiones en diferencias finitas para la derivada
consiste (suponiendo que la función es suficientemente derivable) en combinar
diferentes desarrollos de Taylor en un entorno de los puntos donde conocemos
76 Derivación Numérica
−f (x + 2h) + 8f (x + h) − 8f (x − h) + f (x − 2h)
f (x) = + O(h4 ) .
12h
Problema 5.5 Considerar la siguiente tabla de valores
x 0 1 2
f (x) 3.121 2.714 4.128
(i) Sabiendo que |f (x)| ≤ 3/4 para todo x ∈ [0, 2], hallar una cota para el
error cometido cuando aproximamos f (1.3) por interpolación polinomial
con la tabla adjunta.
(ii) Aproximar f (1) por derivación interpolatoria y, teniendo en cuenta el
apartado anterior, obtener una cota del error cometido.
(iii) ¿Se obtendrı́a la misma aproximación que la efectuada en el apartado an-
terior mediante alguna fórmula de derivación numérica obtenida a través
de desarrollos de Taylor?
donde ξx ∈ (0, 2) y π(x) = x(x − 1)(x − 2). Puesto que |f (ξx )| ≤ 3/4, particu-
larizando para x = 1.3 se tiene la siguiente acotación del error Δ cometido en
5.2 Problemas resueltos 77
de manera que
Como ξ1 ∈ (0, 2) se tiene la cota |f (ξ1 )| ≤ 3/4. Además, como π(x) = x(x −
1)(x − 2), obtenemos π (1) = −1. En definitiva la cota del error es
3/4 1
Δ ≤ − = .
3! 8
(iii) Puesto que los nodos de la tabulación son equiespaciados, seguro que
existe alguna fórmula de derivación numérica procedente de desarrollos de Taylor
78 Derivación Numérica
t (s) 0 1 2
x (m) 0 1.1752 3.6268
El polinomio P2 (t) que interpola a la función x(t) según los datos de la tabla
es
P2 (t) = 1,1752t + 0,6382t(t − 1) .
La aproximación mediante derivación interpolatoria es
8 3
f (x + 2h) = f (x) + 2hf (x) + 2h2 f (x) + h f (ξ2 ) ,
3!
2
f (ξx )
Δ(x) = w2 (x) , w2 (x) = (x − xi ) , ξx ∈ (x0 , x2 ) .
3! i=0
con lo que
−3f (x − h) − 10f (x) + 18f (x + h) − 6f (x + 2h) + f (x + 3h)
f (x) = +O(h4 ) .
12h
5.2 Problemas resueltos 83
F (h) − F (2h)
a1 h4 = + O(h5 ) .
24 − 1
Entonces
F (h) − F (2h) 1
f (x) = F (h) + 4
+ O(h5 ) = (16F (h) − F (2h)) + O(h5 ).
2 −1 15
Desarrollando esta expresión se llega a
1
f (x) = [3f (x − 2h) − 96f (x − h) − 310f (x) + 576f (x + h)
360h
−210f (x + 2h) + 32f (x + 3h) + 6f (x + 4h) − f (x + 6h)]
+O(h5 ) .
|f (n+1) (ξ)|
n
Δ(n) = |f (xj ) − Pn (xj )| = |xj − xi | ,
(n + 1)! i=0
i=j
siendo ξ ∈ (0, π/2). Puesto que todos los nodos pertenecen al intervalo [0, π/2],
es claro que |xj − xi | ≤ π/2, de modo que
n π n
|xj − xi | ≤ .
i=0
2
i=j
84 Derivación Numérica
Por otra parte, como f (x) = sin x, se verifica que |f (n+1) (ξ)| ≤ 1. Entonces,
obtenemos la acotación
π n
2
Δ(n) ≤ .
(n + 1)!
(ii) Debido a la cota obtenida en el apartado anterior se tiene que
π n
2
0 ≤ Δ(n) ≤ .
(n + 1)!
Como
π n
2
lı́m =0,
n→∞ (n + 1)!
se tiene que lı́mn→∞ Δ(n) = 0. Entonces, en el problema de aproximar f (x) por
Pn (x) en el intervalo [0, π/2] no se produce el fenómeno Runge.
Capı́tulo 6
Integración Numérica
siendo ! n
αk := φk (t) dt , (6.3)
0
85
86 Integración Numérica
con
n
t−j
φk (t) = , k = 0, . . . , n . (6.4)
j=0
k −j
j=k
máximo grado de precisión. Para verlo basta con tomar el polinomio Π2 (x) =
n 2
i=0 (x − xk ) de grado 2n + 2 y comprobar que
! b
n
0< w(x)Π2 (x) dx = αk Π2 (xk ) = 0 .
a k=0
2k+1
siendo xk = cos 2n+2 π .
Teorema 6.3 Sea f ∈ C 2n+2 [a, b]. Entonces existe un ξ ∈ (a, b) tal que el error
Δ cometido en la fórmula de integración gaussiana con n + 1 puntos viene dado
por
!
n !
b
f (2n+2) (ξ) b
2
Δ := w(x)f (x) dx − αk f (xk ) = w(x)ψn+1 (x) dx ,
a (2n + 2)! a
k=0
donde los nodos xk con k = 0, 1, . . . , n son las raı́ces del polinomio ortogonal
mónico ψn+1 (x) de grado n + 1 en el intervalo [a, b] respecto al peso w(x).
! 1
sin x f0 f5
dx ≈ Th (f ) := h + f1 + f2 + f3 + f4 +
0 x 2 2
f (0) f (1)
= 0,2 + f (0,2) + f (0,4) + f (0,6) + f (0,8) +
2 2
1 0,841
= 0,2 + 0,993 + 0,973 + 0,941 + 0,896 +
2 2
= 0,944 .
sin x
lı́m f (x) = lı́m =0.
x→0 x→0 x
h f0 f10
Th/2 (f ) := + f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 +
2 2 2
f (0)
= 0,1 + f (0,1) + f (0,2) + f (0,3) + f (0,4) + f (0,5)+
2
f (1)
f (0,6) + f (0,7) + f (0,8) + f (0,9) +
2
1
= 0,1 + 0,998 + 0,993 + 0,985 + 0,973 + 0,958+
2
0,841
0,941 + 0,920 + 0,896 + 0,870 + = 0,945
2
! 1
sin x (2) 1 " (1) (1)
#
dx ≈ Th (f ) := 4Th/2 (f ) − Th (f )
0 x 3
1
= [4 × 0,945 − 0,944] = 0,945
3
92 Integración Numérica
6 √
Problema 6.2 (i) Aproximar el valor de la integral 1 2 + sin(2 x) dx uti-
lizando la regla de Simpson compuesta con 10 subintervalos.
b
(ii) Se sabe que el error absoluto Δ cometido al aproximar la integral a
f (x) dx
por la regla de Simpson compuesta Sh (f ) viene dado por
!
b (b − a)f (iv) (ξ)h4
Δ := f (x) dx − Sh (f ) = ,
a 180
donde ξ ∈ (a, b). Calcular el número de puntos n que se han de tomar del
intervalo [2, 7] para obtener mediante la regla de Simpson compuesta una
7
aproximación de la integral 2 1/x dx con un error menor que 5 × 10−9 .
(ii) Sea la función f (x) = 1/x. La derivada cuarta de esta función viene
dada por f (iv) (x) = 24/x5 . Como f (iv) (x) es una función positiva y decreciente
en el intervalo [2, 7] podemos concluir que para cualquier ξ ∈ [2, 7] se puede
realizar la acotación
(iv) (iv) 3
f (ξ) ≤ f (2) = .
4
De esta forma se consigue la siguiente acotación del error
(b − a)f (iv) (ξ)h4 (7 − 2)f (iv) (ξ)h4 (7 − 2) 34 h4 h4
Δ = = ≤ =
180 180 180 48 .
6.2 Problemas resueltos 93
(5/n)4 625
Δ≤ = .
48 48n4
625
Δ≤ < 5 × 10−9 ,
48n4
de donde se puede despejar n para obtener
4 125
n> × 109 = 225,9 .
48
Como n ∈ N y además ha de ser par tomaremos n = 226.
h
L ≈ (f0 + 4f1 + 2f2 + 4f3 + f4 )
3
20/4
= (f (0) + 4f (5) + 2f (10) + 4f (15) + f (20))
3
5 √
= 2 + 4 × 1,03945 + 2 × 1,30539 + 4 × 1,25584 + 1,08006
3
= 23,8104
94 Integración Numérica
1 √
Problema 6.4 Considerar la integral I = 0
f (x) dx siendo f (x) = cos x/ x.
(i) ¿Qué problema existe si se pretende aproximar la integral I por alguna regla
de Newton-Cotes?
(ii) Utilizando que f (x) = f (x) − √1 + √1 , descomponer la integral I en dos
x x
1
integrales de la forma I = I1 + I2 , siendo I1 = 0 √1x dx. Calcular la
integral I mediante el cálculo exacto de I1 y la aproximación numérica de
I2 por la fórmula compuesta de Simpson con 5 abscisas.
I ≈ 2 − 0,19116 = 1,80884 .
donde xk son las abscisas de Chebychev, es decir los ceros del polinomio de
Chebychev Tn+1 (x) = cos[(n + 1) arc cos x]. En concreto
2k + 1
xk = cos π , k = 0, 1, . . . , n .
2n + 2
donde se conocen la función f (x) y el llamado núcleo K(x, t), ası́ como las
constantes a y b. Se pretende aproximar la función y(x) dando su tabulación
y(xk ) con k = 0, 1, . . . , n, donde a = x0 < x1 < · · · < xn−1 < xn = b es una
partición equiespaciada del intervalo [a, b]. Para conseguirlo, se aproximarán las
integrales que aparecen mediante la regla de Simpson.
(i) Si a = 0, b = 2 y el núcleo K(x, t) = x2 − t, hallar el sistema lineal de
ecuaciones que deben verificar las incógnitas y(0), y(1) e y(2).
(ii) La tabulación y(0), y(1) e y(2) se pretende realizar varias veces, donde
en cada tabulación lo único que varı́a es la función f (x). Razonar cual
es el método más eficaz para resolver los distintos sistemas lineales que
aparecen.
(iii) Resolver por el método descrito en el apartado (ii) el sistema lineal del
apartado (i) sabiendo que estamos realizando la tabulación en el caso en
que f (x) = x.
a+b
x0 = a < x1 = < x2 = b .
2
Utilizando la regla de Simpson para calcular las integrales que aparecen
cuando se realiza la tabulación se obtiene
! b
b−a
y(a) = f (a) + K(a, t)y(t) dt ≈ f (a) + [K(a, a)y(a)
a 6
+4K(a, (a + b)/2)y((a + b)/2) + K(a, b)y(b)] ,
! b
y((a + b)/2) = f ((a + b)/2) + K((a + b)/2, t)y(t) dt ≈ f ((a + b)/2)
a
b−a
+ [K((a + b)/2, a)y(a) + 4K((a + b)/2, (a + b)/2)
6
y((a + b)/2) + K((a + b)/2, b)y(b)] ,
6.2 Problemas resueltos 97
! b
b−a
y(b) = f (b) + K(b, t)y(t) dt ≈ f (b) +
[K(b, a)y(a)
a 6
+4K(b, (a + b)/2)y((a + b)/2) + K(b, b)y(b)] ,
⎛ ⎞ ⎛ ⎞
− 43 − 23
−1 −1 − 43 − 23
1 ⎠ ⎝
A = ⎝ 1
3 −1 − 3 ∼ 0 − 13
9 − 59 ⎠
4 1 20
3 4 −3 0 9 − 11
9
⎛ 4 2
⎞
−1 − 3 −3
∼ ⎝ 0 − 13 9 − 59 ⎠ = U ,
0 0 − 27
13
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
−1 − 43 − 23 1 0 0 −1 − 43 − 23
A=⎝ 1
3 −1 − 13 ⎠ = LU = ⎝ − 13 1 0 ⎠ ⎝ 0 − 13
9 − 59 ⎠
4
3 4 0 − 43 − 20
13 1 0 0 − 27
13
⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
α 0 y(0) α
L ⎝ β ⎠ = ⎝ −1 ⎠ , U ⎝ y(1) ⎠ = ⎝ β ⎠ .
γ −2 y(2) γ
(i) Una propiedad notable de dicha función es que Γ(n) = (n−1)! para todo n ∈
N. Utilizar este hecho para aproximar Γ(5/2) mediante una interpolación
polinomial cúbica por el método de las diferencias divididas de Newton.
(ii) Aproximar el valor de Γ(5/2) utilizando el método de integración Gauss-
Laguerre con 2 puntos.
Indicación 1: Los polinomios de Laguerre son ortogonales en el intervalo
[0, ∞) respecto de la función peso w(x) = e−x . Las raı́ces del polinomio
de Laguerre de segundo grado L2 (x) son x0 = 0,58578, x1 = 3,41421.
∞ ∞
Indicación 2: 0 e−x dx = 0 xe−x dx = 1.
Solución. (i) La propiedad Γ(n) = (n−1)! para todo n ∈ N nos permite realizar
una tabulación de la función Γ(x) en los números naturales. Como queremos
aproximar Γ(5/2) mediante una interpolación polinomial cúbica, realizamos la
siguiente tabulación
xi 1 2 3 4
Γ(xi ) 1 1 2 6
! b
n
w(x)f (x) dx ≈ αk f (xk ) ,
a k=0
! ∞
e−x f (x) dx ≈ α0 f (0,58578) + α1 f (3,41421) . (6.11)
0
! ∞
1 = e−x dx = α0 + α1 ,
!0 ∞
1 = xe−x dx = 0,58578α0 + 3,41421α1 .
0
! ∞
e−x f (x) dx ≈ 0,853551f (0,58578) + 0,146449f (3,41421) .
0
! ∞
Γ(5/2) = x3/2 e−x dx ≈ 0,853551 × 0,585783/2 +
0
0,146449 × 3,414213/2 = 1,30657 .
6.2 Problemas resueltos 101
xk 5 7 9
J2 (xk ) 0.0465 -0.3014 0.1448
(i) Aproximar el valor de J2 (6) por interpolación polinomial con todos los datos
de la tabla adjunta y utilizando el formalismo de las diferencias divididas
de Newton y acotar el error cometido.
(ii) Si se pudiera tabular la función J2 (x) en tres nodos cualesquiera del in-
tervalo [5, 9], explicar detalladamente en qué nodos se tabuları́a sabiendo
que posteriormente dicha tabulación se va a usar para realizar una inter-
polación polinomial.
(iii) Aproximar el valor de J2 (6) mediante una extrapolación de Richardson
según el método de Romberg.
(iv) Aproximar el valor de J2 (7) mediante una derivación interpolatoria a par-
tir de la tabulación dada y acotar el arror absoluto asociado. Sin realizar
ningún cálculo ¿se conoce alguna fórmula de aproximación de la deriva-
da J2 (x) por un truncamiento de series de Taylor de la función J2 (x)
alrededor de x que proporcione la misma aproximación de J2 (7) que hemos
obtenido en el cálculo anterior?
(v) Aproximar, a partir de la expresión integral de J2 (x), el valor de J2 (7)
mediante el método de Simpson.
5 0.0465
-0.1739
7 -0.3014 0.0992
0.2231
9 0.1448
con ξ ∈ (5, 9). A continuación obtendremos una cota de |J2 (ξ)|. En primer
lugar, a partir de la expresión integral de J2 (x), calculamos J2 (x). Derivando
sucesivamente respecto de x se tiene
!
−1 π
J2 (x) = sin t sin[x sin t − 2t] dt ,
π 0
! π
−1
J2 (x) = sin2 t cos[x sin t − 2t] dt ,
π 0
!
1 π 3
J2 (x) = sin t sin[x sin t − 2t] dt .
π 0
b b
Utilizando que a f (t) dt ≤ a |f (t)| dt para cualquier función f integrable
Riemann en el intervalo [a, b] y que | sin θ| ≤ 1 para todo θ ∈ [0, 2π) se obtiene
la siguiente cadena de acotaciones
!
1 π
|J2 (x)| ≤ | sin3 t| | sin[x sin t − 2t]| dt
π 0
!
1 π
≤ | sin3 t| dt .
π 0
Puesto que sin t ≥ 0 para 0 ≤ t ≤ π, se puede quitar el valor absoluto de la
última integral, de manera que
!
1 π 3
|J2 (x)| ≤ sin t dt .
π 0
Esta integral se calcula de manera exacta utilizando la identidad
Se tiene pues que, a partir de la ecuación (6.12), el error absoluto Δ está acotado
de la forma
|J (ξ)| 2
Δ= 2 ≤ = 0,212207 .
2 3π
(ii) Se sabe que, si la interpolación se efectúa en el intervalo −1 ≤ x ≤ 1,
entonces las abscisas de Chebychev definidas de la forma
2k + 1
xk = cos ∈ (−1, 1) , k = 0, 1, . . . , n ,
2n + 2
π π
Tπ(1) (f ) =[f (0) + f (π)] = [1 + 1] = π = 3,14159 ,
2 2
π f (0) π f (π) π
1 1
(1)
Tπ/2 (f ) = +f + = − 0,96017 +
2 2 2 2 2 2 2
= 0,06256 .
104 Integración Numérica
(iv) Una aproximación de J2 (7) mediante derivación interpolatoria con todos
los datos de la tabla adjunta consiste en realizar la aproximación J2 (7) ≈ P2 (7)
siendo P2 (x) el polinomio de interpolación calculado en el apartado (i). Puesto
que P2 (x) = 0,0465 − 0,1739(x − 5) + 0,0992(x − 5)(x − 7), se tiene
Si definimos Δ = |J2 (7) − P2 (7)| como el error absoluto cometido en la anterior
aproximación de la derivada, se sabe que
J2 (ξ) 2
Δ = (7 − 5)(7 − 9) = |J2 (ξ)| , ξ ∈ (5, 9) .
3! 3
En el apartado (i) hemos obtenido una cota de |J2 (x)| para cualquier valor de
4
x ∈ R. En concreto se tiene |J2 (x)| ≤ 3π . Por lo tanto tenemos la siguiente cota
del error Δ cometido en la derivación interpolatoria efectuada anteriormente
2 4 8
Δ ≤ = = 0,282942 .
3 3π 9π
Por último, recordemos que las fórmulas de derivación interpolatoria con
nodos equiespaciados xk = x0 +kh, siendo h la longitud de paso, son totalmente
equivalentes a las fórmulas que aproximan la derivada obtenidas a partir de un
truncamiento de series de Taylor. Puesto que los nodos {x0 = 5, x1 = 7, x2 =
9} donde se ha realizado la interpolación del apartado (i) son equiespaciados
con paso h = 2 y además, en este apartado se ha aproximado la derivada de
la función J2 (x) en el nodo central x1 se tiene que la fórmula de derivación
aproximada equivalente obtenida por un método de Taylor debe ser la fórmula
de diferencias centradas
J2 (x + h) − J2 (x − h)
J2 (x) ≈ .
2h
La comprobación de este resultado es la siguiente. Aplicando la fórmula de
diferencias centradas se obtiene
J2 (9) − J2 (5)
J2 (7) ≈ = 0,0245 ,
4
6.2 Problemas resueltos 105
que coincide con la aproximación J2 (7) ≈ P2 (7) realizada anteriormente.
(v) En este apartado se pretende aproximar el valor de J2 (7) por el método
de Simpson. En consecuencia se debe de expresar J2 (7) mediante una integral
de Riemann. Esta expresión integral la hemos obtenido en el apartado (i). En
concreto !
−1 π
J2 (7) = g(t) dt , g(t) = sin t sin[7 sin t − 2t] .
π 0
Realizando la aproximación de la anterior integral por el método de Simpson se
obtiene
−1 −1 π/2 " π #
J2 (7) ≈ Sπ/2 (g) := g(0) + 4g + g(π)
π π 3 2
−1
= [0 − 4 × 0,656987 + 0] = 0,437991 .
6
Existe una forma alternativa de resolver el problema como se explica a con-
tinuación. En primer lugar se aproxima, por el método de Simpson, el valor de
J2 (x). Puesto que
!
1 π
J2 (x) = f (t) dt , f (t) = cos[x sin t − 2t] ,
π 0
se tiene
1 1 π/2 " π #
J2 (x) ≈ Sπ/2 (f ) := f (0) + 4f + f (π)
π π 3 2
1 1
= [1 + 4 cos(x − π) + 1] = [1 + 2 cos(x − π)] .
6 3
Finalmente, derivamos la expresión anterior
d 1 −2
J2 (x) ≈ [1 + 2 cos(x − π)] = sin(x − π) ,
dx 3 3
y la evaluamos en x = 7 para obtener
−2
J2 (7) ≈ sin(7 − π) = 0,437991 .
3
Problema 6.9 Aproximar el valor de ln 2 a partir de la expresión integral
! 2
1
ln 2 = dx ,
1 x
(ii) Utilizar el resultado del apartado anterior para aproximar el calor total Q
absorbido por 23 moles de nitrógeno en un proceso isobárico cuasi-estático
sabiendo que su temperatura inicial es 350 K y la final 800 K.
Indicación: El calor Q absorbido por n moles en un proceso isobárico
cuasi-estático con temperatura inicial Ti y final Tf viene dado por la ex-
T
presión Q = n Tif cp dT .
donde en el último paso se han tenido en cuenta las propiedades del valor abso-
luto y de la integral de Riemann. Puesto que P2 (T ) es el polinomio que interpola
a cp en los nodos T0 , T1 y T2 , se tiene que existe un valor ξT ∈ (Ti , Tf ) tal que
cp (T ) − P2 (T ) = c
p (T )π(T )/3!, siendo π(T ) = (T − T0 )(T − T1 )(T − T2 ). En
definitiva obtenemos que
!
n Tf
Δ≤ |cp (T )||(T − T0 )(T − T1 )(T − T2 )| dT ,
3! Ti
de manera que, teniendo en cuenta que |c
p (T )| ≤ 10
−13
para toda T ∈ [350, 800]
y substituyendo los datos del enunciado obtenemos
!
23 800
Δ ≤ 10−13 |(T − 300)(T − 1500)(T − 2100)| dT
3! 350
! 800
−13 23
= 10 (T − 300)(T − 1500)(T − 2100) dT
6 350
! 800
23
= 10−13 (−945000000 + 4230000T − 3900T 2 + T 3 ) dT
6 350
23
= 10−13 158048437500 = 0,0605852 .
6
108 Integración Numérica
siendo 2
1 x
f (x) = √ exp − .
2π 2
Observemos que la función f (x) no admite primitiva expresable mediante fun-
ciones elementale, de modo que, su aproximación numérica está perfectamente
justificada. Utilizando ahora la fórmula de Simpson compuesta Sh (f ) con n = 4
subintervalos se tiene que la longitud de paso vale h = (b − a)/n = 5/2, de modo
que
5/2
P (10 ≤ Z ≤ 20) ≈ S5/2 (f ) = (f0 + 4f1 + 2f2 + 4f3 + f4 )
3
5
= [f (10) + 4f1 (25/2) + 2f (15) + 4f (35/2) + f (20)]
6
5 1 " −102 /2 2 2
= √ e + 4e−(25/2) /2 + 2e−15 /2
6 2π
2 2
#
+4e−(35/2) /2 + e−20 /2
= 0,332452[1,92875 × 10−22 ] = 6,41217 × 10−23 .
Problema 6.12 Utilizar el método de los coeficientes indeterminados para ha-
llar una fórmula de integración numérica del tipo
! 3
f (x) dx ≈ af (0) + bf (2) + cf (4) ,
1
con a, b, c ∈ R, que sea exacta para polinomios del máximo grado posible. Dar
dicho grado máximo.
6.2 Problemas resueltos 109
de modo que
lı́m f (z) = 1 . (6.14)
z→0+
de donde
16 & '
a1 h4 = Sh/2 (f ) − Sh (f ) + O(h5 ) .
15
Sustituyendo este término en (6.15) se obtiene la siguiente fórmula mejorada
! b
f (x) dx = F̃ (h) + O(h5 ) ,
a
siendo
16 & ' 1 & '
F̃ (h) = Sh (f ) + Sh/2 (f ) − Sh (f ) = 16Sh/2 (f ) − Sh (f ) .
15 15
2κλ
V (r, θ) = √ K(κ) ,
Rr sin θ
siendo λ = Q/(2πR) la densidad lineal de carga en el anillo,
4Rr sin θ
κ2 = ,
R2 + r2 + 2Rr sin θ
y K(κ) la integral elı́ptica de primera especie
! π/2
dφ
K(κ) = .
0 1 − κ2 sin2 φ
(i) Utilizando el método de los trapecios compuesto con una división de 3 subin-
tervalos, aproximar el potencial eléctrico en cualquier punto de la circun-
ferencia (r, ϕ, θ) = (1, ϕ, π/2) con 0 ≤ ϕ < 2π concéntrica al anillo.
Datos: Q = 10C, R = 5m.
(ii) ¿Es necesario utilizar métodos numéricos para calcular el potencial eléctrico
en cualquier punto del eje z?
Finalmente
2
V (1, π/2) = K( 5/9) ≈ 0,4041 .
3π
(ii) Nótese que los puntos del eje z verifican θ = 0. De este modo, en dichos
puntos se tiene κ = 0 con lo cual el potencial eléctrico en cualquier punto del
eje z será
4λ
V (r, 0) = √ K(0) ,
r + R2
2
x
0.5 1 1.5 2
-1
-2
-3
-4
-5
h
Cv ≈ 67,8472Sh (f ) = 67,8472 × (f0 + 4f1 + 2f2 + 4f3 + f4 )
3
0,2583
= 67,8472 × (f (0) + 4f (0,2583) + 2f (0,5166) + 4f (0,7749)
3
+f (1,0333))
Entonces,
0,2583
Cv ≈ 67,8472 × (0 + 4 × 0,0663 + 2 × 0,2610 + 4 × 0,5713
3
+0,9775) = 67,8472 × 0,348709 = 23,659 .
(ii) Se sabe que, si f ∈ C 4 [a, b], el error absoluto Δ cometido en la aproxi-
mación ! b
f (x) dx ≈ Sh (f ) ,
a
siendo Sh (f ) la fórmula de Simpson compuesta con una longitud de paso h
aplicada a la función f viene dado por
4
h
Δ = (b − a)f (iv) (ξ) , ξ ∈ (a, b) .
180
En nuestro caso, el error absoluto cometido en la aproximación
! 1,0333
f (x) dx ≈ S0,2583 (f ) ,
0
Teniendo en cuenta la gráfica de la función f (iv) , es claro que |f (iv) (ξ)| < 2 para
todo ξ ∈ (0, 1,0333). Entonces
Δ < 0,000025 × 2 = 0,000051 .
Finalmente, como la aproximación efectuada en el apartado (i) ha sido
! 1,0333
Cv = 67,8472 f (x) dx ≈ 67,8472S0,2583 (f ) ,
0
es impropia de segunda especie puesto que la función integrando f (x) no está aco-
tada en el intervalo de integración [0, 1]. En particular, la gráfica de f (x) tiene
una ası́ntota vertical en x = 1. Realizando el cambio de variable x = sin t en la
integral I se obtiene
! π/2 ! π/2
sin(sin t)
I= cos t dt = sin(sin t) dt .
0 1 − sin2 t 0
Notemos que, tras el cambio de varible, la integral I deja de ser impropia puesto
que la función integrando g(t) = sin(sin t) está acotada en el intervalo de inte-
gración [0, π/2].
A continuación aproximamos el valor de I mediante la regla de los trapecios
compuesta con n = 4 subintervalos, es decir, con una longitud de paso h =
(π/2 − 0)/4 = π/8.
g g4
0
I ≈ Th (g) = h + g1 + g2 + g3 +
2 2
π g(0) g(π/2)
= + g(π/8) + g(π/4) + g(3π/8) +
8 2 2
π 0 0,841471
= + 0,373411 + 0,649637 + 0,797946 + = 0,880325 .
8 2 2
b
(ii) El error absoluto Δ en la aproximación a
g(t) dt ≈ Th (g) viene dado
por 2
h
Δ = (b − a)g (ξ) , ξ ∈ (a, b) .
12
En nuestro caso se tiene
(π/8)2 π 3
Δ= (π/2 − 0)|g (ξ)| = |g (ξ)| ,
12 1536
siendo ξ ∈ (0, π/2). Derivando 2 veces la función g(t) = sin(sin t) se tiene
g (t) = − cos(sin x) sin x − cos2 x sin(sin x). Para hallar el valor máximo (es
decir, la menor de las cotas superiores) de Δ se ha de obtener el valor
1
= (f (0) + 4f (1/4) + 2f (1/2) + 4f (3/4) + f (1))
12
1
= (0 + 4 × 0,0624593 + 2 × 0,247404 + 4 × 0,533303 + 0,841471)
12
= 0,309944 .
1/44 (iv)
Δ= |f (ξ)| ,
180
siendo f (iv) (ξ) = 4[(4ξ 4 − 3) sin(ξ 2 ) − 12ξ 2 cos(ξ 2 )]. Aplicando las propiedades
del valor absoluto y teniendo en cuenta que las funciones seno y coseno están
acotadas entre 1 y -1 se tiene la siguiente acotación
|f (iv) (ξ)| ≤ 4[|4ξ 4 − 3|| sin(ξ 2 )| + 12ξ 2 | cos(ξ 2 )|] < 4[|4ξ 4 − 3| + 12ξ 2 ]
< 4[3 + 12] = 60,
1/44
Δ< 60 = 0,00130208 .
180
Problema 6.20 La teorı́a de estática de fluidos dice que la presión P que ejerce
un fluido de densidad ρ a una profundidad x de su superficie viene dada por
P (x) = ρgx, siendo g la intensidad gravitatoria. Entonces, la fuerza total F
que ejerce un fluido estancado una profundidad H sobre la compuerta de una
H
presa es F = ρg 0 xL(x) dx donde L(x) es la anchura de la compuerta a una
profundidad x.
(i) Sabiendo que el fluido es agua (ρ = 103 ), H = 8 y que se han realizado las
siguientes mediciones de la función L(x)
x 0 2 4 6 8
L(x) 150 20 2.5 0.5 0
x 0 2 4 6 8
f (x) 0 40 10 3 0
Notemos que las abscisas de la tabulación son equiespaciadas con una longitud
de paso h = 2. Entonces, podemos aproximar mediante la fórmula de Simpson
compuesta la integral
! 8
h
f (x) dx ≈ (f0 + 4f1 + 2f2 + 4f3 + f4 )
0 3
2
= (0 + 4 × 40 + 2 × 10 + 4 × 3 + 0) = 128 .
3
En definitiva, tomando la gravedad g = 9,8, se tiene
0 150
-65
2 20 14.062
-8.75 -2.027
4 2.5 1.9 0.2177
-1 -0.285
6 0.5 0.1875
-0.25
8 0
e−1/t k 5/2 ∞
lı́m f (t) = lı́m = lı́m k = =0,
t→0+ t→0+ t5/2 k→∞ e ∞
donde la indeterminación ∞/∞ se ha resuelto teniendo que el numerador es un
infinito de orden inferior al del denominador. También se puede ver aplicando
tres veces la regla de l’Hopital. En definitiva, la integral I se ha transformado
en una integral propia y por lo tanto convergente. Ahora podemos utilizar para
su aproximación la fórmula de los trapecios compuesta con 4 subintervalos:
! 1
f0 f4
f (t) dt ≈ h + f1 + f2 + f3 + ,
0 2 2
obtiene que
! 1 " #1
1 π
I= 2
dz = arctan z = arctan 1 = .
0 1+z 0 4
Aproximar el área que genera la curva y = tan x al rotar alrededor del eje
de abscisas desde x = 0 hasta x = 1 mediante una única extrapolación de
Richardson según el método de Romberg.
Notar que se ha quitado el valor absoluto puesto que tan x ≥ 0 para todo
x ∈ [0, 1]. Vamos a aproximar la integral A mediante una extrapolación de
Richardson según el método de Romberg. Entonces, partimos de las aproxima-
(1) h
ciones proporcionadas
" por la fórmula# de los Trapecios Th (f ) = 2 [f (0) + f (1)]
(1)
y Th/2 (f ) = h f (0) 1
2 +f 2 + 2
f (1)
siendo h = 1. Se tiene pues que
(1) 1 1
T1 (f ) = [f (0) + f (1)] = [0 + 5,55761] = 2,77881 ,
2
2
(1) 1 f (0) 1 f (1) 1 0 5,55761
T1/2 (f ) = +f + = + 0,89533 +
2 2 2 2 2 2 2
= 1,83707 .
Ecuaciones Diferenciales
Ordinarias
y = f (x, y) , y(x0 ) = y0 ,
tiene una única solución y(x) definida para |x − x0 | ≤ mı́n{α, β/M } , siendo M
el máximo de |f (x, y)| en el rectángulo R.
125
126 Ecuaciones Diferenciales Ordinarias
x0 x1 x2 ··· xn
y0 y1 y2 ··· yn
k1 = f (ti , xi ) ,
h h
k2 = f ti + , xi + k1 ,
2 2
h h
k3 = f ti + , xi + k2 ,
2 2
k4 = f (ti + h, xi + hk3 ) .
7.2 Problemas resueltos 127
Problema Lineal de Contorno: Si p(x), q(x), r(x) ∈ C[a, b] con q(x) > 0
en [a, b], entonces el problema lineal de frontera y (x) = p(x)y + q(x)y + r(x)
con y(a) = α, y(b) = β tiene una única solución y(x) con x ∈ [a, b]. Además
dicha solución viene dada por y(x) = u(x) + Cv(x), siendo C = (β − u(b))/v(b),
donde u(x) y v(x) son la única solución de los problemas de Cauchy (1) u (x) =
p(x)u + q(x)u + r(x) con u(a) = α, u (a) = 0; (2) v (x) = p(x)v + q(x)v con
v(a) = 0, v (a) = 1.
Realizando esta iteración con el valor inicial (x0 , y0 ) = (0, 0) obtenemos la tabla
siguiente
128 Ecuaciones Diferenciales Ordinarias
i xi yi
0 0 0
1 0.5 0.5
2 1 1.07452
3 1.5 1.91421
4 2 3.16173
donde la tercera columna es una tabulación aproximada de la función y(x).
Problema 7.2 La función de error erf(x), utilizada en estadı́stica y en proba-
bilidad, está definida de la forma
! x
2
erf(x) = √ exp(−t2 ) dt .
π 0
i xi yi
0 0 0
1 0.5 0.56419
2 1 1.00358
3 1.5 1.21113
4 2 1.2706
2 1 4
lı́m φ(xi , yi ; h) = f (xi , yi ) + f (xi , yi ) + f (xi , yi ) = f (xi , yi ) ,
h→0 9 3 9
y (xi )
Δ(xi , yi ; h) = f (xi , yi ) + h + O(h2 ) − φ(xi , yi ; h)
2!
y (xi ) 2 2 1 4
= k1 + h + O(h ) − k1 + k2 + k3
2! 9 3 9
7 y (xi ) 1 4
= k1 + h + O(h2 ) − k2 + k3 .
9 2! 3 9
130 Ecuaciones Diferenciales Ordinarias
con valor inicial (t0 , v0 ) = (0, 0). Se obtiene de esta forma la tabla siguiente
i ti vi
0 0 0
1 0.5 3.43333
2 1 7.00591
Solución. Sea y(x) la única solución del problema de valor inicial y = f (x, y)
con y(x0 ) = y0 . Como la fórmula de Adams-Moulton que hemos de demostrar
es de la forma yi+1 = yi + h(T ERM IN O), entonces separamos variables e
integramos la ecuación diferencial de la forma
! y(xi+1 ) ! xi+1
dy = f (x, y(x)) dx ,
y(xi ) xi
es decir ! xi+1
y(xi+1 ) = y(xi ) + f (x, y(x)) dx .
xi
xi+1 xi yi x0 1/2
= + , = .
yi+1 yi −xi − yi3 y0 1/4
i xi yi ti
0 1/2 1/4 0
1 0.75 -0.2656 1
2 0.484 -0.9968 2
3 -0.512 -0.4905 3
4 -1.003 0.1399 4
5 -0.863 1.1403 5
6 0.2772 0.5205 6
(i) Aproximar a partir de (7.9) el valor de erf(2) utilizando una vez la extrapo-
lación de Richardson según el método de Romberg.
(ii) Utilizar un paso del método clásico de Runge-Kutta de cuarto orden√para
aproximar el valor de x(t) e y(t) solución de (7.6) y (7.7) para t = 2 2 −
1 = 1.82843.
(iii) ¿Es posible comparar los resultados obtenidos en los anteriores apartados?
(1)
Solución. (i) Definamos la función integrando f (s) = exp(−s2 ). Sea Th (f ) la
regla de los trapecios aplicada a la función f con longitud de paso h. Entonces,
aplicando la regla de los trapecios simple se tiene
! 2
2 2 (1)
erf(2) = √ f (s) ds ≈ √ T2 (f )
π 0 π
2 2 2
= √ [f (0) + f (2)] = √ [1 + 0.0183156] = 1.14905 .
π2 π
Utilizando de nuevo la regla de los trapecios, pero ahora con longitud de paso
la mitad que la anterior, se obtiene
! 2
2 2 (1)
erf(2) = √ f (s) ds ≈ √ T1 (f )
π π
0
2 f (0) f (2)
= √ + f (1) +
π 2 2
2 1 0.0183156
= √ + 0.367879 + = 0.989631
π 2 2
7.2 Problemas resueltos 135
A partir de las dos estimaciones de erf(2) realizadas se obtiene una mejor aproxi-
mación por extrapolación de la forma
2 (2) 4 × 0.989631 - 1.14905
erf(2) ≈ √ T2 (f ) = = 0.936491 .
π 4−1
(ii) El método RK4 con una longitud de paso h aplicado a la solución numéri-
ca del problema de Cauchy
ẋ = f (t, x, y)
ẏ = g(t, x, y) ,
donde
r1 = f (ti , xi , yi ) ,
k1 = g(ti , xi , yi ) ,
h h h
r2 = f ti + , xi + r1 , yi + k1 ,
2 2 2
h h h
k2 = g ti + , xi + r1 , yi + k1 ,
2 2 2
h h h
r3 = f ti + , xi + r2 , yi + k2 ,
2 2 2
h h h
k3 = g ti + , xi + r2 , yi + k2 ,
2 2 2
r4 = f (ti + h, xi + hr3 , yi + hk3 ) ,
k4 = g (ti + h, xi + hr3 , yi + hk3 ) .
de modo que
1.82843 1.82843
r2 = f , 0.483941 + 0.198748 , 0.682689
2 2
1.82843
+ 0.483941 = 1.48807 ,
2
136 Ecuaciones Diferenciales Ordinarias
1.82843 1.82843
k2 = g , 0.483941 + 0.198748 , 0.682689
2 2
1.82843
+ 0.483941 = -0.362958 ,
2
y en consecuencia
1.82843 1.82843
r3 = f , 0.483941 + 1.48807 , 0.682689
2 2
1.82843
− 0.362958 = -1.17272 ,
2
1.82843 1.82843
k3 = g , 0.483941 + 1.48807 , 0.682689
2 2
1.82843
− 0.362958 = 1.52359 ,
2
de manera que
y está claro que se verifica que tanto f como ∂f /∂y son continuas en un entorno
de cualquier punto (x0 , y0 ) ∈ R2 y en particular en un entorno del punto (0, 1).
Se concluye pues que el problema de Cauchy y = x2 + y, y(0) = 1 tiene solución
única.
h2
y(x + h) = y(x) + hy (x) + y (x) + O(h3 ) .
2
Calculando las derivadas
∂f ∂f
y (x) = f (x, y(x)) , y (x) = (x, y(x)) + (x, y(x))f (x, y(x)) ,
∂x ∂y
se concluye que el método de Taylor para el problema de Cauchy del enunciado
consiste en la aplicación de la iteración siguiente
h2 " ∂f ∂f #
yi+1 = yi + hf (xi , yi ) + (xi , yi ) + (xi , yi )f (xi , yi ) , i = 0, 1, . . . .
2 ∂x ∂y
h2 " #
yi+1 = yi + h(x2i + yi ) + 2xi + x2i + yi , i = 0, 1, . . . .
2
La aplicación de esta iteración permite realizar la tabla siguiente
n xn yn
0 0.0 1.0
1 0.1 1.105
2 0.2 1.22307
3 0.3 1.35769
∗ 1
yi+1 = yi + yi cos xi , i = 0, 1, 2.
2
1 ∗
yi+1 = yi + [yi cos xi + yi+1 cos xi+1 ] ,
4
que da lugar a la tabulación
i xi yi
0 0 1
1 1/2 1.57909
2 1 2.23242
donde se han incorporado las condiciones iniciales del problema de Cauchy en los
lı́mites inferiores de las integrales. En definitiva, la solución general del problema
de Cauchy viene dada por y(x) = esin x .
Finalmente, el error absoluto Δ cometido en la aproximación y(1) ≈ y2 es
Δ = |y(1) − y2 | = esin 1 − 2,23242 = |2,31978 − 2,23242| = 0,08736 .
Problema 7.11 Obtener el error absoluto que comete el método de Euler cuan-
do aproxima y(1/2) mediante 2 iteraciones siendo y(x) la solución del problema
de Cauchy y = (x − 1)y 2 con la condición inicial y(0) = 2.
En nuestro caso se tiene (x0 , y0 ) = (0, 2). Como además se quiere aproximar
y(1/2) mediante 2 iteraciones, es decir, y(1/2) ≈ y2 , se concluye que la condición
x2 = x0 +2h debe ser 1/2 = 0+2h de modo que la longitud de paso vale h = 1/4.
De este modo, la iteración de Euler que se debe realizar es
1
yi+1 = yi + hf (xi , yi ) = yi + (xi − 1)yi2 .
4
La tabla que se obtiene es la siguiente.
i xi yi
0 0 2
1 1/4 1
2 1/2 0.8125
con la condición inicial (x0 , v0 ) = (1, 0). Se genera de esta forma la siguiente
sucesión
n tn xn vn
0 0 1 0
1 0.5 1 -2.32042
2 1 -0.160208 -4.65339
3 1.5 -2.4869 5.06518
4 2 0.0456857 7.39007
T ≈ t4 = 4h = 4 × 0,5 = 2 .
Problema 7.13 Un cohete tiene una velocidad v(t) respecto de un cierto sis-
tema de referencia inercial y sus gases de escape son expulsados a una veloci-
dad constante V respecto del cohete. Supóngase que 0 << v/c < 1, siendo
c = 3 × 108 m/s la velocidad de la luz en el vacı́o, es decir para describir el
movimiento del cohete es necesaria la teorı́a de la relatividad especial de Ein-
stein. En estas condiciones se puede demostrar que la ecuación del movimiento
del cohete es 3/2
dv dm0 v2
m0 +V 1− 2 =0,
dt dt c
donde m0 (t) es la masa variable del cohete en su sistema de referencia en reposo.
(i) Supongamos que se conoce la siguiente tabulación en el sistema interna-
cional de la función m0 (t)
ti 0 10 20 30 40
m0 (ti ) 1200 1100 1005 910 900
Calcular una tabulación aproximada de dm0 /dt para los tiempos ti = 10i,
con i = 0, 1, 2, 3, 4, utilizando de entre las fórmulas de diferencias hacia
adelante, hacia atrás y centradas, la más adecuada en cada caso.
(ii) Utilizando el resultado del apartado anterior, calcular mediante el método
de Euler con una longitud de paso adecuada, una aproximación de v(20)
sabiendo que V = 103 m/s y v(0) = 104 .
142 Ecuaciones Diferenciales Ordinarias
Solución. (i) Las fórmulas de diferencias hacia adelante, hacia atrás y centradas
son
dm0 m0 (ti + h) − m0 (ti )
(ti ) = + O(h) ,
dt h
dm0 m0 (ti ) − m0 (ti − h)
(ti ) = + O(h) ,
dt h
dm0 m0 (ti + h) − m0 (ti − h)
(ti ) = + O(h2 ) ,
dt 2h
respectivamente, siendo h la longitud de paso. Siempre que sea posible (es decir,
en todos los nodos interiores) utilizaremos la fórmula de diferencias centradas
puesto que es la que contiene un menor error de truncamiento. Para el primer
y último nodo sólo se pueden usar las fórmulas de diferencias hacia adelante y
hacia atrás respectivamente. En definitiva, tomando la longitud de paso h = 10,
se tiene
dm0 m0 (10) − m0 (0) 1100 − 1200
(0) ≈ = = −10 ,
dt 10 10
dm0 m0 (20) − m0 (0) 1005 − 1200 39
(10) ≈ = =− ,
dt 20 20 4
dm0 m0 (30) − m0 (10) 910 − 1100 19
(20) ≈ = =− ,
dt 20 20 2
dm0 m0 (40) − m0 (20) 900 − 1005 21
(30) ≈ = =− ,
dt 20 20 4
dm0 m0 (40) − m0 (30) 900 − 910
(40) ≈ = = −1 .
dt 10 10
(ii) La función v(t) es la solución del problema de Cauchy
dv
= f (t, v) , v(0) = 104 ,
dt
siendo
3/2
V dm0 v2
f (t, v) = − (t) 1 − 2 .
m0 (t) dt c
Teniendo en cuenta que c = 3 × 108 , V = 103 y los valores m0 (ti ) y dm 0
dt (ti )
para los tiempos ti = 10i, con i = 0, 1, 2, 3, 4, se pueden obtener del apartado
anterior, podemos aplicar el método de Euler
vi+1 = vi + hf (ti , vi ) , i = 0, 1, 2, 3 ,
3/2
103 dm0 (104 )2
= 104 − 10 × (0) 1 −
m0 (0) dt (3 × 108 )2
3/2
4 103 (104 )2
= 10 − 10 × × (−10) 1 − = 10083,3 ,
1200 (3 × 108 )2
v2 = v1 + 10f (t1 , v1 ) = 10083,3 + 10f (10, 10083,3)
3/2
103 dm0 (10083,3)2
= 10083,3 − 10 × (10) 1 −
m0 (10) dt (3 × 108 )2
3/2
103 39 (10083,3)2
= 10083,3 − 10 × − 1− = 10171,9
1100 4 (3 × 108 )2
k
v̇ = f (t, v) = −g − v = −9,8 − 2v
m
El problema de Cauchy está formado por la anterior ecuación diferencial acom-
pan̄ada de la condición inicial v(0) = v0 = 20. Aplicaremos el método de Heun
para aproximar la solución v(t). Dicho método consiste en realizar la iteración
144 Ecuaciones Diferenciales Ordinarias
predictor–corrector
∗ 1
vi+1 = vi + hf (ti , vi ) = vi + (−9,8 − 2vi ) ,
2
h& ∗
'
vi+1 = vi + f (ti , vi ) + f (ti+1 , vi+1 )
2
1& ∗
'
= vi + −19,6 − 2(vi + vi+1 ) ,
4
donde los tiempos ti son equiespaciados con longitud de paso h, es decir, ti =
ih = i/2 con i = 0, 1, 2, . . .. De este modo se obtienen las aproximaciones v(ti ) ≈
vi dadas en la siguiente tabla
i ti vi
0 0 20
1 0.5 7.55
2 1 1.325
3 1.5 -1.7875
de modo que
m k k
v(t) = −g + g + v0 exp − t .
k m m
Con los datos del problema se tiene
1
v(t) = [−9,8 + 49,8 exp (−2t)] ,
2
de modo que, la ecuación v(τ ) = 0 tiene por solución
m v0 k
τ= ln 1 + = 0,8128 .
k mg
El error absoluto cometido en el apartado anterior es
Δ = |0,8128 − 1,25| = 0,4372 .
7.2 Problemas resueltos 145
Problema 7.15 Una viga horizontal, colocada a lo largo del eje x se pandea
bajo la influencia de diversas cargas verticales. Según la teorı́a de la elasticidad,
la curva de flexión y(x) se puede determinar, para pequen̄as flexiones, mediante
la ecuación diferencial lineal de segundo orden
d2 y
EI = M (x) ,
dx2
siendo M (x) el momento flector en x, E el módulo de Young del material e I
el momento de inercia de una sección transversal de la viga respecto a su eje
central.
(i) Se tiene una viga con un extremo (situado en x = 0) empotrado en una
dy
pared, es decir, y(0) = dx (0) = 0 y tal que M (x) = 10(1 − x) con EI = 2.
Aproximar la flexión de la viga en x = 2 mediante el método de Runge–
Kutta de orden 4 con una longitud de paso h = 1.
(ii) Hallar la flexión exacta y(x) de la viga del apartado anterior. ¿Se cal-
culó una buena aproximación en el apartado (i)?
Solución. (i) Con los datos del problema, la flexión y(x) de la viga verifica la
ecuación diferencial
d2 y
= 5(1 − x) ,
dx2
dy
con las condiciones iniciales y(0) = dx (0) = 0. Realizando el cambio de variable
dy
usual dx = z, reescribimos el anterior problema de Cauchy de segundo orden
como un problema de Cauchy de primer orden de la forma
k1 = f (xi , wi ) ,
h h
k2 = f xi + , wi + k1 ,
2 2
h h
k3 = f xi + , wi + k2 ,
2 2
k4 = f (xi + h, wi + hk3 ) .
146 Ecuaciones Diferenciales Ordinarias
i xi wi = (yi , zi )
0 0 (0,0)
1 1 (5/3 , 5/2)
2 2 (10/3 , 0)
d2 y
= 5(1 − x) ,
dx2
dy
con las condiciones iniciales y(0) = dx (0) = 0. Podemos hallar la solución y(x)
de manera exacta integrando dos veces la ecuación diferencial respecto de x.
Integrando una vez se tiene que
!
dy x2
=5 (1 − x) dx + C1 = 5 x − + C1 ,
dx 2
!
2
x2 x x3
y(x) = 5 x− dx + C2 = 5 − + C2 ,
2 2 6
siendo C2 una constante arbitraria que calculamos con la condición inicial y(0) =
0, es decir, C2 = 0. En resumen, la flexión exacta de la viga es
5
y(x) = (3 − x)x2 ,
6
de modo que y(2) = 10/3. Entonces, el método RK4 del apartado (i) calculó y(2)
de manera exacta, es decir, con un error absoluto cero.
7.2 Problemas resueltos 147
Q
LQ̈ + RQ̇ + = V (t) .
C
Si L = R = 2, C = 1 y V (t) = 3 sin(t), aproximar mediante el método de Euler
con una longitud de paso h = 1 la carga del condensador para tiempo t = 2,
sabiendo que inicialmente tenı́a una carga Q(0) = 3. Recordar que la intensidad
de corriente I(t) que circula es, por definción, I = Q̇ siendo I(0) = 0.
Solución. Con los datos del enunciado, la carga eléctrica Q(t) almacenada en el
condensador satisface la ecuación diferencial de segundo orden 2Q̈ + 2Q̇ + Q =
3 sin(t). Esta ecuación se escribe como un sistema de primer orden de la forma
1
Q̇ = I , I˙ = [3 sin(t) − 2I − Q] .
2
Se satisface la condición inicial (Q(0), I(0)) = (3, 0). En definitiva, definiendo el
vector x = (Q, I) ∈ R2 y la función f : R3 → R2 de la forma
1
f (t, x) = I, [3 sin(t) − 2I − Q] ,
2
Utilizando el método de Euler con una longitud de paso h = 1, se tiene que las
aproximaciones x(tn ) = (Q(tn ), I(tn )) ≈ xn = (Qn , In ) ∈ R2 para los tiempos
equiespaciados tn = 0 + nh = n con n ∈ N vienen dadas de forma recurrente
por la iteración xn+1 = xn + hf (tn , xn ). Se tiene pues
Qn+1 Qn In Q0 3
= + 1 , = .
In+1 In 2 [3 sin(tn ) − 2In − Qn ] I0 0
n tn (Qn , In )
0 0 (3,0)
1 1 (3, -3/2)
2 2 (1.5, -0.237794)
xx2 x
y (x) − f (x, y(x)) = 2x − 2 2 4
= 2x − = 0 .
3x − (x ) 3 − x6
(iii) Sea y(x) la solución del problema de Cauchy y = f (x, y), y(x0 ) = y0 .
Discretizando equiespaciadamente la variable independiente con una longitud de
paso h, definimos xi = x0 + ih con i = 0, 1, 2, . . .. Si aproximamos y(xi ) ≈ yi , el
7.2 Problemas resueltos 149
i xi yi∗ yi
0 1 1 1
1 1.5 1.25 1.23379
2 2 1.44255 1.43221
Solución. Teniendo en cuenta que la posición en función del tiempo del punto
de suspensión es (x(t), y(t)) = (t2 + t, 7t), su aceleración es (ẍ(t), ÿ(t)) = (2, 0).
Entonces, introduciendo el dato = 1, la función ϕ(t) satisface el problema de
Cauchy de segundo orden
i ti (ϕi , wi )
0 0 (3,1)
1 1 (4, 1.59701)
2 2 (5.59701, 10.321)
Solución de Ecuaciones no
Lineales
Teorema 8.1 (Bolzano) Sea f una función continua en [a, b] tal que cambia
de signo en los extremos, es decir, f (a)f (b) < 0. Entonces, existe al menos un
ξ ∈ (a, b) tal que f (ξ) = 0.
El algoritmo de bisección se basa en el teorema de Bolzano y consiste en
ir estrechando de manera sistemática el intervalo en el cual una función con-
tinua cambia de signo. De esta forma conseguiremos obtener un intervalo ar-
bitrariamente pequen̄o que contenga el cero de la función. El método de bisec-
ción para una función f ∈ C[a0 , b0 ] procede de la forma siguiente: (a0 , b0 ) ⊃
(a1 , b1 ) ⊃ · · · ⊃ (an , bn ) ⊃ · · ·. La sucesión de puntos medios {cn }∞
n=0 siendo
cn = (an + bn )/2, verifica lı́mn→∞ cn = x∗ con |x∗ − cn | ≤ (b − a)/2n+1 .
Método de Newton-Raphson:
f (xn )
xn+1 = xn − .
f (xn )
153
154 Solución de Ecuaciones no Lineales
Método de la Secante:
xn − xn−1
xn+1 = xn − f (xn ) .
f (xn ) − f (xn−1 )
x = (x1 , x2 , . . . , xn )T .
(1/2)8 1
|x∗ − x8 | ≤ |x1 − x0 | = |x1 − x0 | = 0,0078125|x1 − x0 |.
1 − 1/2 128
P (xn ) x5 + 5xn + 8
xn+1 = xn −
= xn − n 4
P (xn ) 5(xn + 1)
tomando x0 = −2 se obtiene
n xn
0 -2
1 -1.58025
2 -1.27709
cosh x - cos x
10
x
0.5 1 1.5 2 2.5 3
Solución. (i) Una condición necesaria para que la función g(x) tenga un ex-
tremo relativo en x∗ es que la función f (x) := g (x) se anule para x = x∗ . Note-
mos que f (x) = cosh x+cos x−3 es una función continua en todo R y por lo tanto
en particular f ∈ C[1, 2]. Como además f (1)f (2) = (−0,916617) × 0,346049 < 0
se puede aplicar el Teorema de Bolzano de manera que se concluye con la existen-
cia de al menos un valor x∗ ∈ (1, 2) tal que f (x∗ ) = 0. Una condición suficiente
para que x∗ sea un extremo relativo de g es que g (x∗ ) = 0. Pero esta condición
está garantizada puesto que si las curvas y = sinh x, y = sin x no tiene ningún
punto de corte en el intervalo [1, 2] entonces la función g (x) = sinh x − sin x no
8.2 Problemas resueltos 157
(ii) El extremo relativo x∗ de g(x) viene dado por la raı́z de la función f (x) en
el intervalo [1, 2]. Recordemos que si f ∈ C 2 [1, 2] y además f (x) tiene signo con-
stante para toda x ∈ [1, 2] entonces la recta tangente trazada desde el extremo
del intervalo [1, 2] adecuada para hallar la raı́z de f por el método de Newton-
Raphson es el extremo p (con p = 1 o p = 2) en el cual sigf (p) = sigf (p). A
partir de la gráfica adjunta en el enunciado se comprueba que f (x) tiene signo
constante en [1, 2]. Además, como f (2) = 0,3460 > 0 y f (2) = 4,1783 > 0
resulta p = 2.
n xn f (xn )
0 2 0.346049
1 1.87266 0.0323507
2 1.8581 0.000382472
dicho punto, es decir, ∇u(x, y) := (ux (x, y), uy (x, y)) = (0, 0). Calculando las
derivadas paciales de u se tiene
0 −1
D∇u (x∗2 , y2∗ ) = ,
−1 2
0 1
D∇u (x∗3 , y3∗ ) = ,
1 0
2 1
D∇u (x∗4 , y4∗ ) = .
1 2
30
20
10
r
1 2 3 4
-10
-20
Figura 8.1: Gráfica de f (r) = [r2 T (r)] para los valores 0 ≤ r ≤ 4 cuando
T (r) = (r − R) exp(sin r).
En definitiva, se tiene
27 53
u(2) ≈ − (−16,8896) + (−37,184) = −0,116041 .
316 1264
Por otra parte, se sabe que el error absoluto Δ cometido al utilizar la regla
b 2
de los Trapecios viene dado por Δ := a f (r) dr − Th (f ) = − h12 (b − a)f (ξ)
con ξ ∈ (a, b). En nuestro problema hemos realizado dos aproximaciones y por
lo tanto se tienen dos errores
! 2
1
Δ1 := f (r) dr − T1 (f ) = |f (ξ1 )| ,
0 6
! 4
1
Δ2 := f (r) dr − T1 (f ) = |f (ξ2 )| ,
0 3
n rn
0 0
1 1.42881
2 2.14499
3 2.34559
Demostrar que existe al menos una raı́z ∗ positiva de la ecuación (8.5) es equi-
valente a demostrar que existe ∗ > 0 tal que f (∗ ) = 0. Para ello nos basaremos
en el Teorema de Bolzano. En primer lugar, observemos que, la expresión del
tiempo T dada por (8.5) sólo tiene sentido fı́sico si ≥ b > 0. Matemáticamente
observamos que si no se verifica la desigualdad ≥ b, debido a la raı́z cuadrada
que aparece en (8.5), la función T no es real.
De esta forma, sólo tiene sentido que estudiemos la función f () para > 5.
Además, bajo esta restricción, los únicos problemas de continuidad de la fun-
ción f , es decir, la raı́z cuadarada y el logaritmo, desparecen ya que sus ar-
gumentos son mayor o igual que cero y positivo respectı́vamente. En definitiva
la función f ∈ C[5, ∞). En particular f ∈ C[5, 150] y verifica f (5)f (150) =
15 × (−1,01723) < 0, de manera que, aplicando el Teorema de Bolzano con-
cluimos que existe ∗ ∈ (5, 150) tal que f (∗ ) = 0.
log 106 6
n > log2 10−6 − 1 = −1= − 1 = 18,93 .
log 2 log 2
Se puede asegurar entonces que, con 19 iteraciones del método de bisección
obtenemos la longitud ∗ con un error menor que 10−5 .
8.2 Problemas resueltos 167
f tiene un único mı́nimo relativo en x = 3/5 = 0,6 de manera que f (x) < 0
para todo x ∈ (0, 0,6) y f (x) > 0 para todo x ∈ (0,6, ∞). Además f ∈ C(0, ∞),
f (0,6) < 0, lı́mx→0+ f (x) = ∞ y f (1) = 0, de manera que, ver Figura 8.2, se ha
demostrado que f tiene un único cero x̃1 ∈ (0, 1).
Otra forma de proceder es la siguiente. Utilizando el Teorema de Bolzano se
puede ver que f tiene al menos un cero x̃1 ∈ (0, 1). Efectivamente, puesto que
170 Solución de Ecuaciones no Lineales
f
8
x
0.5 1 1.5 2
f es una función continua para todo x > 0, en particular f ∈ C[0,1, 0,5] por
ejemplo. Además f (0,1)f (0,5) = f (0,1)f (0,5) < 0. Entonces, se concluye que
f tiene al menos un cero x̃1 ∈ [0,1, 0,5]. Para ver que dicho cero es único en
este intervalo la función f es decreciente y por lo tanto f tiene un único cero en
x̃1 ∈ [0,1, 0,5].
(ii.2) Una vez se sabe que f tiene un único cero x̃1 ∈ [0,1, 0,5], se utilizará el
método de Newton-Raphson con el objetivo de aproximarlo, es decir
Como f ∈ C 2 [0,1, 0,5], f (0,1)f (0,5) < 0 y tanto f como f tienen signo con-
stante en el intervalo [0,1, 0,5], se tomará, con el objetivo de asegurar la con-
vergencia de la sucesión {xn }n∈N hacia x̃1 , el punto inicial x0 de la iteración de
Newton-Raphson verificando x0 ∈ {0,1, 0,5} tal que sign(f (x0 )) = sign(f (x0 )).
Es evidente que el valor inicial es x0 = 0,1. De esta forma, la sucesión {xn }n∈N
generada mediante la iteración de Newton-Raphson (8.13) con el valor inicial
x0 = 0,1 tiene los siguientes primeros términos de la sucesión.
n xn
0 0.1
1 0.19631
2 0.280509
3 0.318437
(ii.3) El periodo T dado por la fórmula integral (8.7) es una integral impropia
de segunda especie ya que la función integrando
1
f (x) = ,
E − U (x)
n tn xn
0 0 0.9
1 0.1 0.853091
8π hc
u(λ) = 5
hc ,
λ exp λKT − 1
(i) Demostar que la longitud de onda λmax para la cual u(λ) es máxima se puede
escribir como λmax = αhc/(KT ), siendo α ∈ R una raı́z de la ecuación
trascendente
e1/α (5α − 1) − 5α = 0 . (8.14)
(ii) Demostrar que la ecuación (8.14) tiene una única raı́z en el intervalo
[0,15, 0,25].
(iii) Aproximar el valor de α mediante dos iteraciones del método de Newton-
Raphson, partiendo de un punto inicial tal que nos asegure la convergencia
hacia la raı́z buscada.
(iv) Supongamos que la temperatura del cuerpo negro es T = 350 grados
Kelvin. Aproximar con qué exactitud se puede calcular la energı́a radi-
ada por unidad de volumen u sabiendo que la longitud de onda vale λ =
5 × 10−6 ± 10−7 metros. Datos: En el sistema internacional c = 3 × 108 ,
K = 1,381 × 10−23 y h = 6,626 × 10−34 .
Solución. (i) Puesto que u(λ) es máxima para λ = λmax se debe verificar que
la derivada du/dλ se anule para λ = λmax . La primera derivada de la función u
respecto de λ es
&
hc '
du 8chπ exp λKT (hc − 5KT λ) + 5KT λ
= &
'2 ,
dλ KT λ7 exp hc − 1
λKT
f (αn )
αn+1 = αn − , n = 0, 1, 2, . . .
f (αn )
n αn
0 0.15
1 0.165583
2 0.181035
30
20
10
-6 -4 -2 2 4 6
-10
-20
-30
n cn f (cn )
0 1.25 -0.259
1 1.125 0.782
2 1.187 0.332
176 Solución de Ecuaciones no Lineales
n xn
0 1
1 0.52134
2 0.44988
ẋ = f (t, x, y) = y ,
ẏ = g(t, x, y) = −y − 2x + F (t) = −y − 2x + t sin t .
xi+1 = xi + hf (ti , xi , yi ),
yi+1 = yi + hg(ti , xi , yi ),
xi+1 = x i + yi ,
yi+1 = −2xi + i sin i,
i ti xi yi
0 0 5 0
1 1 5 -10
2 2 -5 -9.15853
Problema 8.13 Los laureados con el premio Nobel de fı́sica E. Fermi, J. Pasta
y S. Ulam descubrieron mediante experimentos numéricos con un ordenador que
un sistema de osciladores no lineales acoplados no verificaban una propiedad
de equilibrio termodinámico llamada equipartición de la energı́a. Posteriormen-
te, M. Toda fue capaz de estudiar el problema de forma analı́tica mediante la
introducción de un potencial de interacción del tipo
α
V (r) = exp(−βr) + αr .
β
Tomando α = 2 y β = 1:
(i) Demostrar que existe un único r∗ > 0 verificando V (r∗ ) = 4. ¿Es posible
hallar r∗ de manera exacta?
(ii) Aplicar 2 iteraciones del algoritmo de bisección con el objetivo de aproximar
r∗ . Iniciar el algoritmo en un intervalo del tipo [m, m+1] ⊂ R con m ∈ N.
(iii) Partiendo del intervalo de menor anchura obtenido en el apartado anterior
tras 2 iteraciones del algoritmo de bisección, refinar la aproximación de r∗
mediante 2 iteraciones del método de Newton-Raphson. Iniciar la iteración
en un punto adecuado que asegure la convergencia hacia r∗ de la sucesión
de Newton-Raphson.
(iv) Hallar una aproximación del error relativo que se comete al evaluar la
función V (r) cuando r = 10 ± 1.
es inmediato comprobar que f (r) > 0 para todo r > 0. Por lo tanto la función
f (r) es creciente en el intervalo (0, ∞). Como además verifica que
(ii) Observemos en primer lugar que la función f (r) ∈ C(R). Entonces, el in-
tervalo [m, m+1] ⊂ R con m ∈ N para iniciar el algoritmo de bisección sólo debe
de verificar la condición f (m)f (m + 1) < 0. Se halla calculando el signo de f (m)
con m = 0, 1, 2, . . . y parando cuando se observe un cambio de signo. En concreto
se tiene f (0) = −2 < 0, f (1) = −1,26424 < 0 y f (2) = 0,270671 > 0 por lo tanto
el intervalo inicial es [a0 , b0 ] = [m, m+1] = [1, 2]. Iniciamos el algoritmo calculan-
do el punto medio c0 = (1 + 2)/2 = 1,5. Como f (c0 ) = f (1,5) = −0,55374 < 0,
entonces el nuevo intervalo es [a1 , b1 ] = [c0 , b0 ] = [1,5, 2] y su punto medio
c1 = (1,5 + 2)/2 = 1,75. Como f (c1 ) = f (1,75) = −0,152452 < 0, entonces el
siguiente intervalo es [a2 , b2 ] = [c1 , b1 ] = [1,75, 2] con punto medio c2 = 1,875.
En definitiva, r∗ ≈ c2 = 1,875.
(iii) Según los resultados anteriores, f ∈ C 2 [1,75, 2], f (1,75)f (2) < 0 y
f (r) > 0 para todo r > 0. Además, como
f (r) = 2 exp(−r) ,
la función f (r) tiene curvatura positiva en todo R puesto que f (r) > 0 para
todo r ∈ R. En definitiva, como el signo de f (2) coincide con el signo de f (2),
tomamos el punto inicial r0 = 2 en la iteración de Newton-Raphson
f (rn ) exp(−rn ) + rn − 2
rn+1 = rn − = rn − .
f (rn ) 1 − exp(−rn )
De esta forma se obtiene la siguiente tabla
n rn
0 2
1 1.8434
2 1.8414
h2
Δ= |f (ξ)| , ξ ∈ (0, 1) .
12
8.2 Problemas resueltos 181
6ξ 2 (10 + ξ 4 )
f (ξ) =
(10 − ξ 4 )5/2
y además ξ ∈ (0, 1), es obvio que f (ξ) > 0. Además, al aumentar ξ aumenta
el numerador y disminuye el denominador de f (ξ) con lo que vemos que la
función f (ξ) es creciente1 en el intervalo ξ ∈ (0, 1). En definitiva, se puede
acotar el error absoluto de la forma
h2 h2 22 11 2
Δ< f (1) = = h .
12 12 81 486
Si imponemos que la longitud de paso h verifique la desigualdad
11 2
h < 10−4
486
nos aseguramos de que Δ < 10−4 como se desea. De esta forma se tiene que
486
h < 10−4 = 0,066 .
11
(ii) En primer lugar observamos que los puntos de retorno x∗ son los ceros
de la función f (x) = U (x) − E. En nuestro caso se tiene f (x) = x2 ex − 0,4.
(ii.1) Hemos de demostrar que la función f (x) tiene un único cero positivo
x∗ . Para verlo, notemos que f (0) < 0 y lı́mx→∞ f (x) = ∞ de modo que, como
la función f (x) es continua, esta tiene al menos un cero x∗ > 0. Para ver que x∗
es único, basta con demostrar que f (x) es una función creciente en el intervalo
(0, ∞). Esto es cierto puesto que f (x) = x(2 + x)ex > 0 para todo x > 0.
(ii.2) Como f ∈ C[0, 1], f (0) = −0,4 < 0 y f (1) = 2,31828 > 0, aplicando
el Teorema de Bolzano sabemos que existe x∗ ∈ [0, 1] tal que f (x∗ ) = 0. De
hecho, por el apartado anterior, sabemos que x∗ es único. Aplicando 3 iteraciones
en el método de bisección, obtendremos una colección de intervalos anidados
(a0 , b0 ) = (0, 1) ⊃ (a2 , b2 ) ⊃ (a3 , b3 ) de anchura del siguiente la mitad que la del
anterior tal que [ai , bi ] x∗ ≈ c3 , siendo ci el punto medio del intervalo i–ésimo.
Un simple cálculo nos muestra la siguiente tabla
i ci f (ci )
0 0.5 0.0121803
1 0.25 -0.319748
2 0.375 -0.195392
3 0.4375 -0.103544
1 Otraforma de ver que f (ξ) es creciente en el intervalo ξ ∈ (0, 1) es calculando f (ξ) =
24x(50 + 35x4 + x8 )/(10 − x4 )7/2 y viendo que, efectivamente, f (ξ) > 0 para todo ξ ∈ (0, 1).
182 Solución de Ecuaciones no Lineales
x2 − y = a , x2 + 4y 2 = 4 .
1.5
0.5
-2 -1 1 2
-0.5
-1
1.5
0.5
-2 -1 1 2
-0.5
-1
i xi yi
0 1 1
1 0.9444 0.8888
2 0.9395 0.8828
0 = f (x) = −4 + x2 + 4x4 .
del método. Para elejir x0 , observemos que f (x) = 2(1 + 24x2 ) > 0 para todo
x ∈ R. En definitiva, aplicando el Teorema de convergencia de la sucesión de
Newton–Raphson, se tiene que, puesto que i) f ∈ C 2 [1/2, 1]; ii) f (1/2)f (1) < 0;
iii) f (x) > 0 y f (x) > 0 para todo x ∈ [1/2, 1], concluimos que una elección
de x0 que asegura la convergencia de la sucesión de Newton–Raphson hacia x∗
es x0 = 1 ya que el signo de f (1) y de f (1) coinciden.
Realizaremos dos iteraciones del método de Newton-Raphson para aproxi-
mar x∗ . La iteración es
f (xi ) −4 + x2i + 4x4i
xi+1 = xi − = xi − ,
f (xi ) 2xi (1 + 8x2i )
i xi f (xi )
0 1 1
1 0.9444 0.0737
2 0.9396 0.0005
-3 -2 -1 1 2
f (xn )
xn+1 = xn − ,
f (xn )
y f (x) > 0. Entonces, tomando x0 ∈ {−2, −1} que verifique signf (x0 ) =
signf (x0 ) se asegura la convergencia de la sucesión de Newton–Raphson hacia
la única raı́z de f en [−2, −1]. En nuestro caso es claro que x0 = −2. Realizando
2 iteraciones de Newton-Raphson se obtiene la siguiente tabla
n xn f (xn )
0 -2 0.135335
1 -1.84348 0.00174769
2 -1.84141 3.40602 ×10−7
500
400
300
200
100
2 4 6 8 10 12
100
Figura 8.7: Gráfica del polinomio cúbico f () = (12 − 2)(20 − 2) − 49,5.
n n
0 0
1 0.20625
2 0.218801
La derivada de f es
f (r) = C(n + 1)rn − Cnrn−1 − anrn−1 = rn−1 [C(n + 1)r − n(C + a)]
= 5000r99 (202r − 203) .
de modo que los puntos de inflexión de f (r) son los ceros de f (r), es decir,
r = 0 y r = 20097/20200 ≈ 0,994901. Tras el análisis efectuado, es fácil ver
que un intervalo [α, β] ⊂ R tal que r∗ ∈ [α, β] y además el signo de f (r) y
de f (r) sea constante en todo [α, β] debe ser de la forma α > 203/202 con
f (α) < 0 y f (β) > 0. Además, para asegurar la convergencia de la sucesión
{rn }n∈N generada por el método de Newton-Raphson hacia r∗ , se debe tomar el
valor inicial r0 = β puesto que el signo de f (β) coincide con el signo de f (β).
Calculando vemos que f (2) > 0 de manera que tomaremos r0 = 2. Entonces, se
obtiene la siguiente tabla.
i ri
0 2
1 1.9804
2 1.96099
3 2 1 1 2 3
2
4
6
8
√
al menos un cero en el intervalo [1, 3]. Finalmente mostraremos√la unicidad de
dicha raı́z puesto que f (x) = 2x − sin x√= 0 para todo x ∈ [1, 3]. De hecho,
se tiene que f (x) > 1 para todo x ∈ [1, 3] debido a que −1 ≤ sin x ≤ 1.
√
(iv) Una vez sabemos que existe un único x∗ ∈ [1, 3] tal que f (x∗ ) = 0,
vamos a aproximarlo mediante dos iteraciones del método de Newton–Raphson
xn+1 = xn − ff(x n)
(xn ) . Para asegurar la convergencia de la sucesión {xn }n∈N ,
hemos
√ de elegir adecuadamente el punto inicial x0 . Tomamos el intervalo [a, b] =
[1, 3]. Por los cálculos realizados en el apartado (iii) sabemos que f ∈ C 2 [a, b]
con f (a)f (b) < 0. Además, f (x) > 0 para todo x ∈ [a, b]. Sólo falta estudiar
el signo de f (x) en el intervalo [a, b]. Puesto que f (x) = 2 − cos x, es claro
que f (x) > 0 para todo x ∈ [a, b]. De este modo, se tomará el valor inicial
x0 ∈ {a, b} que verifique signf (x0 ) = signf (x0 ). En definitiva, x0 = a = 1.
Realizando la iteración de Newton–Raphson
n xn f (xn )
0 1 -0.459698
1 1.39679 0.12416
2 1.32815 0.00425146
Aproximación de Funciones
f − f ∗ = ı́nf f − g .
g∈G
195
196 Aproximación de Funciones
ti (s) 0 1 2 3
T (ti ) (K) 27 34 25 38
siendo ϕ1 (t) = 1, ϕ2 (t) = sin t y ϕ3 (t) = cos2 t. Las ecuaciones normales para
obtener los parámetros α, β y γ según la teorı́a de mı́nimos cuadrados son
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
< ϕ 1 , ϕ1 > < ϕ 1 , ϕ2 > < ϕ 1 , ϕ3 > α < T, ϕ1 >
⎝ < ϕ2 , ϕ1 > < ϕ2 , ϕ2 > < ϕ2 , ϕ3 > ⎠ ⎝ β ⎠ = ⎝ < T, ϕ2 > ⎠ ,
< ϕ 3 , ϕ1 > < ϕ 3 , ϕ2 > < ϕ 3 , ϕ3 > γ < T, ϕ3 >
donde <, > denota el producto escalar discreto euclideano. Se tiene pues que
4
4
< ϕ 1 , ϕ1 > = ϕ21 (ti ) = 1=4,
i=1 i=1
198 Aproximación de Funciones
4
4
< ϕ 2 , ϕ2 > = ϕ22 (ti ) = sin2 (ti ) = 0 + 0,708073 + 0,826822
i=1 i=1
+0,0199149 = 1,55481 ,
4 4
< ϕ3 , ϕ 3 > = ϕ23 (ti ) = cos4 (ti ) = 1 + 0,0852211 + 0,0299907
i=1 i=1
+ 0,960567 = 2,07578 ,
4 4
< ϕ 1 , ϕ2 > = ϕ1 (ti )ϕ2 (ti ) = sin(ti ) = 0 + 0,841471 + 0,909297
i=1 i=1
+0,14112 = 1,89189 ,
4 4
< ϕ 1 , ϕ3 > = ϕ1 (ti )ϕ3 (ti ) = cos2 (ti ) = 1 + 0,291927 + 0,173178
i=1 i=1
+0,980085 = 2,44519 ,
4 4
< ϕ 2 , ϕ3 > = ϕ2 (ti )ϕ3 (ti ) = sin(ti ) cos2 (ti ) = 0,541428 ,
i=1 i=1
4
< T, ϕ1 > = T (ti )ϕ1 (ti ) = 27 + 34 + 25 + 38 = 124 ,
i=1
4
< T, ϕ2 > = T (ti )ϕ2 (ti ) = 56,705 ,
i=1
4
< T, ϕ3 > = T (ti )ϕ3 (ti ) = 78,4982 .
i=1
tiene
⎛ ⎞ ⎛ ⎞
1 sin 0 cos2 0 1 0 1
⎜ 1 sin 1 cos2 1 ⎟ ⎜ 1 0,841471 0,291927 ⎟
A=⎜ ⎟ ⎜ ⎟ ,
⎝ 1 sin 2 cos2 2 ⎠ = ⎝ 1 0,909297 0,173178 ⎠
1 sin 3 cos2 3 1 0,14112 0,980085
de manera que las ecuaciones normales adoptan la forma
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
4 1,89189 2,44519 α 124
⎝ 1,89189 1,55481 0,541428 ⎠ ⎝ β ⎠ = ⎝ 56,705 ⎠ ,
2,44519 0,541428 2,07578 γ 78,4982
siendo su solución es (α, β, γ) = (−86,7386, 102,583, 113,234).
Observar que los dos métodos utilizados llevan a las mismas ecuaciones nor-
males, aunque debido a los errores de redondeo en los cálculos efectuados exista
algún dı́gito diferente en ellas.
0 27
7
1 34 -8
-9
2 25
El polinomio interpolador es
P (t) = T [t0 ] + T [t0 , t1 ](t − t0 ) + T [t0 , t1 , t2 ](t − t0 )(t − t1 ) = 27 + 7t − 8t(t − 1) ,
que desarrollado resulta P (t) = 27 + 15t − 8t2 .
siendo ξ ∈ (0, 2). Como tanto t como ti pertenecen al intervalo [0, 2] se tiene que
2
|t − ti | ≤ 2 y por lo tanto i=0 |t − ti | ≤ 23 = 8. Además el enunciado asegura
que |T (3) (t)| ≤ 0,1 para todo tiempo t ∈ [0, 2] de manera que |T (3) (ξ)| ≤ 0,1.
Con estas acotaciones se tiene que
0,1
Δ = |T (t) − P (t)| ≤ 8 = 0,13333 .
3!
(iic) La aproximación requerida es
! 2 ! 2 ! 2
188
T (t) dt ≈ P (t) dt = [27 + 7t − 8t(t − 1)] dt = ≈ 62,66 .
0 0 0 3
200 Aproximación de Funciones
Tk (◦ C) 1 10 38
μk (centipoises) 50 10 5
ln μ(T ) = ln μ0 + k ln T . (9.2)
ln Tk 0 2.30259 3.63759
ln μk 3.91202 2.30259 1.60944
sistema lineal por la matriz traspuesta de los coeficientes del sistema. De esta
forma se obtiene el sistema lineal de ecuaciones normales
3 5,94018 log μ0 7,82405
= , (9.3)
5,94018 18,534 k 11,1564
cuya solución es (log μ0 , k) = (3,87571, −0,64023). Por lo tanto
μ0 = exp(3,87571) = 48,217 ,
y, en definitiva, la mejor aproximación por mı́nimos cuadrados viene dada por
μ(T ) = 48,217 T −0,64023 .
Existe otra forma de resolver el problema, basada en la obtención de las
ecuaciones normales (9.3) a partir del producto escalar euclideano discreto que
denotaremos por < ., . >. En concreto, tomando la función f (T ) = ln μ(T ) y las
funciones ϕ1 (T ) = 1, ϕ2 (T ) = ln T , se tiene que f (T ) es una combinación lineal
de ϕ1 (T ), ϕ2 (T ), es decir, f (T ) = ln μ0 ϕ1 (T ) + kϕ2 (T ). De esta forma, se sabe
que las ecuaciones normales adoptan la forma
< ϕ1 , ϕ1 > < ϕ1 , ϕ2 > log μ0 < f, ϕ1 >
= . (9.4)
< ϕ2 , ϕ1 > < ϕ2 , ϕ2 > k < f, ϕ2 >
Calculando los productos escalares euclideanos discretos asociados a la tabla del
enunciado e involucrados en este sistema se tiene
2
< ϕ1 , ϕ1 > = ϕ21 (Tk ) = 12 + 12 + 12 = 3 ,
k=0
2
< ϕ2 , ϕ2 > = ϕ22 (Tk )
k=0
= ln2 T0 + ln2 T1 + ln2 T2 = ln2 1 + ln2 10 + ln2 38 = 18,534 ,
2
< ϕ 1 , ϕ2 > = ϕ1 (Tk )ϕ2 (Tk ) = ln 1 + ln 10 + ln 38 = 5,94018 ,
k=0
2
< f, ϕ1 > = f (Tk )ϕ1 (Tk )
k=0
= ln μ(1) + ln μ(10) + ln μ(3) = ln 50 + ln 10 + ln 5
= 7,82405 ,
2
< f, ϕ2 > = f (Tk )ϕ2 (Tk ) = ln 50 ln 1 + ln 10 ln 10 + ln 5 ln 38
k=0
= 11,1564 .
De estos cálculos se desprende que las ecuaciones normales (9.4) coinciden con
las ecuaciones (9.3) y por lo tanto el problema sigue de igual manera.
202 Aproximación de Funciones
T − T1 T − T2 T − 10 T − 38
0 (T ) = = ,
T0 − T1 T0 − T 2 1 − 10 1 − 38
T − T0 T − T2 T − 1 T − 38
1 (T ) = = ,
T1 − T 0 T1 − T 2 10 − 1 10 − 38
T − T0 T − T1 T − 1 T − 10
2 (T ) = = .
T2 − T 0 T2 − T1 38 − 1 38 − 10
lleva los valores x ∈ [−1, 1] a valores T ∈ [1, 38]. Por lo tanto, los nodos que
tomaremos son
39 37 39 37 2k + 1
Tk = + xk = + cos π , k = 0, 1, . . . , n .
2 2 2 2 2n + 2
Finalmente, como n = 2, las temperaturas Tk donde se realizarán las mediciones
de la viscosidad serán
39 37 π 39 37 √3 √
78 + 37 3
T0 = + cos = + = = 35,5215◦ C ,
2 2 6 2 2 2 4
39 37 π 39
T1 = + cos = = 19,5◦ C ,
2 2 2 2 √ √
39 37 5π 39 37 3 78 − 37 3
T2 = + cos = − = = 3,47853◦ C .
2 2 6 2 2 2 4
Figura 9.1: Gráfica de la intensidad radioactiva I(t) en función del tiempo su-
perpuesta a la nube de puntos experimentales.
ti 0 1 2 3 4
x(ti ) 1.1 0.3 -6.5 -21 -39
Solución. Imponiendo que la función x(t) = x0 + v0 t + 12 at2 pase por todos los
puntos de la tabulación del enunciado obtenemos el sistema lineal sobredeter-
minado ⎛ ⎞ ⎛ ⎞
1 0 0 ⎛ ⎞ 1,1
⎜ 1 1 1/2 ⎟ x0 ⎜ 0,3 ⎟
⎜ ⎟ ⎜ ⎟
⎜ 1 2 2 ⎟ ⎝ v0 ⎠ = ⎜ −6,5 ⎟
⎜ ⎟ ⎜ ⎟
⎝ 1 3 9/2 ⎠ a ⎝ −21 ⎠
1 4 8 −39
9.2 Problemas resueltos 205
⎛ ⎞T ⎛ ⎞ ⎛ ⎞T ⎛ ⎞
1 0 0 1 0 0 ⎛ ⎞ 1 0 0 1,1
⎜ 1 1 1/2 ⎟ ⎜ 1 1 1/2 ⎟ x ⎜ 1 1 1/2 ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟ 0 ⎜ ⎟ ⎜ 0,3 ⎟
⎜ 1 2 2 ⎟ ⎜ 1 2 2 ⎟ ⎝ v0 ⎠ = ⎜ 1 2 2 ⎟ ⎜ −6,5 ⎟.
⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ 1 3 9/2 ⎠ ⎝ 1 3 9/2 ⎠ a ⎝ 1 3 9/2 ⎠ ⎝ −21 ⎠
1 4 8 1 4 8 1 4 8 −39
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
5 10 15 x0 −65,1
⎝ 10 30 50 ⎠ ⎝ v0 ⎠ = ⎝ −231,7 ⎠ ,
15 50 177/2 a −419,35
cuya solución es
t
1 2 3 4
-10
-20
-30
-40
Problema 9.5
La siguiente tabla muestra el volumen v de un gas en un cilindro cuando la
presión p en un pistón es aumentada
pk (Pa) 5 6 7
vk (m3 ) 2.5 1 0.5
(i) Aproximar, con todos los datos de la tabla adjunta, la función v(p) mediante
una interpolación polinómica de Lagrange.
(ii) Posteriormente se averigua que la evolución termodinámica del gas es
adiabática y por lo tanto verifica la ecuación pv γ = c con c y γ con-
stantes. Calcular dichas constantes por el método de mı́nimos cuadrados
a partir de las mediciones experimentales dadas.
(iii) Supongamos que se realiza una tabulación del volumen v del gas en los
n + 1 nudos dados por las presiones 3 = p0 < p1 < · · · < pn−1 < pn =
5 y posteriormente se calcula el polinomio de interpolación de Lagrange
Pn (p). Se sabe que se tienen N = 15 moles de gas a temperatura constante
T = 300 K y que además suponemos válida la aproximación de gas ideal,
es decir, se verifica la ecuación de estado pv = N RT con R = 8,314
J/K mol. Acotar el error absoluto Δn que se comete cuando se realiza la
aproximación v(p) ≈ Pn (p) para cualquier valor p ∈ [3, 5] en el caso de
que n = 3. Calcular (únicamente en función de n) la anterior acotación
de Δn para un valor arbitrario de n y obtener lı́mn→∞ Δn .
(iv) Un modelo de gas más realista donde se tienen en cuenta las interacciones
entre sus moléculas ası́ como su taman̄o mediante las constantes a y b
caracterı́sticas del gas viene dado por la ecuación de estado de van der
Waals
N 2a
p+ 2 (v − N b) = N RT .
v
Si el gas que se estudia es el Oxı́geno, se tiene a = 137,4, b = 0,031.
Suponiendo que T = 300, N = 15 y p = 2000
(iv.a) Reescribir la ecuación de estado de van der Waals en la forma
f (v) = 0, siendo f un polinomio cúbico.
(iv.b) Demostrar que el volumen v ∗ ocupado por el Oxı́geno está en el
intervalo [17, 20]. Partiendo de dicho intervalo, calcular el número
de iteraciones necesarias en el algoritmo de bisección para aproximar
v ∗ con un error absoluto menor que 10−5 .
(iv.c) Aproximar el valor de valor de v ∗ mediante 2 iteraciones con el
método de la secante partiendo de los valores iniciales v0 = 17 y
v1 = 20.
En definitiva se tiene
p − p1 p − p2 p−6 p−7 1
0 (p) = = = (p − 6)(p − 7) ,
p0 − p1 p0 − p2 5−6 5−7 2
p − p0 p − p2 p−5 p−7
1 (p) = = = −(p − 5)(p − 7) ,
p1 − p0 p1 − p2 6−5 6−7
p − p0 p − p1 p−5 p−6 1
2 (p) = = = (p − 5)(p − 6) ,
p2 − p0 p2 − p1 7−5 7−6 2
de manera que
1 1
P2 (p) = 2,5 × (p − 6)(p − 7) − (p − 5)(p − 7) + 0,5 × (p − 5)(p − 6)
2 2
1
= (50 − 14p + p2 ) .
2
(ii) En primer lugar observamos que el parámetro γ no interviene de forma
lineal en la ecuación pv γ = c. Por consiguiente utilizaremos, en primer lugar,
una estrategia de linearización de los datos. Para conseguirlo tomamos logar-
itmos neperianos en la anterior expresión y obtenemos ln c − γ ln v = ln p, de
manera que los parámetros ln c y γ intervienen de forma lineal. En segundo
lugar calculamos la nueva tabulación
Veamos como acotar esta expresión. En primer lugar, puesto que ξp ∈ (3, 5),
es evidente que 1/ξpn+2 < 1/3n+2 para cualquier n ∈ N. En segundo lugar,
como p ∈ [3, 5] y pk ∈ [3, 5] para cualquier k = 0, 1, . . . , n, se deduce que
n
|p − pk | ≤ 5 − 3 = 2 para k = 0, 1, . . . , n y por lo tanto k=0 |p − pk | ≤ 2n+1 .
En definitiva, hemos probado la siguiente cadena de desigualdades
37413 37413
n n
37413
Δn = |p − p k | < |p − pk | ≤ n+2 2n+1 .
ξpn+2 k=0 3n+2
k=0
3
n vn
0 17
1 20
2 18.3111
3 18.3516
La estabilidad del punto de equilibrio (θ, θ̇) = (0, 0) depende de dos con-
tribuciones: intuitivamente, el sistema sugiere que si el muelle es blando
(k pequen̄a) y la masa m es grande entonces el punto de equilibrio es es-
table y viceversa, es decir, k grande y m pequen̄a convertirá el punto de
equilibrio en inestable.
Si = k = 1, m = 100 y la masa se encuentra inicialmente en reposo
y formando un ángulo θ(0) = π/4, usar un paso del método RK4 con
una longitud de paso h = 1/5 para verificar las hipótesis de estabilidad
establecidas.
(ii) Suponiendo que el muelle verifica la ley de Hooke, geometrı́a elemental
muestra que la fuerza F que ejerce el muelle es
θ
F = k 2 cos − 1 .
2
θ̇ = f (t, θ, w) = w ,
1 θ θ
ẇ = g(t, θ, w) = −9,8 sin θ + 2 cos − 1 sin ,
100 2 2
con la condición inicial (θ(0), w(0)) = (π/4, 0). Puesto que f = f (w) sólo de-
pende de w y g = g(θ) sólo de θ, el método RK4 reduce a
θi+1 θi h r1 + 2r2 + 2r3 + r4
= + ,
wi+1 wi 6 k1 + 2k2 + 2k3 + k4
siendo
r1 = f (wi ) , k1 = g(θi ) ,
h h
r2 = f wi + k1 , k2 = g θi + r1 ,
2 2
h h
r3 = f wi + k2 , k3 = g θi + r2 ,
2 2
r4 = f (wi + hk3 ) , k4 = g (θi + hr3 ) .
Imponiendo la tabulación
es decir,
3,469k = 10,9108 ,
con lo que k = 3,14523.
donde <, > denota el producto escalar euclideano continuo en el intervalo [−π, π].
Se tiene pues que
! π
< ϕ 1 , ϕ1 > = dx = 2π ,
−π
! π
< ϕ2 , ϕ2 > = cos2 (x)dx = π
−π
! π
< ϕ3 , ϕ3 > = sin2 (x)dx = π
−π
! π
< ϕ 1 , ϕ2 > = cos(x)dx = 0
−π
! π
< ϕ1 , ϕ3 > = sin(x)dx = 0
−π
! π
< ϕ2 , ϕ3 > = sin(x) cos(x)dx = 0
−π
! π
< f, ϕ1 > = 1/2 xdx = 0
−π
! π
< f, ϕ2 > = 1/2 x cos(x)dx = 0
−π
! π
< f, ϕ3 > = 1/2 x sin(x)dx = π .
−π
Las integrales que aparecen en los anteriores productos escalares son, o bien in-
mediatas, o bien (segunda y tercera) se calculan con las identidades trigonométri-
cas
1 − cos(2x) 1 + cos(2x)
sin2 (x) = , cos2 (x) = ,
2 2
excepto la última que se calcula mediante la fórmula de integración por partes.
De hecho hay muchas de esas integrales (desde la quinta hasta la octava inclu-
idas) que valen cero puesto que la función integrando es impar y el intervalo de
integración es simétrico respecto del origen.
t 0 3 5 7
c 1 20 22 23
1
c(t) = .
α + β exp(−3t)
1
= α + β exp(−3t) .
c(t)
cuya solución es
α = 0,0462715 , β = 0,953729 .
En definitiva, la mejor aproximación viene dada por
1
c(t) = .
0,0462715 + 0,953729 exp(−3t)
Introducción a la
Programación
Do[
Print["Para T = ", T, "la velocidad del sonido es ", v[T]],
{T, 20, 35, 1}
]
217
218 Introducción a la Programación
PROGRAM VELOCIDAD
REAL T,DT,V
T=20.
DT=1.
10 V=331.*SQRT(1.+T/273.)
WRITE(*,*) ’ Para T= ’,T,’ La velocidad del sonido es ’, V
T=T+DT
IF(T .LT. 36.0) THEN
GO TO 10
ELSE
WRITE(*,*) ’Cálculo terminado’
STOP
END IF
END
ax2 + bx + c = 0 .
Sabemos que la naturaleza de las raı́ces depende del valor del discriminante
δ = b2 − 4ac. En concreto tenemos
1 √
(i) si δ > 0 =⇒ dos raı́ces reales y diferentes x± = −b ± δ .
2a
b
(ii) si δ = 0 =⇒ una única raı́z real doble x+ = x− = − .
2a
1 √
(iii) si δ < 0 =⇒ dos raı́ces complejas conjugadas x± = −b ± i −δ .
2a
El objetivo de esta práctica es realizar un programa para calcular la raı́z real
(si existe) de menor valor absoluto de una ecuación de segundo grado.
Un posible código Fortran para calcular la raı́z real de menor valor abso-
luto de una ecuación de segundo grado puede ser el dado en la Figura 10.2. El
análogo código Mathematica es, por ejemplo, el siguiente.
10.2 Ecuación de Segundo Grado 219
a = 1; b = 8; c = 5;
PROGRAM CUADRATICO
REAL A,B,C,DISCR,DENOM,T1,T2,XMIN
WRITE(*,*) ’ Introduce los coeficientes separados por comas ’
READ(*,*) A,B,C
DISCR=B*B-4.*A*C
IF(DISCR .LT. 0) THEN
WRITE(*,*) ’Raı́ces Complejas’
STOP
ELSE
DENOM=2.*A
T1=-B/DENOM
T2=SQRT(DISCR)/DENOM
IF(B .GE. 0) THEN
XMIN=T1+T2
ELSE
XMIN=T1-T2
END IF
END IF
WRITE(*,*)’ La raı́z real de menor valor absoluto es: ’,XMIN
STOP
END
Figura 10.2: Raı́z real de menor valor absoluto de una ecuación de segundo
grado.
220 Introducción a la Programación
mente a Heron de Alejandrı́a, quien vivió en el siglo I, pero probablemente es anterior. Natu-
ralmente, no fue obtenido a partir de lo que hoy se conoce como el método de Newton-Raphson
(descrito en secciones posteriores), sino por un razonamiento que muestra bien √ la estructura
del pensamiento matemático griego: si xn es una aproximación por exceso de a, entonces
a/xn es una aproximación por defecto y recı́procamente. Tomando la media aritmética entre
una aproximación por defecto y una por exceso se obtiene una aproximación que es mejor que
las anteriores.
2 Recordemos que el método de reducción al absurdo consiste en suponer cierto lo contrario
de lo que se quiere demostrar y comprobar que esta suposición implica algún tipo de absurdo.
10.3 Método de Newton: Cálculo de Raı́ces Cuadradas 221
a = 2;
x = (a + 1)/2;
PROGRAM NEWTON
REAL A,X
INTEGER I
WRITE(*,*) ’ Introduce el número cuya raı́z queremos hallar ’
READ(*,*) A
I=0
X=(A+1.)/2.
1 X=0.5*(X+A/X)
I=I+1
IF(I. GT. 100) THEN
WRITE(*,*) ’ No hay convergencia en 100 iteraciones ’
STOP
ELSE IF (ABS(X*X-A) .LT. 1.E-6) THEN
WRITE(*,*) ’ La raı́z cuadrada de ’, A, ’ es ’ , X
STOP
ELSE
GO TO 1
END IF
END
1 1 1
pk := Pk = 2k k = 2k 2 sin θk = 2k sin θk , (10.3)
2 2 2
10.4 Algoritmo de Arquimedes y el Número π 223
siendo el ángulo θk = 2πk , ver las Figuras 10.4 y 10.5. De la ecuación (10.3)
observamos que si doblamos el número de lados, es decir k pasa a k + 1, se tiene
θk
pk+1 = 2k+1 sin . (10.4)
2
Realizando ahora un poco de trigonometrı́a
θk 1 1
sin2 = (1 − cos θk ) = (1 − 1 − sin2 θk ) ,
2 2 2
k = 2;
P = 2*Sqrt[2];
PROGRAM PI
REAL P
INTEGER K
K=2
P=2*SQRT(2)
DO WHILE(K .LE. 15)
WRITE(*,*) K,P
P=SQRT(2)*2**K*SQRT(1-SQRT(1-(P/2**K)**2))
K=K+1
END DO
END
k pk
2 2.828427
3 3.061467
4 3.121445
5 3.136549
6 3.140331
7 3.141277
8 3.141514
9 3.141573
10 3.141588
11 3.141592
12 3.141593
13 3.141593
14 3.141593
15 3.141593
10.5. El problema 3x + 1
Vamos a describir un problema no resuelto en matemáticas conocido como el
problema 3x + 1. A partir de un número natural x ∈ N cualquiera, se construye
una sucesión de números naturales siguiendo el siguiente procedimiento:
function s=sucesion(x)
s=x;
while(x>1)
if(rem(x,2)==0)
x=x/2;
else
x=3*x+1;
endif
s=[s, x];
endwhile
endfunction
Capı́tulo 11
11.1. Introducción
Esta sección está dedicada a la presentación de una práctica en concreto con
el objetivo de que el alumno disponga de un material preliminar.
227
228 El Periodo del Péndulo Simple
1 2 g
θ̇ + (1 − cos θ) = C ,
2
siendo C una constante que se determina con las condiciones iniciales del movimien-
to del péndulo. Tomando θ̇ = 0 para tiempo t = 0 en el cual el péndulo forma
un ángulo θ = θ0 se tiene, despejando de la anterior ecuación
dθ 2g
= θ̇ = cos θ − cos θ0 .
dt
Recordando que el periodo T es cuatro veces el tiempo que tarda el péndulo en
ir de θ = 0 hasta θ = θ0 , separando variables e integrando se tiene
* !
θ0
dθ
T =4 √ .
2g 0 cos θ − cos θ0
sin(θ/2)
sin ξ = ,
sin(θ0 /2)
siendo ! π/2
dξ
K(k) := , (11.1)
0 1 − k 2 sin2 ξ
la integral elı́ptica completa de primera especie.
11.2 Algoritmo de Integración Gaussiana con 5 Nodos 229
PROGRAM INTGAUSS5
REAL L, K
COMMON THETA0
DIMENSION X(0:1), ALPHA(0:2)
PARAMETER(PI=3.141592, G=9.806)
X(0)= -0.906179
X(1)= -0.538469
ALPHA(0)= 0.236926
ALPHA(1)= 0.478628
ALPHA(2)= 0.568888
A=0.
B=PI*0.5
K=ALPHA(2)*F((A+B)*0.5)
11.3 Resultados Numéricos 231
DO I=0,1
U=((B-A)*X(I)+A+B)*0.5
V=(-(B-A)*X(I)+A+B)*0.5
K=K+ALPHA(I)*(F(U)+F(V))
END DO
K=0.5*(B-A)*K
T=4.*SQRT(L/G)*K
STOP
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
FUNCTION F(XI)
COMMON THETA0
F=1./SQRT(1.-(SIN(THETA0*0.5)*SIN(XI))**2)
RETURN
END
que
! π/2 ! π
2 −
dξ dξ
K(1) = = lı́m
2 cos ξ
0 1 − sin ξ →0+ 0
" # π2 −
= lı́m ln(sec ξ + tan ξ) =∞.
→0+ 0
Sistemas Lineales de
Ecuaciones: Eliminación
Gaussiana
Una vez hemos triangularizado la matriz de los coeficientes del sistema con el
233
234 Sistemas Lineales de Ecuaciones: Eliminación Gaussiana
Se puede realizar un único programa que realize los dos pasos (eliminación
Gaussiana y sustitución hacia atrás) o bién partirlo en subroutinas.
para i=n, 1, -1
xi = bi /uii
para j=i+1, n
xi = xi − uij xj /uii
finpara
finpara
PROGRAM GAUSS
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 M
PARAMETER(NMAX=100)
DIMENSION A(NMAX,NMAX), M(NMAX,NMAX), B(NMAX), X(NMAX)
C-----TRIANGULARIZACION
50 DO K=1,N-1
IF(A(K,K) .EQ. 0.) THEN
WRITE(*,*)’Hay un pivote nulo: MATRIZ SINGULAR’
STOP
END IF
END DO
END DO
C-----OPCIONES DE SALIDA
60 WRITE(*,*)’ Solucion por (1) PANTALLA, (2) FICHERO’
READ(*,*) L
IF ((L .NE. 1) .AND. (L .NE. 2)) GO TO 60
IF(L .EQ. 1) GO TO 70
C-----ADVERTENCIAS
100 WRITE(*,*)’Error de lectura en el fichero de datos’
STOP
200 WRITE(*,*)’Error de lectura en el fichero de datos’
STOP
END
n = Dimensions[A][[1]];
M = IdentityMatrix[n];
L = M; (* Matrices L y U *)
For[i = 1, i <= n, i++,
For[j = 1, j <= n, j++,
If[i >j, U[[i, j]] = 0, U[[i, j]] = A[[i, j]] ]
]
]
Print["L=",MatrixForm[L]] (*salida*)
Print[‘‘U=", MatrixForm[U]]
Print["x=", MatrixForm[x]]
Print["FIN"]
MatrixForm[AA]
MatrixForm[L]
MatrixForm[U]
MatrixForm[AA - L.U]
x
LinearSolve[AA, bb]
LinearSolve[U, b]
13.1. Introducción
Consideremos un material de coeficiente de conductividad térmica κ, calor
especı́fico c y densidad ρ tal que su anchura es despreciable respecto de su
longitud L (por ejemplo un alambre). Es fácil demostrar que, si el material
está aislado lateralmente y se verifica la ley de Fourier, entonces la ecuación que
govierna la evolución de la temperatura T (x, t) en un punto de posición x para
tiempo t es la llamada ecuación del calor
∂T ∂2T
(x, t) = σ 2 2 (x, t) , (13.1)
∂t ∂x
donde σ 2 = cρ κ
es el llamado coeficiente de difusividad. Este es un ejemplo
de ecuación en derivadas parciales parabólica. Para hallar una única solución
T (x, t) se necesita una condición inicial, la distribución inicial de temperatura
en el material
T (x, 0) = f (x) , 0 ≤ x ≤ L , (13.2)
y una condición de contorno, temperatura en los extremos del material
T (0, t) = T0 , T (L, t) = TL , 0 ≤ t < ∞ . (13.3)
Las constantes T0 y TL ası́ como la función f (x) son datos del problema.
241
242 La Ecuación del Calor Unidimensional
(xi , tj ) ∈ {(x, t) : 0 ≤ x ≤ L , 0 ≤ t ≤ τ } .
∂T T (x, t + k) − T (x, t)
(x, t + k/2) = + O(k 2 ) , (13.4)
∂t k
∂2T 1
(x, t + k/2) = [T (x − h, t + k) − 2T (x, t + k) + T (x + h, t + k)
∂x2 2h2
+T (x − h, t) − 2T (x, t) + T (x + h, t)] + O(h2 ) . (13.5)
h2
k= . (13.7)
σ2
Conocidos los T j se calculan a partir del sistema (13.6) los T j+1 . Además, hay
que tener en cuenta las condiciones de contorno T0j = T0j+1 = T0 y Tnj = Tnj+1 =
TL .
13.3 Separación de Variables: Series de Fourier 243
Impongamos que T (x, t) = ξ(x)η(t) sea una solución de la ecuación del calor
(13.1). Entonces se obtiene que
ξ(x)η̇(t) = σ 2 ξ (x)η(t) .
0 = η̇(t) + σ 2 λ2 η(t) ,
0 = ξ (x) + λ2 ξ(x) ,
Por lo tanto, según (13.15), la solución exacta T (x, t) viene dada por la serie
∞
nπ nπ 2 2 '
T (x, t) = αn sin x exp [−( σ t . (13.13)
n=1
L L
13.3 Separación de Variables: Series de Fourier 245
Sólo falta hallar los coeficientes αn , que vendrán determinados por la distribu-
ción inicial de temperaturas. En concreto, puesto que
∞
nπ
T (x, 0) = f (x) = αn sin x ,
n=1
L
FUNCTION F(X)
F=SIN(X)
RETURN
END
(iii) Puesto que todos los sistemas lineales (13.8) tienen la misma matriz de
coeficientes y lo único que varı́a son los términos independientes, un buen
método para su resolución consiste en hallar la descomposición LU de
la matriz de coeficientes del sistema (13.8). Además para hallar tal des-
composición LU , la eliminación Gaussiana será numéricamente estable.
Aprovechar además la estructura tridiagonal de la matriz de coeficientes
de (13.8).
(iv) Definir en una sentencia PARAMETER los valores máximos de n y m
asignándoles, por ejemplo, los valores N M AX = 100, M M AX = 100.
Dimensionar posteriormente las matrices L, U ∈ Mn−1 (R). Los elementos
T (xi , tj ) ≈ Tij almacenarlos también en una matriz que se habrá dimen-
sionado adecuadamente y que en el programa serán utilizados de la forma
T (I, J).
(v) Nota: Llamar en el programa LL a la longitud de la barra para no tener
confusiones con la matriz triangular inferior L.
(vi) La salida del programa será un fichero, por ejemplo temperaturas.dat, que
contenga las parejas de valores (xi , Tij ) ordenados de la forma
x0 , T0
246 La Ecuación del Calor Unidimensional
x1 , T10
..
.
0
xn−1 , Tn−1
xn , TL
x0 , T0
x1 , T11
..
.
1
xn−1 , Tn−1
xn , TL
..
.
x0 , T0
x1 , T1m
..
.
m
xn−1 , Tn−1
xn , TL
3. Extraer de la lista temp las parejas de puntos (xi , Tij ) para cada val-
or de j, es decir, las posiciones y temperaturas de la barra par tiempo tj con
j = 0, 1, . . . , m e incluirlos en sublistas llamadas t[j]:
s nπ nπ 2 2 '
T (x, t) ≈ αn sin x exp [−( σ t , (13.15)
n=1
L L
f1[x ] = x
f2[x ] = 1-x
lambda[n ]=n*Pi/L
Taprox[x ,t ]=Sum[T[x,t,n],{n,1,s}]
Interpolación de Newton
14.1. Introducción
El problema de la interpolación tiene una relevancia especial en muchas
aplicaciones a la ingenierı́a. Los métodos de interpolación tratan de encontrar
una función (polinómica, trigonométrica, exponencial,...) que satisfaga ciertas
condiciones en un conjunto finito de puntos. Estas condiciones pueden darse
sobre los valores de la función, y también sobre sus derivadas. En diferentes
disciplinas cientı́ficas y técnicas se plantea muy a menudo este problema. De
hecho, son muchas y muy distintas las situaciones en las que aparecen series
de datos o resultados de mediciones experimentales de los cuales sólo se conoce
cómo se comportan en una cierta cantidad finita de puntos y para los cuales
se necesita encontrar una “ley general”. Dicha ley no es mas que una función
que tome unos valores predeterminados. Éste es, precisamente, el cometido de
la interpolación. La función interpoladora servirá para sustituir a la función
desconocida, tanto para evaluarla en puntos en los que no se conoce su val-
or (interpolación en sentido estricto), como para conocer su tasa de variación
(diferenciación numérica) o su distribución acumulativa (integración numérica).
Además, la interpolación sirve para fundamentar una amplia gama de métodos
para el tratamiento numérico de ecuaciones diferenciales.
La función interpoladora debe ser, por tanto, fácil de evaluar, derivar e in-
tegrar.
Nos limitaremos a la interpolación polinómica que consiste en lo siguiente:
249
250 Interpolación de Newton
donde
i−1
ϕ0 (x) := 1 , ϕi (x) = (x − xk ) para i = 1, 2, . . . , n.
k=0
Los coeficientes ci se conocen con el nombre de diferencias divididas y su
notación histórica es
ci := f [x0 , x1 , . . . , xi ] .
Se puede demostrar que las diferencias divididas vienen definidas de la forma
siguiente
f [xi ] = fi ,
f [xi+1 , xi+2 , . . . , xi+j ] − f [xi , xi+1 , . . . , xi+j−1 ]
f [xi , xi+1 , . . . , xi+j ] = .
xi+j − xi
Un ejemplo de como esquematizar el proceso descrito para el cálculo de las
diferencias divididas en el caso n = 3 es el siguiente
x0 f [x0 ]
f [x0 , x1 ]
x1 f [x1 ] f [x0 , x1 , x2 ]
f [x1 , x2 ] f [x0 , x1 , x2 , x3 ]
x2 f [x2 ] f [x1 , x2 , x3 ]
f [x2 , x3 ]
x3 f [x3 ]
A continuación, ver Figura 14.1, damos un algoritmo para el cálculo de las
diferencias divididas.
Este algoritmo permite el cálculo de los coeficientes ci que no son más que
f [x0 , x1 , . . . , xi ]. Una forma de ahorrar memoria en la computadora consiste en
dimensionar un vector c = (c0 , c1 , . . . , cn ). Al introducir los datos de entrada
f0 , f1 , . . . , fn , se almacenan en el vector c, de manera que c0 ya es el primer
coeficiente del polinomio interpolador. Después se computa la primera columna
del esquema de diferencias divididas y colocamos estos valores en c1 , c2 , . . . , cn ,
de manera que c0 conserva su valor y el nuevo c1 es el segundo coeficiente del
polinomio interpolador y ası́ sucesivamente.
14.3 Fenómeno Runge 251
para i=0,n
ci = fi
finpara
para j=1,n
para i=n,j
ci = (ci − ci−1 )/(xi − xi−j )
finpara
finpara
n
Pn (x) = ci ϕi (x)
i=0
= c0 + c1 (x − x0 ) + · · · + cn (x − x0 ) · · · (x − xn−1 ) .
Como resultado final de esta práctica se tiene que apreciar los efectos de
borde: el polinomio interpolador tiende a oscilar en los extremos del intervalo de
interpolación [−5, 5].
1 En esta práctica tomamos esta forma del polinomio aunque recordemos que, la forma más
conveniente desde el punto de vista del número de operaciones necesarias para de evaluar
un polinomio es la de Horner. Por ejemplo el polinomio de cuarto grado P4 (x) = 4i=0 ai xi
adopta la forma de Horner P4 (x) = a0 + x(a1 + x(a2 + x(a3 + aa4 ))).
252 Interpolación de Newton
inicializar x
valor=c0
p=1
para i=0,n-1
p = p(x − xi )
valor=valor+ci p
finpara
salida x, valor
siendo
n
Πn (x) = (x − xi ) .
i=0
Teorema 14.1 (Chebychev) Dada una función f ∈ C n+1 [−1, 1], el menor
error en la norma del máximo en la interpolación de f se comete cuando se
toman las siguientes abscisas de interpolación {x0 , x1 , . . . , xn } ⊂ [−1, 1]
(2j + 1)π
xj = cos , j = 0, 1, . . . , n .
2(n + 1)
Notas:
15.1. Introducción
El fenómeno Runge mostrado en la práctica anterior nos está indicando que
la interpolación polinomial sobre un número grande de abscisas no es siempre
adecuada debido a la tendencia a oscilar que tienen los polinomos. Una solu-
ción adecuada para este problema es la interpolación mediante funciones spline1 .
sirve para trazar curvas suaves. Precisamente, la propiedad de adaptarse bien a formas dadas
que tienen las funciones spline es lo que hace que se les dé tal nombre
255
256 Interpolación por Splines Cúbicos Naturales
1 zi hi hi 1
Ai := (zi+1 −zi ) , Bi := , Ci := − zi+1 − zi + (fi+1 −fi ) , (15.3)
6hi 2 6 3 hi
siendo
hi := xi+1 − xi , zi := S (xi ) . (15.4)
Observar que lo único que nos queda por hacer para finalizar el problema de
la interpolación de splines es calcular los valores z0 , z1 , . . . , zn .
x − xn−1 , x − xn−2 , . . . , x − x1 ,
Un código para resolver el sistema lineal (15.5) puede ser el dado en la Figura
15.1.
15.2 Realización de la Práctica 257
entrada n, xi , fi
th1 = 2(h0 + h1 )
tf1 = b1 − b0
para i=2, n-1
thi = 2(hi + hi−1 ) − h2i−1 /thi−1
tfi = bi − bi−1 − hi−1 tfi−1 /thi−1
finpara
zn = 0
para i=n-1, 1, -1
zi = (tfi − hi zi+1 )/thi
finpara
z0 = 0
salida zi
Figura 15.1: Algoritmo para generar las segundas derivadas de splines cúbicos
naturales calculadas en los nodos.
Mostrar
√ (con un software adecuado) las gráficas superpuestas de la función
x y de la función spline S(x) obtenida en el apartado anterior. Por
ejemplo, con el software Mathematica:
spline=ListPlot[datos, PlotJoined→True]
√
(iii) Dibujar una gráfica de la función x en el intervalo [0, 1]
PROGRAM SPLINE
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER(NMAX=100)
DIMENSION X(0:NMAX), F(0:NMAX), Z(0:NMAX)
DIMENSION FT(0:NMAX), H(0:NMAX), HT(0:NMAX), B(0:NMAX)
C-----INICIALIZACION DE PARAMETROS
15.3 Solución en Fortran 259
50 DO I=0,N-1
H(I)=X(I+1)-X(I)
B(I)=6*(F(I+1)-F(I))/H(I)
END DO
DO I=1,N-1
HT(I)=2*(H(I)+H(I-1))
FT(I)=B(I)-B(I-1)
END DO
C-----INICIALIZACION DE PARAMETROS
DO I=2,N-1
HT(I)=HT(I)-H(I-1)**2/HT(I-1)
FT(I)=FT(I)-H(I-1)*FT(I-1)/HT(I-1)
END DO
C-----SPLINES NATURALES
Z(0)=0.
Z(N)=0.
C-----CALCULO DE LA APROXIMACION
55 AA=(Z(I+1)-Z(I))/(6.*H(I))
BB=Z(I)/2.
CC=-H(I)*Z(I+1)/6.-H(I)*Z(I)/3.+(F(I+1)-F(I))/H(I)
S=F(I)+(XX-X(I))*(CC+(XX-X(I))*(BB+(XX-X(I))*AA))
WRITE(*,*)’ La aproximacion calculada vale’, S
STOP
C-----ADVERTENCIAS
260 Interpolación por Splines Cúbicos Naturales
n = Dimensions[tabul][[1]];
z = LinearSolve[A, bb];
(* resuelvo sistema lineal e impongo splines naturales *)
z = AppendTo[z, 0];
z = PrependTo[z, 0];
grafspline = { };
For[i = 1, i <= n - 1, i++, AppendTo[grafspline, graf[i]]];
grafspline = Show[grafspline,
DisplayFunction ->Identity];
0.5
x
4 6 8 10
-0.5
Capı́tulo 16
Derivación Numérica
16.1. Introducción
Para la aproximación de la derivada de una función f (x) en un punto x0 ; los
paquetes de cálculo simbólico, por ejemplo el comando D[f[x],x] de Mathemat-
ica, proporcionan la expresión de la derivada f (x) y, por lo tanto, sólo queda
evaluar dicha función en el punto dado para obtener f (x0 ). Sin embargo, en
el caso en que tan sólo se conocen los valores de la función en algunos puntos,
aquı́ las fórmulas de derivación aproximada si serán de gran utilidad.
f (x) − f (x − h)
f (x) = + O(h) . (16.6)
h
Restando la ecuación (16.2) de la (16.3) se deduce la fórmula de diferencias
centradas
f (x + h) − f (x − h)
f (x) = + O(h2 ) . (16.7)
2h
Combinando (16.1), (16.2), (16.3) y (16.4) podemos obtener
−f (x + 2h) + 8f (x + h) − 8f (x − h) + f (x − 2h)
f (x) = + O(h4 ) . (16.8)
12h
Otra forma de obtener la fórmula (16.8) es mediante extrapolación de Richard-
son a partir de (16.7).
READ(*,*) CRND
Para modificar los valores de la función, utilizar el siguiente bucle
DO I=0,N
RND=RND*CRND
RND=DMOD(1.6807D4*RND, 2.147483647D9)
RAND=RND*4.65661287524579692D-10
F(I)=F(I)+0.01*RAND
END DO
De esta forma, la variable RAND toma valores aleatorios entre 0 y 1.
Capı́tulo 17
Integración Romberg
17.1. Introducción
Uno de los problemas matemáticos más antiguos es el del cálculo del área
que encierra una curva. Este problema desembocó en lo que hoy se conoce como
cálculo integral.
Las fórmulas de integración numéricas, también conocidas como fórmulas de
cuadratura, tienen como objetivo aproximar el valor de la integral de una fun-
b
ción en un intervalo a f (x) dx. El cálculo simbólico mediante un manipulador
algebraico es de gran ayuda, ver por ejemplo el comando Integrate[f[x], {x,a,b}]
del programa Mathematica. Sin embargo esta ayuda sólo es posible cuando
se conoce la expresión explı́cita de f (x) y además cuando dicha función admite
una primitiva. En el resto de los casos, por ejemplo cuando sólo se conoce los
valores de la función en algunos puntos o cuando la primitiva de la función no
es posible expresarla mediante funciones elementales, habrá que acudir a las
fórmulas de integración numérica.
267
268 Integración Romberg
! b
f (x) dx = R(2h, k − 1) + a1 (2h)2k + a2 (2h)2k+2 + · · · ,
a
T (1, 1)
T (2, 2)
T (2, 1) T (3, 3)
T (3, 2) T (4, 4)
T (3, 1) T (4, 3)
T (4, 2)
T (4, 1)
(i) Si se programa en Fortran, dimensionar una matriz 100 × 100 para alma-
cenar los valores de T (j, k). En Octave no es necesario.
(ii) Las entradas del programa tienen que ser
El intervalo de integración [a, b].
Dos criterios para finalizar las iteraciones: (1) El número máximo
permitido para j (por ejemplo 20). (2) La toleracia (por ejemplo
10−5 ). Pararemos las iteraciones cuando se verifique la condición
|T (j, j) − T (j − 1, j − 1)| < tolerancia .
17.2 Realización de la Práctica 269
(iii) La salida del programa debe ser la escritura en un fichero de las parejas
j, T (j, j) para j = 0, 1, 2, . . .
Pb Na Ag Cu Al Be
ΘD 86◦ K 160◦ K 220◦ K 310◦ K 380◦ K 980◦ K
Ecuaciones Integrales:
Transmisión de Calor
18.1. Introducción
Consideremos dos placas de longitud infinita, anchura y separadas una
distancia d. Supongamos que las placas están a temperaturas constantes T1 y
T2 , siendo sus emisividades 1 y 2 respectivamente. Tomemos dos sistemas de
referencia unidimensionales x1 y x2 con el origen en un punto cualquiera del
centro de las placas y dirigidos hacia sus extremos finitos respectivamente. Se
puede demostrar que las radiosidades locales B1 y B2 , es decir, la suma de las
cantidades de calor (radiación) emitida y reflejada que abandonan una unidad de
área por unidad de tiempo, verifican el siguiente sistema de ecuaciones integrales
! /2
4
B1 (x1 ) = 1 σT1 + (1 − 1 ) B2 (x2 )K(x1 , x2 , d) dx2 , (18.1)
−/2
! /2
B2 (x2 ) = 2 σT24 + (1 − 2 ) B1 (x1 )K(x1 , x2 , d) dx1 , (18.2)
−/2
271
272 Ecuaciones Integrales: Transmisión de Calor
debido a la radiación emitida por todos los puntos de la otra placa. De forma
totalmente análoga se tiene la irradiosidad en la otra placa
! /2
I2 (x2 ) = (1 − 2 ) B1 (x1 )K(x1 , x2 , d) dx1 . (18.5)
−/2
Los caudales netos de calor por unidad de longitud Q1 y Q2 que se deben propor-
cionar a dichas placas con el objetivo de mantenerlas a temperatura constante
son
! /2
Q1 = (B1 (x1 ) − I1 (x1 )) dx1 , (18.6)
−/2
! /2
Q2 = (B2 (x2 ) − I2 (x2 )) dx2 . (18.7)
−/2
tirá suponer que B1 (x1 ) e I1 (x1 ) son funciones simétricas respecto del centro de la placa
x1 = 0 y, de forma análoga, B2 (x2 ) e I2 (x2 ) son simétricas respecto de x2 = 0.
18.3 Existencia de Solución para Ecuaciones de Fredholm 273
que nos permite estimar según (18.1) y (18.2) las aproximaciones de orden
uno para B1 (x1 ) y B2 (x2 ). Estos valores son introducidos nuevamente
en las integrales (18.4) y (18.5) para posteriormente obtener mediante
(18.1) y (18.2) las aproximaciones de orden dos para B1 (x1 ) y B2 (x2 ) y
ası́ sucesivamente. El proceso iterativo terminará o bien cuando se llegue
a un número máximo ITMAX de iteraciones permitido por el usuario, o
bien cuando los valores de B1 y B2 no cambien más que un cierto valor
de toleracia TOL en dos iteraciones consecutivas.
(iv) Cada vez que se aproximen, según el algoritmo descrito anteriormente, los
valores de las integrales I1 (x1 ) e I2 (x2 ) dadas por (18.4) y (18.5) se re-
alizará una llamada a la subroutina o función SIMPSON(A,B,N,F). Nótese
que el número de aplicaciones de la regla de Simpson compuesta es n y
que a su vez es también el número de abscisas en cada semiancho de las
placas.
! b
yn (x) = f (x) + α K(x, t, yn−1 (t)) dt ,
a
genera una sucesión de funciones {yn (x)} que converge uniformemente hacia una
solución y(x) de (18.8) cuando n → ∞ para cualquier valor de α que verifique
1 f¯
|α| < , M := máx L, K̄ 1 + .
M (b − a) |α|K̄(b − a)
Aproximación
Trigonométrica y
Coeficientes de Fourier
19.1. Introducción
Se define el espacio vectorial real Tn como el generado por el conjunto de
funciones
1
, cos x, sin x, cos 2x, sin 2x, cos 3x, sin 3x, . . . , cos nx, sin nx , (19.1)
2
a0 n n
gn (x) = + am cos mx + bm sin mx , (19.2)
2 m=1 m=1
f − f ∗ = mı́n f − gn 2 ,
gn ∈Tn
275
276 Aproximación Trigonométrica y Coeficientes de Fourier
mı́nimos cuadrados.
! 2π
< sin x, sin mx > = sin x sin mx dx
0
! 2π
1
= [cos( − m)x − cos( + m)x] dx ,
2 0
! 2π
< cos x, cos mx > = cos x cos mx dx
0
! 2π
1
= [cos( − m)x + cos( + m)x] dx ,
2 0
! 2π
< sin x, cos mx > = sin x cos mx dx
0
! 2π
1
= [sin( − m)x + sin( + m)x] dx .
2 0
Por lo tanto
0 si = m ,
< sin x, sin mx >=< cos x, cos mx >=< sin x, cos mx >=
π si =m.
para todo , m ∈ N.
Sea
a∗0 n n
f ∗ (x) = + a∗m cos mx + b∗m sin mx , (19.3)
2 m=1 m=1
Distribución de Velocidades
en un Fluido
20.1. Introducción
Se ha deducido1 una expresión para la distribución de velocidades en un
fluido en movimiento por una tuberı́a circular lisa mediante una modificación de
la llamada longitud de mezcla de Prandtl. En concreto, la velocidad adimensional
axial u es una función de la distancia adimensional y desde la pared interior de
la tuberı́a dada por la expresión integral siguiente
! y
−1 + 1 + 4c(y)d(y)
u(y) = dy , (20.1)
0 2c(y)
siendo
" y
#2
c(y) = 0,36y 1 − e−φ ȳ , (20.2)
y
d(y) = 1− , (20.3)
ȳ
donde ȳ es la distancia adimensional del interior de la tuberı́a a su eje, cuya
expresión viene dada por
NRe f¯
ȳ = , (20.4)
2 2
y φ es una función de ȳ y de constantes empı́ricas dada por
ȳ − 60
φ= . (20.5)
22
En la expresión (20.4), NRe y f¯ representan el número de Reynolds y el factor
de fricción de Fanning respectivamente, siendo estos parámetros adimensionales
1 W.N. Gill y M. Schear, A modification of the momentum transport hypothesis,
279
280 Distribución de Velocidades en un Fluido
siendo
b−a b+a
f (z) = 1 − z2 g z+ .
2 2
20.3 Realización de la Práctica 281
La Condensación de
Bose-Einstein
21.1. Introducción
El premio Nobel de Fı́sica del an̄o 2001 fué otorgado a los cientı́ficos Carl
E. Wieman (Universidad de Colorado, EEUU), Eric A. Cornell (National Insti-
tute of Standards & Technology, Colorado, EEUU) y Wolfgang Ketterle (Mas-
sachusetts Institute of Technology, Cambridge, Massachusetts, EEUU). Los galar-
donados crearon, de forma independiente, en el an̄o 1995, por primera vez en el
laboratorio un condensado de Bose-Einstein. De hecho, las relaciones existentes
entre la condensación de Bose-Einstein y otros fenómenos de la fı́sica de bajas
temperaturas como es la superfluidez (descubierta en el 4 He lı́quido por Kapitza
en 1937) ha sido un tópico continuo en las investigaciones realizadas en los últi-
mos 60 an̄os.
283
284 La Condensación de Bose-Einstein
1
n = " # , (21.1)
−μ
exp KB T −1
1
n = " # . (21.2)
z −1 exp
KB T −1
1
N = N0 + " # , (21.3)
>0 z −1 exp
KB T −1
1 Ver por ejemplo K. Huang, Statistical Mechanics, 2a ed. Wiley, New-York, 1987.
21.2 Un Gas de Bosones Ideal 285
Puede ocurrir que la temperatura T sea tan grande que N0 = 0. Pero, por
el contrario, si la temperatura es suficientemente baja entonces N0 puede ser
muy grande. En el caso extremo en que z → 1 se puede definir una temperatura
crı́tica Tc definida por
2/3
2π2 N
Tc = (21.7)
m B3/2 (1)V
de manera que, si T > Tc entonces no hay una ocupación del estado fundamental
que sea macroscópicamente apreciable. Por el contrario, si T < Tc , un número
macroscópico de partı́culas pueden ocupar el estado fundamental. Obsérverse
que el fenómeno que se acaba de describir implica que el gas puede ser conside-
rado como una mezcla de dos fases; una fase constituida por aquellas partı́culas
que se encuentran en el estado fundamental y la otra fase por el resto de las
partı́culas. Una forma interesante de ver el fenómeno consiste en reescribir la
ecuación (21.6) de la forma
3/2
N0 T
=1− , (21.8)
N Tc
con lo que la Tc puede interpretarse como la temperatura de transición entre
las dos fases del gas: para T << Tc prácticamente todas las partı́culas del gas
están en el estado fundamental y se ha formado la llamada condensación de
Bose-Einstein3 .
2 Ver el apéndice.
3 Nótese las semejanzas con la clásica condensación gas-lı́quido. En realidad, la diferencia
básica entre el cambio de fase de la clásica condensación gas-lı́quido y la condensación de
Bose-Einstein reside en que la primera es debida a la interacción entre las partı́culas mientras
que la segunda es un fenómeno puramente cuántico y se puede dar en un gas ideal sin ningún
tipo de interacción entre partı́culas.
286 La Condensación de Bose-Einstein
Las integrales de Bose (21.9) han sido estudiadas por J.E. Robinson, Phys.
Rev. 83, pag. 678 (1951).
Capı́tulo 22
Ecuaciones Diferenciales
Ordinarias: Método RK4
22.1. Introducción
Una de las herramientas matemáticas más efectivas de que disponemos ac-
tualmente para modelizar y explicar la naturaleza son, sin lugar a dudas, las
ecuaciones diferenciales (o más generalmente los sistemas dinámicos). Una re-
solución exacta de las ecuaciones diferenciales no es posible en general. De esta
forma, son necesarias vias alternativas tales como los métodos numéricos o la
teorı́a cualitativa para poder estudiar las ecuaciones diferenciales. Por lo que
respecta a los métodos numéricos, el algoritmo que presentamos es uno de los
más populares para aproximar la solución del problema de valor inicial
ẋ = f (t, x, y) , x(t0 ) = x0 ,
(22.1)
ẏ = g(t, x, y) , y(t0 ) = y0 ,
t0 t1 t2 ··· tn
x0 x1 x2 ··· xn
y0 y1 y2 ··· yn
donde xi e yi sean una aproximación numérica de los valores exactos x(ti ) e y(ti )
respectivamente. Además, siempre consideraremos una discretización equiespa-
ciada de la variable independiente ti = t0 + ih para i = 0, 1, . . . , n, siendo h > 0
287
288 Ecuaciones Diferenciales Ordinarias: Método RK4
la longitud de paso.
siendo
r1 = f (ti , xi , yi ) ,
h h h
r2 = f ti + , xi + r1 , yi + k1 ,
2 2 2
h h h
r3 = f ti + , xi + r2 , yi + k2 ,
2 2 2
r4 = f (ti + h, xi + hr3 , yi + hk3 ) ,
k1 = g(ti , xi , yi ) ,
h h h
k2 = g ti + , xi + r1 , yi + k1 ,
2 2 2
h h h
k3 = g ti + , xi + r2 , yi + k2 ,
2 2 2
k4 = g (ti + h, xi + hr3 , yi + hk3 ) .
k1 = f (ti , xi ) ,
h h
k2 = f ti + , xi + k1 ,
2 2
22.3 Ciclos Lı́mite 289
h h
k3 = f ti + , xi + k2 ,
2 2
k4 = f (ti + h, xi + hk3 ) .
El número de puntos n.
El vector de condiciones iniciales x0 ∈ Rd .
La salida de la función debe ser:
d2 x dx
2
+ f (x) + g(x) = 0 . (22.5)
dt dt
Nótese que una partı́cula cuya ecuación del movimiento esté regida por la
ecuación (22.5) no puede tener soluciones periódicas cuyas órbitas estén en
la región en la que f (x) sea de signo constante puesto que fı́sicamente estas
zonas corresponden a regiones de amortiguamiento positivo o negativo. Solu-
ciones periódicas, es decir cuyas trayectorias formen curvas cerradas en el plano
de las fases, desempen̄an un papel importante en muchos fenómenos fı́sicos. Una
solución periódica a menudo representa una clase de “estado final”, hacia el cual
todas las soluciones vecinas tienden. La existencia de estas soluciones periódi-
cas (ciclos lı́mite) sólo pueden explicarse por la consideración de términos no
lineales en las ecuaciones que gobiernan el movimiento.
d2 x dx
2
+ μ(x2 − 1) +x=0 , (22.7)
dt dt
siendo μ un parámetro real positivo por razones fı́sicas.
Además, puede demostrarse que el ciclo lı́mite C es estable, es decir, todas las
órbitas excepto el punto de equilibrio (0, 0) tienden asintóticamente hacia el
ciclo lı́mite C cuando t → ∞, ver Figura 22.1.
de radio-ingenierı́a, inició el estudio de la ecuación que lleva su nombre en los an̄os veinte,
estimulando con ello a Liénard y otros a investigar la teorı́a matemática de las oscilaciones
automantenidas en mecánica no lineal. En su ecuación, la variable x denota la intensidad de
corriente eléctrica en un oscilador de triodo. Experimentalmente se observa que el flujo de
corriente es periódico. Este flujo periódico sólo puede ser explicado mediante la consideración
de una ecuación no lineal.
22.4 El Problema del Centro-Foco 291
dx dy
= −y , =x. (22.8)
dt dt
Este sistema es trivial de integrar puesto que dy/dx = −x/y. Entonces, sepa-
rando las variables e integrando, obtenemos que las órbitas en el plano de las
fases x − y vienen dadas por circunferencias concéntricas x2 + y 2 = C 2 . Como el
origen (0, 0) posee un entorno totalmente relleno de órbitas cerradas del sistema
diferencial, se dice que el punto crı́tico (0, 0) es un centro.
Se sabe que, las órbitas en el plano de las fases del sistema perturbado (22.9)
pueden seguir siendo periódicas, de manera que el punto crı́tico (0, 0) continua
siendo un centro, o bien espiralan alrededor del punto crı́tico en cuyo caso se
292 Ecuaciones Diferenciales Ordinarias: Método RK4
dice que el punto crı́tico (0, 0) es un foco, ver las Figuras 22.2 y 22.3. Averiguar
si el punto (0, 0) es un centro o un foco para el sistema (22.9) es un problema
clásico y difı́cil.
2 R.E. Berg and T. S. Marshall, Wilberforce pendulum oscillations and normal modes, Am.
arriba por conducción. La tendencia del fluido (menos denso) en la parte inferior
a elevarse se cancela debido a la pérdida de calor hacia el entorno y el rozamiento
debido a la viscosidad. En estas condiciones estacionarias, la temperatura del
fluido varı́a linealmente respecto de la posición vertical.
Por el contrario, si la diferencia de temperaturas ΔT aumenta suficiente-
mente, las fuerzas flotantes llegan a ser suficientemente grandes como para
vencer al rozamiento viscoso y se desarrollan corrientes estacionarias dentro
del fluido. En estas condiciones, la transferencia de calor desde la parte inferior
hacia la parte superior de la capa de fluido se realiza por convección, es decir,
debida al movimiento de la masa del fluido. Como resultado se obtiene una
circulación de fluido que es estable en el tiempo.
Con un aumento posterior de la diferencia de temperaturas ΔT , tanto las
corrientes de fluido como las diferencias de temperaturas dentro del fluido dejan
de ser estacionarias y comienzan a variar con el tiempo.
Si r > 1, pero no muy grande (por ejemplo menor que 13) se puede ver
que toda solución que comience cercana al punto (0, 0, 0) se aleja de este
y se acerca (con un tiempo suficiente) a alguno de los puntos
A± = (± b(r − 1), ± b(r − 1), r − 1)
Raı́ces de Ecuaciones:
Método de
Newton-Raphson
23.1. Introducción
f (xn )
xn+1 = xn − . (23.1)
f (xn )
299
300 Raı́ces de Ecuaciones: Método de Newton-Raphson
nRT
P = , (23.2)
V
siendo n el número de moles del gas, T su temperatura y R = 8,3143JK −1 mol−1
la constante de los gases ideales.
nRT n2 n3 n4
P = + A 2 + B 3 + C 4 + ··· . (23.3)
V V V V
La ecuación (23.3) se puede considerar como la ecuación de estado de un gas
real. A, B, C, ... son cantidades caracterı́sticas de cada gas y son llamados se-
gundo, tercer, cuarto, etc.. coeficientes viriales. Estos coeficientes son funciones
de la temperatura T y dependen de la intensidad de las fuerzas moleculares.
23.2 Ecuaciones de Estado de los Gases 301
Es fácil ver que un gas real que sigue la ecuación de van der Waals
n2 a
P + 2 (V − nb) = nRT , (23.5)
V
donde a y b son las constantes de van der Waals caracterı́sticas del gas, tiene
un módulo de elasticidad
V 2 (V − nb)2
κT = . (23.6)
nRT V 3 − 2n2 a(V − nb)2
Problema 23.4 Se sabe que el Oxı́geno tiene unas constantes de van der Waals
a = 137,4 y b = 0,031. Aproximar mediante el método el Newton-Raphson el
volumen que ocupan 30 moles de Oxı́geno a temperatura T = 310K sabiendo que
su módulo de elasticidad es κT = 2,74P a−1 y el volumen ocupado es superior a
200000m3 .
302 Raı́ces de Ecuaciones: Método de Newton-Raphson
Es fácil probar que la única solución ρ(x) del problema de valor final formado
por las ecuaciones (23.7) y (23.9) viene dada por
!
ρ(x) = − (s − x)a(s) ds . (23.10)
x
x = (x1 , x2 , . . . , xn )T .
d2 u du
a(x) + b(x) + c(x)u(x) = g(x) (24.1)
dx2 dx
definida para x perteneciente al domino [0, L] ⊂ R. Supondremos que las fun-
ciones a(x), b(x), c(x) y g(x) son continuas en dicho dominio. Si la función a(x)
tiene signo constante en el intervalo [0, L], es decir, (24.1) es una ecuación elı́pti-
ca, entonces, (24.1) puede reescribirse como una ecuación diferencial autoadjunta
de la forma
d du
− a1 (x) + a0 (x)u = f (x) , (24.2)
dx dx
siendo
! x
b(t)
a1 (x) := exp dt ,
0 a(t)
−c(x)a1 (x)
a0 (x) := ,
a(x)
−g(x)a1 (x)
f (x) := .
a(x)
307
308 Método de los Elementos Finitos
du
u(0) = Ū0 , a1 (L)(L) = QL , (24.4)
dx
du du
−a1 (0) (0) = Q0 , a1 (L) (L) = QL , (24.5)
dx dx
du
u(0) = Ū0 , a1 (L) (L) + βu(L) = γ , (24.6)
dx
n
Uk (x) := ukj φkj (x) . (24.8)
j=1
Definiendo
du du
QA := −a1 (xA ) (xA ) , QB := a1 (xB ) (xB ) , (24.10)
dx dx
310 Método de los Elementos Finitos
siendo las funciones aproximadoras φkj (x) los polinomios básicos de Lagrange
x − xk+1 x − xk
φk1 (x) = , φk2 (x) = .
xk − xk+1 xk+1 − xk
Observar que se verifica Uk (xk ) = uk1 y Uk (xk+1 ) = uk2 para todo k = 1, . . . , N .
# 2
φk1 (x)f (x) dx − φk1 (xkj )Qkj ,
j=1
⎛ ⎞ ⎛ ⎞
! xB " 2 2
dφk2 ⎝ dφkj ⎠
0 = a1 (x) ukj + a0 (x)φk2 ⎝ ukj φkj (x)⎠ −
xA dx j=1
dx j=1
# 2
φk2 (x)f (x) dx − φk2 (xkj )Qkj ,
j=1
para cada elemento Ωk . Aquı́ se han denotado por xkj , con j = 1, 2, los nodos
del elemento Ωk y por Qk1 = QA , Qk2 = QB .
24.2 Solución Numérica 311
Teniendo en cuenta que los polinomios basicos de Lagrange φki (x) verifican
la propiedad φki (xkj ) = δij , siendo δij la delta de Kronecker, es decir, δij = 1
si i = j y δij = 0 para i = 0, podemos escribir la ecuación i−ésima del método
de Galerkin mediante un reordenamiento como
2 ! xB 2 ! xB
dφki dφkj
0 = a1 (x) dx ukj + a0 (x)φki φkj dx ukj
j=1 xA dx dx j=1 xA
! xB
− a0 (x)f (x)φki dx − Qki , i = 1, 2 .
xA
donde la matriz de coeficientes K [k] , llamada matriz de rigidez, viene dada por
[k] [k] [k,1] [k,0]
K [k] = (Kij ) ∈ M2 (R) , Kij = Kij + Kij , (24.14)
! xB
[k,1] dφki dφkj
Kij = a1 (x) dx , (24.15)
xA dx dx
! xB
[k,0]
Kij = a0 (x)φki φkj dx , (24.16)
xA
y la otra por el vector Q[k] = (Qk1 , Qk2 )t . Observar que la matriz de rigidez
[k] [k]
K [k] es simétrica puesto que Kij = Kji .
U1 = u11 ,
Ui = ui−1,2 = ui1 , para i = 2, . . . , N , (24.24)
UN +1 = uN 2 .
que, en el caso de que Qk2 + Qk+1,1 = 0, esta condición implica en dicho nodo
de conexión que la derivada de la solución exacta du(x)/dx es continua debido
a la definición de Qki , ver (24.10).
KU = F , (24.25)
⎛ [1] [1]
⎞
K11 K12
⎜ [1] [1] [2] [2] ⎟
⎜ K21 K22 + K11 K12 ⎟
⎜ [2] [2] [3] [3] ⎟
⎜ K21 K22 + K11 K12 ⎟
⎜ ⎟
⎜ .. .. .. ⎟
K = ⎜ . . . ⎟ ,
⎜ ⎟
⎜ ⎟
⎜ .. .. .. ⎟
⎜ . . . ⎟
⎜ [N −1] [N −1] [N ] [N ] ⎟
⎝ K21 K22 + K11 K12 ⎠
[N ] [N ]
K21 K22
314 Método de los Elementos Finitos
⎛ ⎞
U1
⎜ U2 ⎟
⎜ ⎟
⎜ U3 ⎟
⎜ ⎟
⎜ .. ⎟
⎜ ⎟
U = ⎜ . ⎟
⎜ ⎟
⎜ .. ⎟
⎜ . ⎟
⎜ ⎟
⎝ UN ⎠
UN +1
⎛ [1]
⎞ ⎛ ⎞
F1 Q11
⎜ [1] [2] ⎟ ⎜ ⎟
⎜ F 2 + F1 ⎟ ⎜ Q12 + Q21 ⎟
⎜ [2] [3] ⎟ ⎜ ⎟
⎜ F 2 + F1 ⎟ ⎜ Q22 + Q31 ⎟
⎜ ⎟ ⎜ ⎟
⎜ .. ⎟ ⎜ .. ⎟
F = ⎜ . ⎟+⎜ . ⎟.
⎜ ⎟ ⎜ ⎟
⎜ .. ⎟ ⎜ .. ⎟
⎜ ⎟ ⎜ ⎟
⎜ . ⎟ ⎜ . ⎟
⎜ [N −1] [N ] ⎟ ⎝ ⎠
⎝ F2 + F1 ⎠ QN −1,2 + QN 1
[N ] QN 2
F2
r(cm) 0 2 4 6 8
T (◦ C) 30 35 40 45 50
μ(Kg/ms) 0.40 0.28 0.20 0.12 0.10
Problemas Variacionales
con Fronteras Fijas
25.1. Introducción
Sea E un espacio de funciones y(x) de clase C 1 [a, b] tomando valores fijos
en los extremos, es decir, y(a) = ya , y(b) = yb . Consideremos el funcional
J : E → R, definido de la forma
! b
J(y) = F (x, y, y ) dx . (25.1)
a
El problema que nos planteamos consiste en hallar una función y(x) ∈ E que
sea un extremal (máximo o mı́nimo) del funcional J(y) definido en (25.1). Por
ejemplo, la función y ∗ (x) minimiza al funcional J si se verifica J(y ∗ ) ≤ J(y) para
todo y ∈ E. Por supuesto, el principal problema con el que nos encontramos y
que lo diferencia radicalmente de las técnicas utilizadas en los problemas clásicos
de hallar máximos y mı́nimos de funciones de varias variables es el hecho de que
la dimensión del espacio de funciones E es infinita.
b−a
xi = a + ih , h = , i = 0, 1, . . . , n + 1 .
n+1
317
318 Problemas Variacionales con Fronteras Fijas
Entonces
∂Jn ∂F yi − yi−1 ∂F yi − yi−1
= xi , y i , h + xi , yi ,
∂yi ∂y h ∂y h
∂F yi+1 − yi
− xi+1 , yi+1 , = 0 , i = 1, . . . , n. (25.2)
∂y h
∂fi ∂fj
= .
∂yj ∂yi
∂fi ∂ 2 Jn
=
∂yi−1 ∂yi ∂yi−1
25.3 El problema de la braquistocrona 319
∂2F yi − yi−1 ∂2F yi − yi−1 1
= − xi , y i , − 2 xi , yi , (25.3)
∂y∂y h ∂y h h
∂fi ∂ 2 Jn ∂2F yi − yi−1 ∂2F yi − yi−1
= = xi , yi , h+2 xi , yi ,
∂yi ∂yi2 ∂y 2 h ∂y∂y h
2
2
∂ F yi − yi−1 1 ∂ F yi+1 − yi 1
+ 2 xi , yi , + xi+1 , yi+1 , (25.4)
∂y h h ∂y 2 h h
2 2
∂fi ∂ Jn ∂ F yi+1 − yi
= =− xi+1 , yi+1 ,
∂yi+1 ∂yi ∂yi+1 ∂y∂y h
2
∂ F yi+1 − yi 1
− 2 xi+1 , yi+1 , , (25.5)
∂y h h
1
mv 2 + mgh = mghA ,
2
320 Problemas Variacionales con Fronteras Fijas
de modo que v 2 = 2g(hA − h). Con el sistema de referencia tomado está claro
que x = hA − h, de modo que, la ecuación diferencial que rige el movimiento de
la masa m es
ds
v := = 2gx ,
dt
siendo v el módulo de la velocidad que tiene la partı́cula cuando se encuentra en
un punto con coordenada
(x, y(x)) y ds el elemento de longitud de la trayectoria.
Sabemos que, ds = 1 + y (x)2 dx, de modo que se obtiene
1 + y (x)2 dx
dt = √ .
2gx
Introducción a
Mathematica
26.1. Introducción
En este apéndice se pretende dar unas breves nociones sobre algunos de los
comandos utilizados por el programa Mathematica1 .
Además, la ayuda en Mathematica para conocer información detallada
sobre un cierto comando puede realizarse escribiendo ?? delante del nombre del
comando o bien acudiendo al menú de ayuda.
321
322 Introducción a Mathematica
Una forma de visualizar los elementos de una lista x en forma de una tabla es
mediante el comando T ableF orm[x]. Por ejemplo, la salida que produce Math-
ematica tras el comando T ableF orm[{{1, 2}, {3, 4}}] es
1 2
3 4
26.3. Funciones
Si se desea definir una función real de variable real f : R → R, por ejemplo,
x3 exp(2x)
f (x) = ,
sin x
la forma de realizarlo es
f [x ] := (x ∧ 3 ∗ Exp[2 ∗ x])/Sin[x] .
Para evaluar la función f en un punto, por ejemplo x0 , simplemente se escribe
f [x0 ].
26.4 Representación Gráfica 323
Por supuesto también es posible definir una función real de varias variables
reales f : Rm → R. Por ejemplo, para definir la función
f (x, y) = x3 − 3xy ,
escribirı́amos la sentencia
f [x , y ] := x ∧ 3 − 3 ∗ x ∗ y
f [{x , y }] := x ∧ 3 − 3 ∗ x ∗ y
escribirı́amos la sentencia
f [x , y ] := {x ∧ 3 ∗ y, Cos[x], E ∧ (x ∗ y)} .
Una instrucción muy util que nos permitirá, dada una nube de puntos, re-
presentarlos gráficamente es el comando ListP lot[]. La forma de utilizarlo es la
siguiente. Consideremos una conjunto de puntos en el plano (xi , yi ) ∈ R2 con
i = 1, . . . , n. Entonces introduciendo
se dibujan los puntos en el plano. P lotStyle− > tiene diversas opciones (ver
el manual del programa), nosotros hemos utilizado una opción que nos permite
controlar el taman̄o de los puntos. Si además se utiliza
ListP lot[{{x1 , y1 }, {x2 , y2 }, . . . , {xn , yn }}, P lotJoined− > T rue]
entonces los puntos representados serán unidos por una recta. También es posi-
ble crear gráficas con colores mediante la opción P lotStyle− > Hue[m], con
0 ≤ m ≤ 1, de manera que a medida que aumenta m se pasa por rojo, amarillo,
verde, azul, etc...
IdentityM atrix[n]
463908967738245731201181971316266541335354
606980986542278147922285826879934715156756
753035191305059570813964428091287791299849
47682780134232434682183959939102665927776
439857317019738245619747364381030569622906
549543129198965152594576297260803180404842
596116964782060639065966089166712338842834
343045860606296811709288710361986642298681
053856193817242616463360000000000000000000
000000000000000000000000000000000
Los números racionales son también tratados de forma exacta. Por ejemplo
(2/5 + 7/4) ∧ 2 − 1/2 devuelve 1649/400. Nótese la diferencia con (2/5 + 7/4) ∧
2 − 0,5 que devuelve la aproximación 4,1225.
Respecto a los números irracionales, debido a sus infinitas cifras decimales,
o bien se trabaja con ellos de forma simbólica o bien son aproximados
√ con
un número finito de decimales. Por ejemplo Sqrt[2] devuelve 2. Sin embargo
N [Sqrt[2]] da la aproximación numérica 1,41421 o bien la aproximación con 20
dı́gitos 1,4142135623730950488. Otro ejemplo es la aproximación del número π
con 30 dı́gitos dada por N [P i, 30] que da
3,14159265358979323846264338328
Por ejemplo, si se requiere que el usuario introduzca un número real y este sea
asignado a la variable z, la sintaxis será z = Input[”introduzca un numero real”].
imprime otra vez por pantalla las primitivas de las funciones xi para i =
1, 2, . . . , 5.
i = 1;
W hile[i <= 5, P rint[Integrate[x ∧ i, x]]; i = i + 1]
de la forma
f [x ] = If [x > 0, Log[x], Sin[x]]
En muchas ocasiones, es necesario incluir en las sentencias condicionales
anteriores algún operador lógico. Estos operadores son
Operador Significado
! NO
&& Y
|| O
El valor de n es 14
calculo = {};
F or[n = 1, n <= 20, n + +, AppendT o[calculo, (1. + 1/10 ∧ n) ∧ (10 ∧ n)];
]
calculo
cero debida a errores de redondeo en la resta entre paréntesis. Veamos este efecto
de los errores por cancelación con el siguiente programa.
calculo = {};
F or[n = 1, n <= 20, n + +,
mal = 10 ∧ n ∗ (Sqrt[10 ∧ n + 1.] − Sqrt[10 ∧ n]);
bien = 10 ∧ n/(Sqrt[10 ∧ n + 1.] + Sqrt[10 ∧ n]);
AppendT o[calculo, {mal, bien}];
]
T ableF orm[calculo, T ableHeadings− > {N one, {”M AL”, ”BIEN ”}}]
M AL BIEN
1 54347 1 54347
4 98756 4 98756
15 8074 15 8074
49 9988 49 9988
158 113 158 113
.. ..
. .
5 02914 106 5 106
1 49012 107 1 58114 107
0 5 107
0 1 58114 108
.. ..
. .
0 5 109
Observar que el primer argumento es una lista que contiene las coordenadas de la
nube de puntos y el segundo argumento es el nombre de la variable del polinomio
interpolador. Por ejemplo, la recta en el plano que pasa por los puntos (0, 1) y
(3, 7) se puede obtener de la forma InterpolatingP olynomial[{{0, 1}, {3, 7}}, x],
siendo la solución 1 + 2x.
La sintaxis a utilizar es
Observar que el primer argumento es una lista que contiene las coordenadas de
la nube de puntos, el segundo argumento es una lista con las funciones base y el
tercero el nombre de la variable. Por ejemplo, la recta de regresión en el plano
que pasa por los puntos (xi , yi ) para i = 0, 1, . . . , n se puede obtener de la forma
26.7.4. Integración
Con el comando N Integrate[f, {x, a, b}] se obtiene una aproximación numéri-
ca de la integral
! b
f (x) dx .
a
devuelve 1 91173.
1 2
Para la integral 0 ex dx, el comando Integrate[E ∧ (x ∧ 2), {x, 0, 1}]
√
de Mathematica devuelve 12 πErf i[1], es decir, utiliza la función error
iz
imaginaria Erf i[] que está definida como Erf i[z] := i√2π 0 exp(−t2 ) dt.
Sin embargo N Integrate[E ∧ (x ∧ 2), {x, 0, 1}] da el valor aproximado
1 46265.
o bien
N Solve[{p1 == 0, . . . , pm == 0}, {x1 , . . . , xm }]
respectivamente. Ası́, los puntos de corte entre la circunferencia unidad x2 +y 2 =
1 y la recta bisectriz del segundo cuadrante y = −x se calculan de forma exacta
mediante Solve[{x∧2+y∧2−1 == 0, x+y √ == 0}, {x, y}]
√ obteniéndose√las pare-
jas de
√ coordenadas en la lista {{x− > −1/ 2, y− > 1/ 2}, {x− > 1/ 2, y− >
−1/ 2}}. De forma numérica se tiene N Solve[{x ∧ 2 + y ∧ 2 − 1 == 0, x + y ==
0}, {x, y}] que produce por resultado la aproximación {{x− > −0707107, y− >
0 707107}, {x− > 0 707107, y− > −0 707107}}.
N Solve[f == 0, x]
F indRoot[{f1 , . . . , fm }, {y1 , . . . , ym }]
dy
= f (x, y) .
dx
Mathematica intenta hallar la solución general exacta y(x) mediante el co-
mando DSolve[y [x] == f [x, y[x]], y[x], x]. Por ejemplo, para hallar la solución
general exacta de la ecuación diferencial lineal de primer orden dy/dx = y + 2x,
mediante la instrucción DSolve[y [x] == y[x] + 2 ∗ x, y[x], x] se obtiene
dy
= f (x, y) , y(x0 ) = y0 , (26.1)
dx
entonces la condición inicial y(x0 ) = y0 se introduce de la forma DSolve[{y [x] ==
f [x, y[x]], y[x0 ] == y0 }, y[x], x]. De este modo, la solución exacta del prob-
lema de valor inicial dy/dx = y + 2x con y(0) = 1 mediante el comando
DSolve[{y [x] == y[x] + 2 ∗ x, y[0] = 1}, y[x], x] es
expresada a través de la función error imaginaria Erf i[] que se ha visto ante-
riormente que es una expresión integral y no es por lo tanto de ayuda. Es en
este tipo de casos donde los métodos numéricos tienen sentido de ser aplicados.
Mathematica devuelve
2.5
1.5
2 4 6 8 10
Introducción a Octave
Las instrucciones exit o bien quit cierra el intérprete mientras que <CTRL>-c
para la ejecución de los comandos. El comando help seguido del nombre de una
función permite acceder a la ayuda asociada a una función.
337
338 Introducción a Octave
La instrucción who o bien whos devuelve una lista con las variables almace-
nadas en la sesión de trabajo. Una variable variable se elimina de la memoria con
el comando clear variable. Si se usa sólo clear se borran todas las variables
de la memoria.
El comando save seguido de un nombre de archivo y un conjunto de vari-
ables sirve para guardar el archivo en el llamado directorio de trabajo (worwing
directory) con el valor de las variables especificadas. Posteriormente, en otra
sesion, utilizando el comando load se pueden recuperar todas las variables al-
macenadas. Por ejemplo, save sesion v1 v2 v3 guarda en el archivo sesion el
valor que tengan almacenadas las variables v1, v2 y v3.
Operación Sintaxis
x±y x±y
xy x*y
x
y x/y
y
x x∧y
Algunas funciones elementales con números reales x e y son:
Función Sintaxis Descripción
|x| abs(x) valor absoluto
sign(x) sign(x) -1 si x < 0; 0 si x = 0; 1 si x > 0
rem(x/y) rem(x,y) resto del cociente x/y
ex exp(x) exponenial
ln x log(x) logaritmo neperiano de x > 0
log10 x log10(x) logaritmo decimal de x > 0
√
x sqrt(x) raı́z cuadrada de x ≥ 0
sin(x) sin(x) seno
cos(x) cos(x) coseno
tan(x) tan(x) tangente
arcsin(x) asin(x) arco seno de x ∈ [−1, 1] da ángulo en [−π/2, π/2]
arc cos(x) acos(x) arco coseno de x ∈ [−1, 1] da ángulo en [0, π]
arctan(x) atan(x) arco tangente da ángulo en [−π/2, π/2]
Recordar que los ángulos en los argumentos de las funciones trigonométricas
han de estar siempre en radianes.
Función Descripción
length(x) Devuelve n si x ∈ Rn
dot(x,y) Produto escalar x.y
norm(x) Norma euclideana x
x’ Traspuesto conjugado
max(x) Valor máximo de las componentes de x
min(x) Valor mı́nimo de las componentes de x
sum(x) Suma todas las componentes de x
prod(x) Multiplica todas las componentes de x
Sea A ∈ Mn×m (R). Por lo que se refiere al orden de las matrices, se dispone
de la función size(A) que devuelve el vector fila [n, m] con las filas y columnas
de A.
Para definir rápidamente matrices que contienen ceros y unos se dispone,
entre otras, de las siguientes funciones:
Función Descripción
zeros(n) Matriz nula cuadrada de orden n
zeros(n,m) Matriz nula de orden n × m
ones(n) Matriz de unos cuadrada de orden n
ones(n,m) Matriz de unos de orden n × m
eye(n) Matriz identidad de orden n
rand(n,m) Matriz de orden n × m con elementos aleatorios
uniformemente distribuidos en el intervalo [0, 1]
Por ejemplo, zeros(size(A)) crea una matriz de elementos nulos del mismo
orden que la matriz A.
Las operaciones básicas entre matrices A y B de órdenes adecuados para
poder realizarse son:
Operación Descripción
a*A Producto del escalar a por la matriz A
A±B Suma y resta
A*B Multiplicación
A∧n Potencia
A’ Traspuesta conjugada
En Matlab–Octave existe lo que se conoce como operaciones entre matrices
elemento a elemento y son de gran utilidad. Por supuesto, las matrices deben
tener el orden adecuado para poder realizarse la operación. Ası́, si A = (aij ), B =
(aij ) ∈ Mn×m (R) y k, l ∈ R se tiene:
Operación Descripción
C=k*A ± l C = (cij ) tal que cij = kaij ± l
C=A.*B C = (cij ) tal que cij = aij bij
C=A./B C = (cij ) tal que cij = aij /bij
b
C=A.∧B C = (cij ) tal que cij = aijij
C=A.∧k C = (cij ) tal que cij = akij
342 Introducción a Octave
Función Descripción
size(A) Vector fila [n, m] si A ∈ Mn×m (R)
det(A) Determinante de A
rank(A) Rango de A
inv(A) Inversa de A
x=A\b Solución del sistema lineal Ax = b
eig(A) Vector columna con los valores propios de A
[P,D]=eig(A) P es una matriz con los vectores propios de A en columnas
D es una matriz diagonal con los valores propios de A en la diagonal
poly(A) Polinomio caracterı́stico de A
expm(A) Matriz exponencial de A
27.1.4. Funciones
Una de las herramientas más potentes de Matlab–Octave es el concepto de
función. A parte de las funciones elementales habituales en un curso de cálculo,
como son las trigonométricas, logaritmos, etc..., Matlab–Octave dispone de una
enorme biblioteca de funciones matemáticas menos conocidas. Por ejemplo la
función gamma de Euler definida de la forma
! ∞
Γ(z) = tz−1 exp(−t) dt ,
0
function [vs1,...,vsm]=nombre(ve1,...,ven)
...comandos...
endfunction
function [M]=media(x)
M=sum(x)/length(x);
endfunction
x3 exp(2x)
f (x) = ,
sin x
la forma de realizarlo es
function y=f(x)
y=(x ∧ 3*exp(2*x))/sin(x);
endfunction
Por supuesto también es posible definir una función real de varias variables
reales g : Rm → R. Por ejemplo, para definir la función
g(x, y) = x3 − 3xy ,
escribirı́amos
function z=g(x,y)
z=x∧3-3*x*y;
endfunction
Para evaluar la función g en un punto, por ejemplo (1, 2), se escribe g(1,2)
y Matlab-Octave devuelve -5.
344 Introducción a Octave
function y=h(x)
y(1)=x(1)∧3*x(2);
y(2)=cos(x(1));
y(3)=exp(x(1)*x(2));
endfunction
f=inline(’sin(x)*sqrt(x)’);
f(1)
f=inline(’sin(x)*sqrt(x)’);
feval(f,1)
27.1.5. Gráficos
Una de las mayores ventajas de Matlab–Otave es su capacidad gráfica para
mostrar resultados. Se generan gráficas 2D y 3D. Los comandos básicos son:
27.1 Introducción a Matlab–Octave 345
fplot para dibujar funciones reales de variable real; plot para dibujar puntos
en R2 y plot3 para dibujar puntos en R3 .
f=inline(’sqrt(x)’);
fplot(f,[0,1])
√
muestra la gráfica de la función f (x) = x en el intervalo [0, 1]
x=[-pi:0.01:pi];y=cos(x); plot(x,y)
Nota: Observar cómo, con el puntero del ratón colocado sobre la ventana
gráfica, apretando el botón izquierdo y desplazando, se obtienen perspectivas
diferentes de la gráfica 3D. Se puede conseguir el mismo efecto con las teclas de
flechas del cursor.
En los gráficos se puede poner un tı́tulo y etiquetar los ejes. Por ejemplo,
el código siguiente genera la gráfica etiquetada de la parábola y = x2 − 2 en el
intervalo [−5, 5].
GRAFICA DE LA PARABOLA
25
line 1
20
15
y=x2-2
10
-5
-6 -4 -2 0 2 4 6
x
27.1.7. Polinomios
Sea p(x) ∈ Rn [x] un polinomio a coeficientes reales de grado n, es decir,
n
f (x) = i=0 ai xi con ai ∈ R. En Matlab–Octave, los polinomios se introducen
como vectores fila de Rn+1 que contienen los coeficientes del polinomio p en
orden descendente. De este modo, se debe introducir p=[an , an−1 , . . . , a1 , a0 ].
Para evaluar el polinomio en un punto, es decir, calcular p(c) con c ∈ R, se
utiliza el comando polyval(p,c). La derivada de un polinomio p (x) se obtiene
mediante la instrucción polyder(p) o polyderiv(p). El polinomio F (x) que es
x
la primitiva del polinomio p(x) tal que F (0) = 0, es decir, F (x) = 0 p(s) ds,
se calcula mediante la instrucción polyinteg(p).
Dados dos polinomios p(x) ∈ Rn [x] y q(x) ∈ Rm [x], el polinomio producto
p(x)q(x) ∈ Rn+m [x] se obtiene de la forma conv(p,q). Veamos un ejemplo con
los polinomios p(x) = x3 − x + 1 y q(x) = −x + 2:
p=[1,0,-1,1]; q=[-1,2];
polyval(p,4)
polyder(p)
polyinteg(p)
conv(p,q)
41
3 0 -1
0.25000 0.00000 -0.50000 1.00000 0.00000
-1 2 1 -3 2
x=1;
disp(’El valor de x es:’)
disp(x)
El valor de x es:
1
Recordemos que, los comandos save y load permiten escribir y leer del disco
duro. Por ejemplo, podemos tabular la función f (x) = sin(x) en el intervalo
[−10, 10] con una longitud de paso h = 1 y guardar en el directorio de trabajo la
tabulación obtenida en un fichero llamado tabulacion.dat de la forma siguiente:
load tabulacion.dat
.....proceso......
endfor
for i = 2:2:8
disp(i), disp(’es multiplo de 2’)
endfor
Matlab-Octave devuelve
2
es multiplo de 2
4
es multiplo de 2
6
es multiplo de 2
8
es multiplo de 2
while(condicion)
.....proceso......
endwhile
a=0;
while(a <3)
disp(a)
a=a+1;
endwhile
0
1
1 También es posible una estructura del tipo do--until.
350 Introducción a Octave
if(condicion)
.....proceso......
endif
if(condicion 1)
.....proceso 1......
elseif(condicion 2)
.....proceso 2......
else
.....proceso 3......
endif
if(n <= 0)
p=0;
elseif(rem(n,2)==0)
p=2;
else
p=1;
endif
Operador Significado
! NO
& Y
| O
Operador Significado
< menor que
> mayor que
<= menor o igual que
>= mayor o igual que
== igual que
! = o bien ∼= distinto que
27.2.3. Errores
Debido a errores de redondeo, se define el epsilon de la máquina como el
mayor número positivo tal que 1 + x = 1 para todo x ∈ (0, ). Si intentamos
averiguar con un programa escrito en el lenguaje de Matlab–Octave podrı́amos
intentar lo siguiente.
x=1;
while(x+1>1)
x=x/2;
endwhile
disp(’El epsilon de la maquina es:’)
disp(2*x)
for i=1:10
n=100∧i; disp((1+1/n)∧n)
352 Introducción a Octave
endfor
la salida es la siguiente:
2.7048
2.7181
2.7183
2.7183
2.7183
2.7185
2.7161
1
1
1
for n=1:10
x=10∧n;
mal=x*(sqrt(x+1)-sqrt(x)); bien=x/(sqrt(x+1)-sqrt(x));
A(n,:)=[mal,bien];
endfor
segundos) de ejecución.
flops(0), x = A \ b; flops
flops(0), x=inv(A)*b; flops
o bien de la forma
tic, x = A \ b; toc
tic, x=inv(A)*b; toc
Existen otras formas de operar. Por ejemplo, vander(x) \ y’, siendo x=[x0 ,
. . . , xn ], y=[y0 , . . . , yn ] los vectores fila con las abscisas y ordenadas, respecti-
vamente, de la nube de puntos. El resultado es el polinomio Pn (x) dado como
un vector columna que contiene los coeficientes de Pn (x) en orden descendente.
354 Introducción a Octave
Por ejemplo, la recta en el plano que pasa por los puntos (0, 1) y (3, 7) es
y = P1 (x) = 1 + 2x. Se puede hallar de la forma vander([0, 3]) \ [1,7]’,
siendo el resultado el vector columna
2
1
10
line 1
line 2
9
1
1 1.5 2 2.5 3 3.5 4 4.5 5
quad(inline(’exp(-x∧2)’),0,1)
y se obtiene 0,74682.
Puede aproximar numéricamente
√ integrales impropios. Por ejemplo, se sabe
∞
que 0 exp(−x2 ) dx = π/2. Entonces
quad(inline(’exp(-x∧2)’),0,inf)
ẋ = f (t, x) , x(t0 ) = x0 ,
function xp=f(x)
mu=1;
xp=zeros(2,1);
xp(1)=x(2);
xp(2)=mu*x(2)*(1-x(2)∧2)-x(1);
endfunction
x0=[2,0];
tau=linspace(0,20,100);
sol=lsode(’f’,x0,tau);
Se obtiene una matriz 100 × 2 llamada sol, donde la fila i-ésima contiene una
aproximación del vector (x(ti ), y(ti )). Podemos dibujar la solución x(t) en el
intervalo t ∈ [0, 20] de la forma
x=sol(:,1);
xlabel(’t’) , ylabel(’x(t)’) ,
title(’Ecuacion de van der Pol con mu = 1’) ,
plot(tau,x)
1.5
0.5
x(t)
-0.5
-1
-1.5
0 5 10 15 20
t
27.3 Métodos Numéricos 357
p=[1,-5,13,-19,10]; x=roots(p)
x=
1.0000 + 2.0000i
1.0000 - 2.0000i
2.0000
1.0000
16
line 1
line 2
14
12
10
-2
0 0.5 1 1.5 2 2.5 3
Ecuaciones Trascendentes
Consideremos a partir de ahora una ecuación trascendente f (x) = 0. El
comando fzero(f,x0 ) intenta hallar una solución aproximada de la ecuación
trascendente f (x) = 0 comenzando la búsqueda en el punto inicial x0 . Si no se
converge hacia una solución, devuelve por respuesta NaN. También se puede uti-
lizar de la forma fzero(f,[x0 , x1 ]) si sabemos que la raı́z buscada se encuentra
en el intervalo [x0 , x1 ] ⊂ R.
Ası́, supongamos que f (x) = exp(x) − sin(x). Observar que, en el ejemplo
propuesto, existen infinitas raı́ces negativas y ninguna positiva mediante un
simple análisis gráfico de las funciones involucradas. Entonces el código
f=inline(’exp(x)-sin(x)’);
raiz=fzero(f,1)
f=inline(’exp(x)-sin(x)’);
raiz=fzero(f,-5)
function y=f(x)
y(1)=sin(x(1))-x(2)+2;
y(2)=x(1)∧2+tan(x(1)*x(2))-1;
endfunction
raiz=fsolve(’f’,[0;1])
y se obtiene
raiz=
1.05524
2.87002
Por ejemplo, la recta de regresión en el plano que pasa por los puntos (0, 1),
(2, 3) y (3, 5) se puede obtener de la forma:
x=[0,2,3]; y=[1,3,5];
p=polyfit(x,y,1)
1.28571
0.85714
360 Introducción a Octave
xx=[0:0.1:3];
pp=polyval(p,xx);
plot(x,y,’*’)
hold on
plot(xx,pp)
5
line 1
line 2
4.5
3.5
2.5
1.5
0.5
0 0.5 1 1.5 2 2.5 3
[11] P. Henrici. Elements of numerical analysis. John Wiley, New York, 1964.
361
362 BIBLIOGRAFÍA