FÍSICA COMPUTACIONAL II (GRADO EN FÍSICAS) CÓDIGO: 61042047

PRIMERA SEMANA, PRUEBAS PRESENCIALES FEBRERO DE 2014,

DURACIÓN: 2 HORAS

INSTRUCCIONES: Puede utilizar todo tipo de material auxiliar y calculadora progra-

mable. No es suficiente poner la solución, en cada ejercicio debe especificar y exponer
los pasos intermedios realizados.

1. (2.5 puntos) Considere el siguiente problema de valores iniciales

2z ′′ (t) + 4z ′ (t) − 2z(t) = sin[z(t)] z(t = 0) = 0 z ′ (t = 0) = 10

Obtenga la solución de esta ecuación utilizando el método de predicción

corrección de Euler hasta t = 0.1 utilizando un paso ∆t = 0.02.

Solución:

Sea y(t) = z ′ (t). Entonces la ecuación diferencial de segundo orden se transforma al

sistema equivalente de dos ecuaciones diferenciales de primer grado

z ′ (t) = y(t) z(t = 0) = 0

sin[z(t)]
y ′ (t) = − 2y(t) + z(t) y(t = 0) = 10
2
Sea la función g(t, y(t), z(t)) = sin[z(t)]
2 − 2y(t) + z(t). El método iterativo parte de las
condiciones
t0 = 0 z(t0 ) = 0 y(t0 ) = 10
A continuación avanzamos el tiempo para el ∆t = 0.02 dado en el enunciado

t1 = t0 + ∆t

y calculamos para este tiempo secuencialmente y de manera alternada, es decir, prime-

ro calcularemos z y posteriormente y, las magnitudes que necesitamos para resolver la
ecuación. Para el método de Euler de predicción-corrección las magnitudes predichas
serán

zp (t1 ) = z(t0 ) + ∆ty(t0 ) = 0 + 0.02 · 10 = 0.2

 
sin[z(t0)]
yp (t1 ) = y(t0) + ∆t − 2y(t0 ) + z(t0 )
2
 
sin[0]
= 10 + 0.02 · − 2 · 10 + 0 = 9.6
2
Pasemos ahora a corregir los valores. Empezamos por la variable z, como z ′ = y, tenemos
y(t0) + yp (t1 ) 10 + 9.6
zc (t1) = z(t0 ) + ∆t = 0 + 0.02 = 0.196
| 2
{z } 2
yav

en la que se ha tomado el valor medio de la derivada z ′ = yav en el punto de partida del

intervalo t0 y la predicción de dicha derivada en el punto final del intervalo yp (t0 ) como
y(t )+y (t )
yav = 0 2 p 1 .

Pasemos ahora a corregir los valores de yc (t1). Recordemos la definición de y ′ (t) =

g(t, y(t), z(t)) con g(t, y(t), z(t)) = sin[z(t)]
2 − 2y(t) + z(t).

yc (t1 ) = y(t0 ) + ∆tgav

Veamos que podemos estimar gav de dos formas distintas:

a) Tomar como valor promedio gav la media aritmética de g(t0 , y(t0), z(t0 )) en el
punto inicial del intervalo y en el punto final los valores predichos g(t1 , yp(t1 ), zp(t1 ))
g(t0, y(t0 ), z(t0)) + g(t1 , yp(t1 ), zp(t1))
yc (t1) = y(t0 ) + ∆t
2

b) Tomar como valor promedio gav la media aritmética de g(t0, y(t0 ), z(t0))
en el punto inicial del intervalo y en el punto final utilizando el mejor valor
disponible, es decir, los valores predichos y el valor corregido para la variable z,
dado que ya lo tenemos calculado, es decir, g(t1 , yp (t1), zc (t1))
g(t0, y(t0), z(t0 )) + g(t1 , yp(t1 ), zc(t1))
yc (t1) = y(t0 ) + ∆t
2

Evidentemente ambos resultados diferirán. Veamos los resultados concretos para este
problema:
a)
g(t0 , y(t0), z(t0)) + g(t1 , yp (t1), zp(t1 ))
yc (t1 ) = y(t0) + ∆t
 2   
∆t sin[z(t0)] sin[zp (t1)]
= y(t0) + − 2y(t0 ) + z(t0 ) + − 2yp (t1) + zp (t1 )
2 2 2
   
0.02 sin[0] sin[0.2]
= 10 + − 2 · 10 + 0 + − 2 · 9.6 + 0.2
2 2 2
0.1986693
= 10 + 0.01(−20 + − 19.2 + 0.2) = 9.61099 (1)
2
A continuación avanzamos el tiempo, es decir t0 = t1 . El proceso se repetirá desde el
cálculo de zp (t1 ) avanzando en el tiempo tantas veces como sea necesario para que t0
llegue a tomar el valor que se nos indica como final (t = 0.1) en el enunciado.

La siguiente tabla de datos recoge los valores que se obtienen al hacer la iteración de la
anterior ecuación hasta el tiempo t = 0.1 con el algoritmo a).

t zp (t) yp (t) zc (t) yc (t)

0.00E+00 0.0000000E+00 0.1000000E+02 0.0000000E+00 0.1000000E+02
0.20E-01 0.2000000E+00 0.9600000E+01 0.1960000E+00 0.9610993E+01
0.40E-01 0.3882198E+00 0.9232421E+01 0.3844341E+00 0.9242833E+01
0.60E-01 0.5692908E+00 0.8884559E+01 0.5657080E+00 0.8894393E+01
0.80E-01 0.7435959E+00 0.8555291E+01 0.7402049E+00 0.8564557E+01
0.10E+00 0.9114960E+00 0.8243523E+01 0.9082857E+00 0.8252237E+01

b)
g(t0 , y(t0), z(t0)) + g(t1 , yp(t1 ), zc(t1 ))
yc (t1) = y(t0 ) + ∆t
 2   
∆t sin[z(t0 )] sin[zc (t1)]
= y(t0 ) + − 2y(t0 ) + z(t0 ) + − 2yp (t1) + zc (t1)
2 2 2
   
0.02 sin[0] sin[0.196]
= 10 + − 2 · 10 + 0 + − 2 · 9.6 + 0.196
2 2 2
0.194747
= 10 + 0.01(−20 + − 19.2 + 0.196) = 9.610933 (2)
2
La siguiente tabla de datos recoge los valores que se obtienen al hacer la iteración de la
anterior ecuación hasta el tiempo t = 0.1 con el algoritmo b).
 t zp (t) yp (t) zc (t) yc (t)
0.00E+00 0.0000000E+00 0.1000000E+02 0.0000000E+00 0.1000000E+02
0.20E-01 0.2000000E+00 0.9600000E+01 0.1960000E+00 0.9610933E+01
0.40E-01 0.3882186E+00 0.9232363E+01 0.3844329E+00 0.9242721E+01
0.60E-01 0.5692873E+00 0.8884451E+01 0.5657046E+00 0.8894234E+01
0.80E-01 0.7435893E+00 0.8555139E+01 0.7401984E+00 0.8564357E+01
0.10E+00 0.9114855E+00 0.8243331E+01 0.9082753E+00 0.8252003E+01

2. (2.5 puntos) a) Sabemos que usando una cuadratura de n puntos se puede

obtener la integral de un polinomio de grado n − 1 de forma exacta en
un determinado intervalo x ∈ [a, b]. Supongamos que queremos fabricar
una cuadratura de tres puntos (x1, x2 y x3) que nos permita obtener las
integrales de 1, x y x2 de forma exacta en el intervalo x ∈ [0, 1] basada
en los puntos x1 = 0.1, x2 = 0.5 y x3 = 0.9. ¿Cuáles serían los valores de
los pesos que deberíamos poner en cada punto? Ayuda: Queremos saber
cuáles son los pesos ω1 , ω2 y ω3 que satisfacen que la integral numérica
de 1, x, x2 sea exacta en ese intervalo x ∈ [0, 1].
b) Con los pesos encontrados en el apartado anterior, obtenga la siguiente
integral, Z 1
1
I= dx
0 1+x

c) Recalcule la integral anterior usando los pesos y puntos usuales de

una cuadratura gaussiana de 3 puntos. (Nota: Los puntos y pesos de una
cuadratura gaussiana están definidos en el intervalo x ∈ [−1, 1] y son:
x1 = −0.77459667, x2 = 0, x3 = 0.77459667, ω1 = 5/9, ω2 = 8/9, ω3 = 5/9.)
d) Calcule el error cometido en cada uno de los métodos con respecto a
la integral resuelta de forma analítica.

Solución:
a) El enunciado pide que expresemos la integral de acuerdo con la siguiente suma:
Z 1
f (x)dx = ω1 f (x1) + ω2f (x2) + ω3 f (x3)
0

dónde: x1 = 0.1, x2 = 0.5 y x3 = 0.9. Podemos calcular los valores de los pesos (ω1 , ω2
y ω3 ) imponiendo la condición de la exactitud de la integral de los tres monomios 1, x y
x2 que son base para cualquier polinomio de grado 2:
Z 1
si f (x) = 1 : 1 dx = 1 = ω1 f (x1) + ω2 f (x2) + ω3f (x3) → ω1 + ω2 + ω3 = 1
0
 Z 1
1 1
si f (x) = x : x dx = = ω1 f (x1)+ω2f (x2)+ω3f (x3) → 0, 1ω1+0, 5ω2+0, 9ω3 =
0 2 2
Z 1
2 1 1
si f (x) = x : = ω1f (x1)+ω2f (x2)+ω3f (x3) → 0, 01ω1+0, 25ω2+0.81ω3 =
1 dx =
0 3 3
Es decir, que para fijar los valores de los 3 pesos, debemos resolver el siguiente sistema
de ecuaciones.

ω1 + ω2 + ω3 = 1
1
0.1ω1 + 0.5ω2 + 0.9ω3 =
2
1
0.01ω1 + 0.25ω2 + 0.81ω3 =
3
Escribimos el sistema de ecuaciones como:
    
1 1 1 ω1 1
 0.1 0.5 0.9   ω2  =  1 
2
1
0.01 0.25 0.81 ω3 3

Y la solución para cada uno de los coeficientes es:

25
ω1 = = 0.2604166656
96
23
ω2 = = 0.47916668
48
25
ω3 = = 0.2604166656
96

b) Utilizando los resultados del apartado anterior, podemos expresar la solución de la

integral como:
Z 1
1
dx ≈ ω1 f (x1) + ω2f (x2) + ω3 f (x3)
0 1 + x
0.2604166656 0.47916668 0.2604166656
= + + = 0.693248
1 + 0.1 1 + 0.5 1 + 0.9

c) Transformamos la integral para que, usando los pesos definidos en en intervalo x ∈

[−1, 1] podamos calcular la integral en x ∈ [0, 1].

t+1 dt
El cambio de variable adecuado es x = 2
, de donde dx = 2
, con lo que la integral
deviene
Z 1 Z 1   Z 1 Z 1
t + 1 dt 1 dt 1
f (x)dx ≡ f = t+1 = dt
0 −1 2 2 −1 1 + 2 2 −1 t + 3
X 1 1 1 1
≈ ωi = ω1 + ω2 + ω3
i
ti + 3 t1 + 3 t2 + 3 t3 + 3
5 1 8 1 5 1
= q + + q = 0.6931216
9− 3 + 3 90 + 3 9 3 + 3
5 5

d) La integral analítica es
Z 1
1
ITeorica = dx = ln(1 + x)|10 = ln 2 − ln 1 = ln 2 ≈ 0.693147180559
0 1+x

El error absoluto que se obtiene al calcular la integral con la cuadratura que se pidió en
el apartado b) es:

∆ε = INum
b − ITeorica = 0.693248 − 0.6931471 = 0.00010

o lo que es igual un error relativo del 0.014 %.

Y el error absoluto que se comete al calcularlo con los puntos que corresponden a la
cuadratura gaussiana del apartado c) es:

∆ε = INum
c − ITeorica = 0.6931216 − 0.6931471 = 0.000025

o que equivale a un error relativo del 0.0036 %.

Recuerde que si ha realizado alguna actividad de evaluación continua

(ya sea PEL o PED) no debe contestar a ninguno de los siguientes
problemas.

3. (2.5 puntos) Resuelva la siguiente ecuación no lineal mediante el método

de Newton comenzando la busqueda en x0 = 0.5.


ln (x) − exp −x2 = 0

Haga las iteraciones que sean necesarias hasta que la diferencia entre los
valores de x sea menor que 10−5.
Solución:

Hacemos

f (x) = ln (x) − exp −x2 .
para aplicar el método de Newton tenemos que calcular
f (xn)
xn+1 = xn − .
f ′ (xn )
siendo
1

f ′ (xn) = + 2xn exp −x2n
xn
por lo que

ln (xn ) − exp −x2n
xn+1 = xn − 1 2)
.
xn + 2x n exp (−x n

Realizamos la iteraciones para obtener:

Iteracion xn xn+1 f (xn)
1 0.5 1.0297062 −0.31708
2 1.0297062 1.2179474 −0.0296991
3 1.2179474 1.2395675 −0.000364294
4 1.2395675 1.2398393 −5.70023·10−8
5 1.2398393 1.2398394 −1.44329·10−15

4. (2.5 puntos) Resuelva el siguiente sistema de ecuaciones mediante el mé-

todo de eliminación gaussiana. Especifique todos y cada uno de los pasos.

2x + 3y + 2z = 0
5x + y − z = 1
x+y+z = 1

Solución:

Escribimos la matriz de los coeficientes aumentada, y realizamos las siguientes transfor-

maciones (debe indicarse explícitamente lo que se está haciendo en cada uno de los pasos,
por ejemplo, tal y como está marcado en la primera flecha)
     
2 3 2 0 F1 2 3 2 0 2 3 2 0
F3 : F3 − 2
 5 1 −1 1  −− −−−−→  5 1 −1 1  −→  0 −13/2 −6 1 
1 1 1 1 0 −1/2 0 1 0 −1/2 0 1
     
2 3 2 0 2 3 2 0 2 3 2 0
−→  0 −13/2 −6 1  −→  0 −13/2 −6 1  −→  0 −13/2 −6 1 
0 0 6/13 12/13 0 0 1 2 0 0 1 2
     
2 3 2 0 2 3 0 −4 2 3 0 −4
−→  0 −13/2 0 13  −→  0 −13/2 0 13  −→  0 −1 0 2 
0 0 1 2 0 0 1 2 0 0 1 2
   
2 0 0 2 1 0 0 1
−→  0 −1 0 2  −→  0 1 0 −2 
0 0 1 2 0 0 1 2

Siendo la solución    
x 1
 y  =  −2 
z 2