0.

1 Ejercicios de Scilab

1. Escriba una función [A, info] = Heron(a, b, c) que calcula, usando la fórmula de Herón de Ale-
jandrı́a, el área de un triángulo cuyos lados miden a, b, c. El parámetro de salida info valdrá 1 si los
datos están bien.
En esta y en todas las funciones que usted escriba:

• Si los datos pueden ser inadecuados o si el proceso puede no acabar satisfactoriamente debe haber
un parámetro de salida que lo indique, por ejemplo info.
• En los comentarios iniciales debe quedar claro lo que hace la función, el método utilizado, el signifi-
cado de los parámetros de entrada (datos) y de salida (resultados).
• Cuando se presente algún error no previsto, además de influir en el parámetro info, la función debe
escribir un aviso indicando el tipo de error, printf(... ).

2. Escriba una función d = dist(x, y) que calcula la distancia entre dos puntos.

3. Escriba una función a = area(x, y, z) que calcula el área de un triángulo cuyos vértices son x, y, y z.

4. ¿Cual es el área del triángulo cuyos vértices están dados por los 6 últimos dı́gitos de su documento de
identidad?

5. Escriba una función [xmin, fmin] = minimo(f, a, b, h) que halla el minimizador de una función en
un intervalo, simplemente comparando todos los valores f (a), f (a + h), f (a + 2h), ... Muestre el resultado
de su utilización en un ejemplo especı́fico.

6. Escriba una función [r, info] = raiz(f, a, b, h, eps) que halla una raı́z de la función en el inter-
valo, simplemente averiguando si entre los los valores f (a), f (a + h), f (a + 2h), ... hay alguno tal que
|f (x)| ≤ ε. Muestre el resultado de su utilización en un ejemplo especı́fico.

7. Escriba una función [a, info] = ordenaSegun(a, j) que ordena de manera ascendente las filas de
una matriz según los elementos de la columna j. Muestre el resultado de su utilización en un ejemplo
especı́fico.

1
0.2 Procesos inestables

1. En un archivo .sci escriba una función function y = redDec(x, n) que redondea a x con n cifras
decimales.

2. Verifique si la siguiente función redondea un número con n cifras significativas.

function y = fun_xyz(x, n)
y = 0.0
if x == 0, return, end

if x >= 0
signo = 1
else
signo = -1
x = -x
end

m = ceil( log10(x) )
k = 10^m
t = x/k
c = 10^n
y = round(t*c)/c
y = y*k*signo
endfunction

Pruébela varias veces, con número positivos, negativos, con cero, con números cercanos a cero, con
número grandes, con números negativos de valor absoluto grande, ...

3. Escriba una función function z = sumaRed(x, y, n) que redondea a x y a y, efectúa la suma x + y

y redondea el resultado con n cifras significativas.

4. Escriba funciones
function z = prodRed(x, y, n)
function z = restaRed(x, y, n)
function z = divRed(x, y, n)
function z = raiz2Red(x n)
que efectúan el producto, la resta, la división y la raiz cuadrada con n cifras significativas.

5. Considere la sucesión definida por

x0 = 1
x1 = 1/3
13 4
xn = xn−1 − xn−2 , n ≥ 2. (1)
3 3

Demuestre por inducción que

1
xn = , n = 0, 1, 2, ... (2)
3n
6. Haga un progama, archivo .sce, que produzca una tabla con cinco columnas:

2
 i) n, de 0 a 20.
ii) xn según (2).
iii) xn según (1) trabajando con todas las cifras que Scilab usa normalmente.
iv) xn según (1) trabajando con 10 cifras significativas.
v) xn según (1) trabajando con 5 cifras significativas.

Los números decimales deben aparecer con 15 cifras decimales (puede usar formato %20.15f.

7. Comente los resultados de la tabla, concluya.

3
0.3 Método de Gauss sin y con pivoteo

El objetivo es escribir una función que implemente el método de Gauss sin pivoteo, una función que im-
plemente el método de Gauss con pivoteo parcial y comparar los resultados. En ambos casos, en todas las
operaciones se redondea con un número fijo de cifras. Puede ser con un número fijo de cifras significativas
o con un número fijo de cifras decimales.

1. Escriba funciones, redond(x, n), sumaRed(x, y, n), multRed(x, y, n), divRed(x, y, n),
restaRed(x, y, n), redMatr(a, n), para redondear un número con n cifras, sumar (redondea los
números, suma y redondea el resultado), multiplicar, dividir, restar y redondear una matriz.

2. Escriba funciones
[x, info] = solTriSupR(a, b, nc)
[a, b, info] = triangR(a0, b0, nc)
[x, info] = GaussR(a, b, nc)
que, respectivamente, resuelve un sistema triangular superior, triangulariza un sistema de ecuaciones y
resuelve un sistema de ecuaciones por el método de Gauss. En las tres funciones las operaciones se hacen
utilizando redondeo con nc cifras. Siempre, por convención, el parámetro de salida o resultado info
valdrá 1 si y solamente si el proceso se desarrolló satisfactoriamente. No olvide incluir en la primera
función (en el sentido de su utilización en el tiempo) el redondeo de los datos.
Puede ensayar la función con solTriSupR con ejemplos construidos “a mano” o con datos aleatorios. Por
ejemplo, puede usar órdenes semejantes a

n = 4;
a = (rand(n,n) - 0.5)*20; a = triu(a);
x0 = (rand(n,1) - 0.5)*20; b = a*x0;
xTS = solTriSup0(a, b); errTS = norm(xTS-x0)

El valor de errTS debe ser cero o casi cero. Para GaussR puede utilizar órdenes semejantes a las
anteriores, excluyendo a = triu(a).
Para escribir las 3 funciones, pueden ser útiles pero no indispensables, los siguientes pasos:

• escribir funciones x = solTriSup0(a, b), [a, b] = triang0(a0, b0), x = Gauss0(a, b) que no

controlan la existencia de pivotes nulos ni la posiblidad de que la matriz no sea invertible.
• escribir funciones [x, info] = solTriSup(a, b), [a, b, info] = triang(a0, b0),
[x, info] = Gauss(a, b) que controlan la existencia de pivotes nulos, que hacen intercambio de
filas cuando es necesario y que consideran la posibilidad de que la matriz sea no invertible. En estas
funciones se pueden utilizar las operaciones matriciales de Scilab.
• escribir funciones [x, info] = solTriSup1(a, b), [a, b, info] = triang1(a0, b0),
[x, info] = Gauss1(a, b), semejantes a las 3 anteriores pero que no utilizan las operaciones
matriciales de Scilab. Por ejemplo, la orden a(i,i+1:n)*x(i+1:n) se puede remplazar por
t = 0
for j = i+1:n
t = t + a(i,j)*x(j)
end

3. Escriba funciones
[a, b, info] = triangPP_R(a0, b0, nc)
[x, info] = GaussPP_R(a, b, nc)
que, respectivamente, triangulariza un sistema de ecuaciones con pivoteo parcial y resuelve un sistema

4
 de ecuaciones por el método de Gauss con pivoteo parcial. En estas funciones las operaciones se hacen
utilizando redondeo con nc cifras.
Puede ser útil escribir funciones previas análogas a las indicadas para las funciones triangR, etc.

4. Construya un ejemplo 2 × 2 o 3 × 3 o 4 × 4, donde:

• no haya diferencia (o muy pequeña) al utilizar GaussR y GaussPP_R con 15 cifras significativas.
• haya diferencia clara al utilizar GaussR y GaussPP_R con un número de cifras entre 4 y 12.

5
0.4 Método de Cholesky

1. Con un valor grande y adecuado de n (este valor de n será común para este paso y los siguientes)
construya aleatoriamente dos matrices cuadadradas. Por medio de las funciones tic() y toc() obtenga
una estimación del tiempo requerido para efectuar el producto.
Es posible que obtenga un aviso semejante a

!--error 17
rand: stack size exceeded (Use stacksize function to increase it)

Entonces mire la ayuda de stacksize y aumente el tamaño de la pila por una orden semejante a
stacksize(8000000) o stacksize(16000000) o ...

2. ¿Cuál es el número de operaciones de punto flotante del producto? ¿Qué pasa con el tiempo cuando se
duplica el valor de n? Explique.

3. Construya aleatoriamente un sistema de ecuaciones. Obtenga una estimación del tiempo requerido para
efectuar la solución del sistema por medio de la orden a\b.

4. Construya aleatoriamente una matriz simétrica definida positiva. Puede ser por ordenes semejantes a
a = rand(n,n); a = a*a’
o, preferiblemente,
a = rand(n,n); a = a+a’ + 2*n*eye(n,n)
¿Porqué?

5. Obtenga una estimación del tiempo requerido para efectuar la factorización de Cholesky por medio de la
función de Scilab chol.

6. Una de las maneras de obtener la factorización de Cholesky U T U = A es mediante las fórmulas

k−1
X
t = akk − u2ik
i=1
si t ≤ ε, entonces A no es definida positiva
√
ukk = t
k−1
X
akj − uik uij
i=1
ukj = , j = k + 1, ..., n
ukk

Escriba una función u = chol0(a) que efectúe la factorización de Cholesky. Es conveniente usar las
operaciones matriciales de Scilab. Por ejemplo, el cálculo de
k−1
X
u2ik
i=1

se puede hacer por u(1:k-1,k)’*u(1:k-1,k) o también por

norm( u(1:k-1,k) )^2

7. Obtenga una estimación del tiempo requerido por la función chol0 .

6
 8. Para resolver una sistema de ecuaciones Ax = b, con matriz definida positiva, es preferible, en lugar del
método de Gauss, utilizar el método de Cholesky, es decir

• Obtener U tal que U T U = A .

• Resolver U T y = b .
• Resolver U x = y .

¿Cual es, aproximadamente, el número de operaciones de punto flotante requeridas para el proceso
anterior?

9. Obtenga una estimación del tiempo requerido para la solución de un sistema con matriz definida positiva
mediante las tres órdenes

u = chol(a)
y = (u’)\b
x = u\y

10. Compare el tiempo anterior con el tiempo de x = a\b . Concluya.

11. Escriba una función x = solTriSup(a, b) que resuelve un sistema de ecuaciones con matriz triangular
superior. Suponga provisionalmente que los tamaños son compatibles, que la matriz sı́ es triangular
superior y que no hay elementos diagonales nulos.

12. Escriba una función x = solTriInf(A, b) que resuelve un sistema de ecuaciones con matriz triangular
inferior.

13. Escriba una función x = sol_Ut_b(U, b) que resuelve el sistema triangular inferior U T x = b donde U T
es triangular inferior, es decir, U es triangular superior. Esta función es muy parecida a la anterior, pero
la información sobre la matriz está en diferente sitio.

14. Escriba la función x = solDP(a, b) que resuelve un sistema de ecuaciones con matriz definida positiva
mediante las órdenes

u = chol(a)
y = sol_Ut_b(U, b)
x = solTris(u, y)

15. Obtenga una estimación del tiempo de la función anterior. Compare con el tiempo de x = a\b . Con-
cluya.

7
0.5 Método de Gauss-Seidel

Se desea resolver un sistema de ecuaciones lineales Ax = b. El algortimo para el método de GS se puede

escribir ası́:

datos: A, b, x0 , ε, maxit
x = x0
para k = 1, ...,maxit
nrmD← 0
para i = 1, ..., n
δi = (bi − Ai· x)/aii
xi ← xi + δi
nrmD←nrmD+|δi |
fin-para i
si nrmD ≤ ε ent x∗ ≈ x, salir
fin-para k

1. Escriba una función [x, k, info] = GS(a, b, x0, eps, maxit) que implementa el método de GS.
Los resultados serán: x aproximación de la solución (posiblemente), k el número de iteraciones completas
realizadas e info que valdrá 1 si se obtuvo una buena aproximación de la solución.

2. Considere una matriz A de la forma

d 1 0 v1 0 0 0 · · ·
 
 0 d 2 0 v2 0 0 
 
 u 1 0 d 3 0 v3 0 
 
 0 u 2 0 d 4 0 v4 
 
A=



 
 

 dn−2 0 vn−2 

 0 dn−1 0 
un−2 0 dn

En palabras, los elementos diagonales de A son d1 , d2 , ..., dn ; en la segunda superdiagonal están los
valores v1 , v2 , ..., vn−2 . en la segunda subdiagonal están los valores u1 , u2 , ..., un−2 . Los otros elementos
(entradas) de A son nulos.
Escriba una funcion A = matriz3(d, u, v) que, dados los tres arreglos d, u y v, construye explı́citamente
la matriz A.

3. Escriba una función [x, k, info] = GS3(d, u, v, b, x0, eps, maxit) que implementa el método
de GS para resolver Ax = b, donde A es una matriz con la forma descrita anteriormente. Esta función
no construye explı́citamente la matriz A. Los resultados parciales (puede agregar temporalmente disp
o printf ) y finales deben ser los mismos por cualquiera de los dos formas siguientes:

• utilizar GS3 .
• con matriz3 construir explı́citamente A y utilizar GS .

Es necesario adaptar el algortimo para tener en cuenta las especificidades de cada fila de A:

8
 para i = 1, ..., n
si i ≤ 2
δi ← · · ·
sino
si i ≤ n − 2
δi ← · · ·
sino
δi ← · · ·
fin-si
fin-si
xi ← · · ·
..
.
fin-para i

4. Utilice GS3 con

d = [ 10 11 12 13 14]’
u = [ 1 2 3]’
v = [ -3 -2 -1]’
b = [11 -16 34 -43 65]’

¿Cual es la solución?

5. Calcule, en función de n, el número de operaciones de punto flotante en cada iteración.

6. Escoja un valor grande de n y utilice GS3 con

d = 100*ones(n,1);
u = -2*ones(n-2,1);
v = 3*ones(n-2,1);
b = 101*ones(n,1);
b(1:2) = 103;
b(n-1:n) = 98;
eps = 1.0e-6

Construya explı́citamente A. Evalúe el tiempo de A\b y de GS3. Justifique. ¿Qué pasa cuando n es
más grande.

7. Dé una razón (teórica) que justifica la convergencia del método GS para los datos anteriores.

8. ¿Qué pasa con el número de iteraciones y con el tiempo si eps = 1.0e-12 ?

9. ¿Cuál es el mayor valor de n para poder resolver Ax = b, construyendo explı́citamente A y usando A\b?
¿Cuál es el mayor valor de n para poder resolver Ax = b, con la información de d, u y v, usando GS3?

9
0.6 Matriz de una cercha

Considere la siguiente cercha plana formada por triángulos rectángulos isósceles:

6 12 7 13 8

5 7 9 11
6 8 10

1 1 2 3 4 5
2 3 4

10 15 25

Los elementos (las barras) están numerados de 1 a 13. Los nudos están numerados de 1 a 8. La cercha está
apoyada en los nudos 1 y 5 y soporta las cargas (en toneladas) indicadas en los nudos 2, 3 y 4. En el nudo
1 la cercha está fijada horizontal y verticalmente y en el nudo 5 la cercha está fijada verticalmente. Sea
fi la fuerza soportada por el elemento i. En los nodos 2, 3, 4, 6, 7 y 8 se pueden plantear dos ecuaciones:
equilibrio horizontal y equilibrio vertical. Por ejemplo, en el nodo 3:

f2 + cf7 = f3 + cf9
cf7 + f8 + cf9 = 15,
√
donde c = 2/2. En el nodo 5 solamente se plantea el equilibrio horizontal. Ası́ se tienen 13 ecuaciones con
13 incógnitas.

1. Plantee y resuelva el sistema de ecuaciones.

2. ¿Cual elemento soporta una mayor fuerza y cual es esta fuerza?

3. Escriba en formato “0X” (0 para los valores nulos y X para los valores no nulos) la matriz del sistema de
ecuaciones. ¿Cual es la densidad de la matriz?

4. Considere ahora las cargas 10d1 , 10d2 y 10d3 donde d1 , d2 y d3 son los tres últimos dı́gitos no nulos de
su documento de identidad. Plantee y resuelva el sistema de ecuaciones.

5. ¿Cual elemento soporta una mayor fuerza y cual es esta fuerza?

6. Considere ahora una sola carga D = 10(d1 + d2 + d3 ) en el nudo 2. ¿Cual elemento soporta una mayor
fuerza y cual es esta fuerza?

7. Misma pregunta para una sola carga D en el nudo 3.

8. Misma pregunta para una sola carga D en el nudo 4.

9. Misma pregunta para una carga D igualmente repartida en los tres nudos (D/3 en cada nudo).

10
0.7 Método de Muller

Este método sirve para hallar raı́ces reales o complejas de polinomios reales p(x) = a0 +a1 x+a2 x2 +...+an xn .
El polinomio p se puede expresar en función de sus raı́ces:

p(x) = an (x − r1 )(x − r2 ) · · · (x − rn ).

Las raı́ces complejas, no reales, siempre vienen por parejas, es decir si r = a + ib, b 6= 0, es una raı́z, entonces
r̄ = a − ib, el conjugado de r, también es raı́z. Para las raı́ces complejas q(x) = (x − r)(x − r̄) divide a p(x).

q(x) = (x − r)(x − r̄) = (x − a − ib)(x − a + ib) = (x − a)2 + b2 = x2 − 2ax + (a2 + b2 ).

Si se halla una raı́z real r entonces q(x) = (x − r) divide a p(x).

En general, si q(x) divide a p(x), entonces existe un polinomio s(x) tal que

p(x) = q(x)s(x),
grado(p) = grado(q) + grado(s).

Entonces para sequir obteniendo las raı́ces de p(x) basta con obtener las raı́ces de s(x), polinomio más
sencillo.
En el método de la secante, dados dos valores x0 y x1 se busca la recta que pasa por los puntos (x0 , f (x0 )),
(x1 , f (x1 )); el siguiente valor, x2 , está dado por el punto donde la recta corta el eje x. En el método de
Muller, en lugar de una recta, se utiliza una parábola. Dados tres valores x0 , x1 y x2 , se construye la
parábola P (x) que pasa por los puntos (x0 , f (x0 )), (x1 , f (x1 )) y (x2 , f (x2 )); el siguiente valor, x3 , está dado
por el (un) punto tal que P (x3 ) = 0.
La parábola se puede escribir de la forma P (x) = a(x − x2 )2 + b(x − x2 ) + c. Las fórmulas que permiten
calcular a, b y c son:

d = (x0 − x1 )(x0 − x2 )(x1 − x2 ),

−(x0 − x2 )(f (x1 ) − f (x2 )) + (x1 − x2 )(f (x0 ) − f (x2 )
a= ,
d (3)
(x0 − x2 )2 (f (x1 ) − f (x2 )) − (x1 − x2 )2 (f (x0 ) − f (x2 )
b= ,
d
c = f (x2 ).

Entonces √
b2 − 4ac
−b ±
x3 − x2 =
2a
Para reducir los errores de redondeo se “racionaliza” el numerador y buscando que el denominador resultante
sea grande (en valor absoluto)

D = b2 − 4ac
√
R= D
(
b+R si |b + R| ≥ |b − R| (4)
δ=
b−R si |b + R| < |b − R|
2c
x3 = x2 − ·
δ

Si en una iteración D = b2 − 4ac < 0 es necesario utilizar, a partir de ahı́, aritmética compleja (Scilab lo
hace automáticamente). Eso hace que los siguientes valores a, b y c no sean necesariamente reales. Muy
posiblemente b2 − 4ac tampoco es real.

11
 I. MÉTODO DE MULLER PARA UNA RAÍZ
datos: p, x0 , x1 , x2 , εf , ε0 , maxit
f0 = p(x0 ), f1 = p(x1 ), f2 = p(x2 )
info= 0
para k = 1, ..., maxit
si |f2 | ≤ εf ent r = x2 , info= 1, parar
d = (x0 − x1 )(x0 − x2 )(x1 − x2 )
si |d| ≤ ε0 ent parar
calcular a, b y c según (3)
D =√ b2 − 4ac
R= D
δ1 = b + R, δ2 = b − R
si |δ1 | ≥ |δ2 | ent δ = δ1 , sino δ = δ2
si |δ| ≤ ε0 ent parar
x3 = x2 − 2c/δ
x0 = x1 , x1 = x2 , x2 = x3 , f0 = f1 , f1 = f2
f2 = p(x2 )
fin-para k

Si el algoritmo anterior acaba normalmente, info valdrá 1 y r será una raı́z, real o compleja. En la anterior
descripción del algoritmo, |t| indica valor absoluto si t es real, e indica módulo si t es complejo.

1. Averigüe cómo efectuar división de polinomios en Scilab.

2. Escriba una función [r, info] = Muller1(coef, x0, x1, x2, epsf, eps0, maxit) que obtiene (si
al final info vale 1) una raı́z de un polinomio cuyos coeficientes reales están, en orden creciente de la
potencia, en el vector coef.

3. Escriba una función [rr, info] = Muller(coef, x0, epsf, eps0, maxit) que obtiene (si al final
info vale 1) todas las raı́ces de un polinomio cuyos coeficientes reales están, en orden creciente de la
potencia, en el vector coef. Las raı́ces quedarán en el vector rr . Esta función utiliza varias veces
Muller1 .

II. MÉTODO DE MULLER

datos: p, x0 , εf , ε0 , maxit
r = x0 , h = 0.5
mientras grado(p) ≥ 3
x0 = r, x1 = x0 + h, x2 = x1 + h
(r, inf o) = M uller1(p, xo , x1 , x2 , εf , ε0 , maxit)
si inf o = 0, ent parar
si |imag(r)| ≤ ε0 ent q(x) = (x − r)
sino q(x) = (x − r)(x − r̄)
p(x) = p(x)/q(x)
fin-mientras
calcular raı́ces de p (de grado no superior a 2)

4. Sean k1 , k2 , ..., k7 los primeros siete dı́gitos del número de su documento de identidad. Halle las seis
raı́ces del polinomio cuyos coeficientes, en orden creciente de la potencia, son k7 + 1, k6 , ..., k1 . Muestre
los resultados intermedios para la primera utilización de Muller1 .

12
0.8 Interpolación polinomial

Dados n puntos (x1 , y1 ), (x2 , y2 ), ..., (xn , yn ), con xi 6= xj si i 6= j, entonces existe un único polinomio p de
grado menor o igual a n − 1 tal que
p(xi ) = yi , i = 1, ..., n.

El polinomio p(x) = a0 + a1 x + a2 x2 + · · · an−1 xn−1 se puede obtener de dos maneras: resolviendo un sistema
de ecuaciones y utilizando polinomios de Lagrange.

1. El sistema de ecuaciones n × n que hay que resolver es

x21 x31 x1n−1

    
1 x1 ··· a0 y1
n−1  
1 x2
 x22 x32 ··· x2   a1   y2 
 

1 x3
 x23 x33 ··· x3n−1 
  a2  =  y3 
   
. ..   ..   .. 
 .. .  .   . 
1 xn x2n x3n · · · xnn−1 an−1 yn

Escriba una función p = polInterpSE(x, y) que obtiene el polinomio de interpolación resolviendo el

sistema de ecuaciones. Las coordenadas de los puntos están en los vectores columna x y y. Como x es
un vector columna, entonces las columnas de la matriz de coeficientes del sistema se pueden construir en
Scilab por órdenes análogas a
x.^k

2. Hay n polinomios de Lagrange, L1 (x), L2 (x), ..., Ln (x),

n
Y
(x − xi )
i=1,i6=k
Lk (x) = n
Y
(xk − xi )
i=1,i6=k

Estos polinomios se pueden construir de manera iterativa

Lk (x) ← 1
d←1
para i = 1, ..., n
si i 6= k
d ← d (xk − xi )
Lk (x) ← Lk (x) (x − xi )
fin-si
fin-para
Lk (x) ← Lk (x)/d

Escriba una función Lk = poliLagr(x, k) que obtiene el polinomio de Lagrange Lk (x). Los valores
xi están en el arreglo x. La siguiente tabla muestra algunos ejemplos de construcción de polinomios en
Scilab.

polinomio Scilab
0 poly([0], ’x’, ’c’)
1 poly([1], ’x’, ’c’)
x−b poly([-b 1], ’x’, ’c’)
x−b poly([b], ’x’, ’r’)

13
3. La fórmula para obtener p utilizando polinomios de Lagrange es
n
X
p(x) = yk Lk (x).
k=1

En seudolenguaje,

p(x) ← 0
para k = 1, ..., n
construir Lk (x)
p(x) ← p(x) + yk Lk (x)
fin-para

Escriba una función p = polInterpPL(x, y) que obtiene el polinomio de interpolación utilizando poli-
nomios de Lagrange. Las coordenadas de los puntos están en los vectores columna x y y.

4. Escriba una función que calcula el factorial de un entero no negativo.

5. Sean x , y dos vectores columna con las coordenadas de n puntos: ( 0, (−1)n−1 (n − 1)! ), (1, 0), (2, 0),
(3, 0), ..., (n − 1, 0) . Considere n = 5 y construya explı́citamente los vectores columna x, y. Obtenga
p, el polinomio de interpolación, usando polInterpSE. Construya un vector columna yse con los valores
del polinomio p evaluado en los valores del vector columna x. Teóricamente los dos vectores, y y yse son
iguales. Obtenga la norma de y-yse.

6. Obtenga p, el polinomio de interpolación, usando polInterpPL. Construya un vector columna ypl con
los valores del polinomio evaluado en los valores del vector columna x. Obtenga la norma de y-ypl.

7. Repita los dos pasos anteriores para n = 3, 4, ..., 20. Construya una tabla de tres columnas, en la primera
los valores de n, en la segunda columna los valores norm(y-yse) y en la tercera columna los valores
norm(y-ypl). Concluya.

14
0.9 Trazadores cúbicos (splines)

Dados n puntos (x1 , y1 ), (x2 , y2 ), ... (xn , yn ), con

x1 < x2 < · · · < xn ,

el trazador cúbico, spline, o función polinomial cúbica por trozos, se define ası́:


 S1 (x) si x ∈ [x1 , x2 ]

S2 (x)

si x ∈ [x2 , x3 ]
S(x) = . (5)
..



Sn−1 (x) si x ∈ [xn−1 , xn ]


En cada uno de los n − 1 intervalos, Si (x) es un polinomio cúbico.

Si (x) = ai (x − xi )3 + bi (x − xi )2 + ci (x − xi ) + di , i = 1, ..., n − 1. (6)

Conocer S(x) quiere decir conocer 4(n − 1) coeficientes: ai , bi , ci , di , para i = 1, ..., n − 1.

Se requiere que S(x) pase por los puntos, y que sea doblemente derivable.

S(xi ) = yi , i = 1, ..., n
Si (xi+1 ) = Si+1 (xi+1 ), i = 1, ..., n − 2
Si′ (xi+1 ) = Si+1
′
(xi+1 ), i = 1, ..., n − 2
Si′′ (xi+1 ) = Si+1
′′
(xi+1 ), i = 1, ..., n − 2

Supongamos además que S tiene condiciones de “frontera libre”:

S1′′ (x1 ) = 0 ,
′′
Sn−1 (xn ) = 0 .

Sea hj = xj+1 − xj , el tamaño del intervalo [xj , xj+1 ]. Las condiciones anteriores se convierten en:

Si (xi ) = di = yi i = 1, ..., n − 1,
Sn−1 (xn ) = an−1 h3n−1 + bn−1 h2n−1 + cn−1 hn−1 + dn−1 = yn
ai h3i + bi h2i + ci hi + di = di+1 i = 1, ..., n − 2,
3ai h2i + 2bi hi + ci = ci+1 i = 1, ..., n − 2,
6ai hi + 2bi = 2bi+1 i = 1, ..., n − 2.

Para resolver estas ecuaciones:

• di = yi , i = 1, ..., n − 1.

• hi = xi+1 − xi , i = 1, ..., n − 1.

15
 • resolver Ab = ζ, sistema tridiagonal n × n. Sea α ∈ Rn la diagonal, β ∈ Rn−1 la superdiagonal,
γ ∈ Rn−1 la subdiagonal,

α = (1, 2(h1 + h2 ), 2(h2 + h3 ), ..., 2(hn−2 + hn−1 ), 1) ,

β = (0, h2 , h3 , ..., hn−1 )
γ = (h1 , h2 , ..., hn−2 , 0)
ζ1 = 0
ζi = 3(yi−1 − yi )/hi−1 + 3(−yi + yi+1 )/hi , i = 2, ..., n − 1,
ζn = 0

Temporalmente, b ∈ Rn .
1 hi
• ci = (yi+1 − yi ) − (2bi + bi+1 ), i = 1, ..., n − 1.
hi 3
bi+1 − bi
• ai = , i = 1, ..., n − 1.
3hi
• suprimir bn .

1. Escriba una función A = matrTrid(d, supd, subd) que construye una matriz tridiagonal. Los parámetros
son vectores columna de n, n − 1 y n − 1 elementos respectivamente. Corresponden a la diagonal, super-
diagonal y subdiagonal.

2. Escriba una función x = solTrid(d, supd, subd, b) que resuelve el sistema tridiagonal Ax = b. Puede
usar un método eficiente para matrices tridiagonales o simplemente construir explı́citamente la matriz
tridiagonal.

3. Sean k1 k2 los dos últimos dı́gitos del número de su cédula. Considere los puntos
(1, 1), (3, 1), (4, (5 + k1 + k2 )/10), (5, 1), (7, 1).
Con ellos construya dos vectores columna x , y .

4. Escriba un función [a, b, c, d] = splineCub(x, y) que calcula los coeficientes del trazador cúbico
S(x). Utilı́cela para los vectores x , y .

5. Escriba una función st = evalSpline(x, a, b, c, d, t) tal que dado los valores xi en x, los vectores
a, b, c y d con los coeficientes de S(x) y un número t, calcula S(t).


1
 si t ≤ x2
j = n − 1 si t ≥ xn−1

i si xi ≤ t ≤ xi+1


S(t) = Sj (t)

6. Grafique los puntos (xi , yi ) .

7. Construya el vector columna T = ( x(1):0.1:x(n) )’ . Construya un vector columna S con los

valores de S en los puntos T(i) .

8. Grafique los puntos cuyas coordenadas están en S y T.

16
0.10 Error local y global del método del
trapecio y de Simpson
1. Escriba una función [a, b] = rectaMinC(X, Y) que obtiene los dos coeficientes de la recta de aproxi-
mación por mı́nimos cuadrados y = ax + b para los puntos cuyas coordenadas están en los dos vectores
columna X y Y .
m = size(X,1)
A = [ones(m,1) X]
aa = (A’*A)\(A’*Y)
b = aa(1)
a = aa(2)
2. Escriba una función s = trapecio(f, a, b, n) que obtiene una aproximación de la integral
definida por medio de la fórmula del trapecio

n−1
!
b
f (a) + f (b) X
Z
f (x) dx ≈ h + f (xi )
a 2
i=1
b−a
h=
n
xi = a + ih

Las órdenes importantes son semejantes a:

h = (b-a)/n
s = ( f(a) + f(b) )/2
for i=1:n-1
s = s + f( a+i*h)
end
s = s*h
3. Escriba una función s = Simpson(f, a, b, n) que obtiene una aproximación de la integral definida
por medio de la fórmula de Simpson (con n par)

 
b n−1 n−2
h
Z X X
f (x) dx ≈ f (a) + f (b) + 4 f (xi ) + 2 f (xi )
a 3
i=1,3,... i=2,4,...
b−a
h=
n
xi = a + ih

4. Escoja una función (diferente de un polinomio de grado inferior a 7) de la cual conoce la (una)
antiderivada, es decir, para la cual se puede obtener fácilmente el valor exacto de una integral definida.
Por ejemplo f (x) = ex . Con a = 0, h = 0.2, b = a + h, n = 1, obtenga el valor aproximado por la
Rb
fórmula del trapecio, el valor exacto de a f (x) dx y e el valor absoluto del error cometido. Estos
valores de h y e deben ser almacenados en arreglos H y E.
5. Repita el paso anterior para h = 0.21, 0.22, ..., 0.25 (de nuevo b = a + h, n = 1).
6. Se puede suponer que el error local es aproximadamente

eloc ≈ chk

17
 entonces

log(eloc ) ≈ log(c) + k log(h) .

Es decir, se tendrı́a la ecuación de una recta. Graficar log(E) contra log(H). ¿Corresponde aproxi-
madamente a una recta?
7. Por medio de la orden [k, logc] = rectaMinC(log(H), log(E)) obtenga el valor de k. ¿Corres-
ponde aproximadamente a lo previsto teóricamente para el error local del método del trapecio?
8. Ahora se va a “mirar” el orden del error global. Sea a = 0, b = 0.6, n = 2, 3, 4, 5, 6, h = (b − a)/n.
Rb
Para esos valores de h, obtenga por la fórmula del trapecio, el valor exacto de a f (x) dx y e el valor
absoluto del error cometido. Estos valores de h y e deben ser almacenados en arreglos H y E.
9. Se puede suponer que el error global es aproximadamente

eglob ≈ dhm

entonces

log(eglob ) ≈ log(d) + m log(h) .

Es decir, se tendrı́a la ecuación de una recta. Graficar log(E) contra log(H). ¿Corresponde aproxi-
madamente a una recta?
10. Por medio de la orden [m, logd] = rectaMinC(log(H), log(E)) obtenga el valor de m. ¿Corres-
ponde aproximadamente a lo previsto teóricamente para el error global para el método del trapecio?
11. Se va a repetir el proceso anterior para la fórmula de Simpson. Para el error local de la fórmula
de Simpson, efectuar pasos análogos a 4, 5, 6 y 7, con a = 0, h = 0.2, 0.21, 0.22, 0.23, 0.24, 0.25,
b = a + h, n = 2 .
12. Para el error global de la fórmula de Simpson, efectuar pasos análogos a 8, 9 y 10, con a = 0, b = 0.6,
n = 2, 4, 6, 8, 10, h = (b − a)/n.

18
0.11 Comparación de métodos de Simpson y Gauss-Legendre

Se desea comparar para una función diferente de un polinomio, con antiderivada conocida (se puede
obtener el valor exacto de una integral definida), el método de Simpson y el método de cuadratura de
Gauss-Legendre, utilizando el mismo número de evaluaciones de la función.

1. Escoja una función, diferente de un polinomio, de la cual conoce la (una) antiderivada. Por ejemplo
f (x) = e−2x + cos(3x). Cada grupo debe escoger una función diferente. Escoja un intervalo
“razonable”, por ejemplo, a = 1, b = 2. Escriba una función y = f08(x) donde se calcula la función
f.
Z b
2. Escoja dos valores a < b y obtenga el valor exacto de f (x) dx
a
Z b
3. Obtenga una aproximación de f (x) dx utilizando intg .
a
Z b
4. Obtenga una aproximación de f (x) dx utilizando la fórmula de Simpson con 6 subintervalos (hay
a
7 evaluaciones de f ).
5. Escriba una función s = Simpson38(f, a, b, n) que obtiene una aproximación de la integral
definida por medio de la fórmula de Simpson 3/8 con n múltiplo de 3.
Z x3
3
f (x) dx ≈ h(y0 + 3y1 + 3y2 + y3 ).
x0 8

Para el intervalo [a, b], con n (múltiplo de 3) subintervalos

b
3
Z
f (x) dx ≈ h(y0 + 3y1 + 3y2 + 2y3 + 3y4 + 3y5 + 2y6 + 3y7 + · · · + 3yn−2 + 3yn−1 + yn ).
a 8

Los pasos principales pueden ser:

h = (b − a)/n
s ← f (a) + f (b)
n−2
X
s←s+3 f (xi )
i=1,4,...
n−1
X
s←s+3 f (xi )
i=2,5,...
n−3
X
s←s+2 f (xi )
i=3,6,...

s ← 3hs/8

donde xi = a + ih. O también,

h = (b − a)/n
s ← f (a) + f (b)
para i = 1, ..., n − 1
y ← f (a + ih)
si 3|i ent s ← s + 2y
sino s ← s + 3y
fin-para
s ← 3hs/8

19
 La notación p|q significa “p divide a q”. En Scilab se puede escribir if modulo(q, p) == 0
Z b
6. Obtenga una aproximación de f (x) dx mediante
a
h = (b-a)/5,
Simpson(f08, a, a+2*h, 2) + Simpson38(f08, a+2*h, b, 3)
(se utilizan 6 evaluaciones de f ).
7. En el método de cuadratura de Gauss-Legendre de orden m el cálculo numérico de la integral se hace
mediante la siguiente aproximación:
Z 1
f (x) dx ≈ wm1 f (xm1 ) + wm2 f (xm2 ) + · · · wmm f (xmm ).
−1

Los valores wmj , los pesos o ponderaciones, y los valores xmj , las abcisas, son fijos y conocidos.

Cuadratura de Gauss-Legendre
m wmj xmj
1 2 0
2 1 ±0.577350269189626
3 0.888888888888889 0
0.555555555555556 ±0.774596669241483
4 0.339981043584856 ±0.652145154862546
0.861136311594053 ±0.347854845137454
5 0.568888888888889 0
0.478628670499366 ±0.538469310105683
0.236926885056189 ±0.906179845938664
6 0.467913934572691 ±0.238619186083197
0.360761573048139 ±0.661209386466265
0.171324492379170 ±0.932469514203152

En Abramowitz y Stegun hay una tabla muy completa. Para integrar en otro intervalo es necesario
hacer un cambio de variable,
b m  
b−aX b−a a+b
Z
f (x) dx ≈ wmj f xmj + .
a 2 2 2
j=1

Rb
Escriba una función s = GaussLeg1(f, a, b, m) que calcula de manera aproximada a f (x) dx por
el método de cuadratura de Gauss-Legendre de orden m. En esta función se definen y se utilizan dos
matrices X y W. La matriz X es triangular inferior de tamaño 6 × 6: en las primeras m posiciones de
la fila m están las m abcisas y, de manera análoga, la matriz W tiene los pesos correspondientes. La
parte principal de esta función puede ser semejante a
X = [ 0 0 0 0 0 0;
//
-0.577350269189626 0.577350269189626 0 0 0 0
//
0 -0.774596669241483 0.774596669241483 0 0 0
//
-0.339981043584856 0.339981043584856 ..
-0.861136311594053 0.861136311594053 0 0
//
0 -0.538469310105683 0.538469310105683 ..
-0.906179845938664 0.906179845938664 0
//

20
 -0.238619186083197 0.238619186083197 ..
-0.661209386466265 0.661209386466265 ..
-0.932469514203152 0.932469514203152
]
W = [2 0 0 0 0 0
1 1 0 0 0 0
0.88888888888889 0.555555555555556 0.555555555555556 0 0 0
//
0.652145154862546 0.652145154862546 ..
0.347854845137454 0.347854845137454 0 0
//
0.568888888888889 0.478628670499366 0.478628670499366 ..
0.236926885056189 0.236926885056189 0
//
0.467913934572691 0.467913934572691 ..
0.360761573048139 0.360761573048139 ..
0.171324492379170 0.171324492379170
]

u = (b-a)/2
v = (a+b)/2
s = 0
for j=1:m
tj = u*X(m,j) + v
s = s + W(m,j)*f(tj)
end
s = u*s
8. Si el intervalo [a, b] es muy grande, éste se puede dividir en n subintervalos y en cada uno de ellos utilizar
cuadratura de Gauss-Legendre de orden m. Escriba una función s = GaussLeg(f, a, b, m, n) que
Rb
calcula de manera aproximada a f (x) dx por el método de cuadratura de Gauss-Legendre de orden
m en cada uno de los n subintervalos. La parte principal de esta función puede ser semejante a
s = 0
h = (b-a)/n
for i=1:n
s = s + GaussLeg1(f, a+(i-1)*h, a+i*h, m)
end
9. Compare ahora los errores obtenidos, utilizando 6 evaluaciones de la función, en los sguientes casos
a) h = (b-a)/5, Simpson(f08, a, a+2*h, 2) + Simpson38(f08, a+2*h, b, 3)
b) promedio del error con Simpson n = 4 y del error con n = 6.
c) GL, m = 1, n = 6
d) GL, m = 2, n = 3
e) GL, m = 3, n = 2
f) GL, m = 6, n = 1
Cambie de intervalo. Cambie de función. Observe y concluya.

21
0.12 Método de Newton en Rn
y cálculo numérico de la matriz jacobiana
1. Considere una función f : Rn → Rn , por ejemplo,

f (x1 , x2 ) = ( 2x21 + 3x1 x2 + 4x2 − 38 − k1 /10, −5x1 x2 + 6x1 − 6x2 + 36 + k2 /10 )

Escriba una función fx = f09(x) que, para un vector columna x, construye el vector columna fx con
las dos componentes de f (x).
2. Encuentre, usando fsolve, un vector x tal que f (x) = 0 para la función definida en f09.
3. Dada una función f : Rn → Rn , se denotará por f ′ (x) la matriz jacobiana de f evaluada en x. En
esta matriz n × n, la fila i está formada por las n derivadas parciales de fi (la i-ésima componente de
f ).
Escriba una función J = jac09(x) que para un vector columna x, construye la matriz jacobiana J de
la función definida en f09.
4. En el método de Newton en Rn (para hallar un x tal que f (x) = 0), en cada iteración hay dos pasos
fundamentales

resolver f ′ (xk ) dk = −f (xk )

xk+1 = xk + dk .

Estos pasos presuponen que para caulquier xk se puede obtener f (xk ) y f ′ (xk ). El método termina
satisfactoriamente si

||f (xk )|| ≤ ε .

Las terminaciones no deseadas, pero posibles, son
• demasiadas iteraciones sin encontrar un x adecuado.
• f ′ (xk ) no es invertible (su determinante es nulo o casi nulo).
Escriba una función [x, info, k] = Newton1(f, jac, x0, eps, maxit) que implementa el método
de Newton. El resultado info valdrá 1 si se obtuvo un x adecuado, 2 si en maxit iteraciones no se
pudo obtener una solución adecuada y valdrá 0 si en alguna iteración la matriz jacobiana es singular
o casi singular. En la función f se calcula f (x), en la función jac se calcula f ′ (x), el vector columna
x0 es la aproximación inicial. El número de iteraciones realizadas estará en k.
5. Verifique la función anterior mediante [x, info, k] = Newton1(f09, jac09, [1; 1], 1.0e-4, 20)
y con otros ejemplos.
6. Averigüe mediante qué función de Scilab puede obtener numéricamente la matriz jacobiana de una
función.
Escriba una función [x, info, k] = Newton2(f, x0, eps, maxit) que implementa el método de
Newton. Aquı́ no es necesario tener una función externa que evalúa la matriz jacobiana, basta con
tener f .
7. Para una función ϕ : R → R la derivada se puede aproximar por
ϕ(x + h) − ϕ(x) ϕ(x + h) − ϕ(x − h)
ϕ′ (x) ≈ y también ϕ′ (x) ≈
h 2h
La derivada parcial se puede aproximar por

∂fi fi (x1 , x2 , ..., xj + h, ..., xn−1 , xn ) − fi (x1 , ..., xn )

(x) ≈
∂xj h

22
 con h cercano a 0. Una mejor aproximación es

∂fi fi (x1 , x2 , ..., xj + h, ..., xn−1 , xn ) − fi (x1 , x2 , ..., xj − h, ..., xn−1 , xn )

(x) ≈ ·
∂xj 2h

Escriba una función Dij = der_fi_xj(f, i, j, x) que calcula la anterior aproximación.

8. Escriba una función J = jacobiana3(f, x) que calcula aproximadamente la matriz jacobiana.
9. Escriba una función [x, info, k] = Newton3(f, x0, eps, maxit) que implementa el método de
Newton aproximando la matriz jacobiana mediante la función anterior.
10. Escriba una función Dj = der_f_xj(f, j, x) que construye un vector columna con las n derivadas
parciales con respecto a xj .
11. Escriba una función J = jacobiana4(f, x) que calcula aproximadamente la matriz jacobiana us-
ando al función anterior.
12. Escriba una función [x, info, k] = Newton4(f, x0, eps, maxit) que implementa el método de
Newton aproximando la matriz jacobiana mediante la función anterior.
13. Con los mismos parámetros de entrada, x0 y eps, observe en número de iteraciones para cada una de
las 4 implementaciones.
14. Una manera más precisa para tratar de medir la eficiencia, es mediante el número de llamados a la
función f . Obtenga este número para cada una de las implementaciones (puede utilizar global).
Concluya.

23
0.13 Solución de ecuaciones e integración numérica

Considere la función
2
f (t) = e−αt .

Encuentre x tal que Z x

f (t) dt = β ,
0.1
donde
1 + k1
α= ,
10
1 + k2
β= ,
10
k1 = penúltimo dı́gito del número de cédula,
k2 = último dı́gito del número de cédula.

24
0.14 Derivación numérica, integración numérica y solución de ecua-
ciones

Sea f : R → R una función derivable, a un real, L un real positivo. Se desea encontrar b tal que la
longitud de la curva y = f (x) entre los puntos (a, f (a)) y (b, f (b)) sea igual a L, es decir, se desea
encontrar b solución de la ecuación

Z bp
ψ(b) = 1 + (f ′ (x))2 dx − L = 0. (7)
a

1. Escriba las funciones

• y = f0506(x) que calcula f (x) = ex .
• y = f0507(x) que calcula f (x) = 3 + x.
• y = f0508(x) que calcula f (x) = 12.3456.
√
• y = f0509(x) que calcula f (x) = 9 − x2 .
2. Escriba la función Df = derivada(f, x) que calcula numéricamente una aproximación de f ′ (x) por
medio de la fórmula
f (x + h) − f (x − h)
f ′ (x) ≈
2h
usando h = 0.01 . Verifique que produce buenas aproximaciones con las funciones f0506, ... para
diferentes valores de x.
p
3. Escriba la función y = g(f, x) que calcula 1 + (f ′ (x))2 . Para el cálculo de f ′ (x) use derivada.
4. Escriba la función longi = longiSimpson(f, a, b, n) que calcula numéricamente una aproximación
de Z bp
1 + (f ′ (x))2 dx
a
por el método de Simpson utilizando n (par) subintervalos. La fórmula de Simpson para aproximar
Rb
numéricamente a ϕ(x) dx con n par es

 
b
h
Z X X
ϕ(x) dx ≈ ϕ(a) + ϕ(b) + 4 ϕ(xi ) + 2 ϕ(xi ) (8)
a 3
i=1,3,...,n−1 i=2,4,...,n−2

Esta función llama varias veces la función g. Si en el caso (8) se calcula ϕ(a), ϕ(b), ..., aquı́ es necesario
g(f, a), g(f, b), ...
5. Obtenga, utilizando razones analı́ticas o geométricas, la longitud exacta de la curva y = f (x) entre
los puntos (a, f (a)) y (b, f (b)), en los siguientes casos
• f (x) = 3 + x, a = 1, b = 2.
• f (x) = 12.3456, a = 1, b = 11.
√ √ √
• f (x) = 9 − x2 , a = −3 2/2, b = 3 2/2.
Compare con los resultados numéricos obtenidos con longiSimpson.
√
6. Utilice longiSimpson para f (x) = 9 − x2 con a = −3 y b = 3, con a = −2.99 y b = 2.99 y con
a = −2.95 y b = 2.95. ¿Qué resultados arroja longiSimpson? ¿porqué?
7. Escriba la función [b, info] = secanteLongi(f, a, L, b0, eps, maxit, nSimpson) que resuelve
(7), la ecuación inicial ψ(b) = 0, por el método de la secante. Esta función llama varias veces la función

25
 longiSimpson. La fórmula general del método de la secante para resolver la ecuación ϕ(x) = 0, par-
tiendo de x0 y x1 es

ϕ(xk )(xk − xk−1 )

xk+1 = xk − , k = 1, 2, ...
ϕ(xk ) − ϕ(xk−1 )

En un programa no es necesario, en cada momento, tener todos los xk . Basta con tener x0 , x1 y x2 .

ϕ(x1 )(x1 − x0 )
x2 = x1 − . (9)
ϕ(x1 ) − ϕ(x0 )

Obviamente es necesario actualizarlos en cada iteración. Es decir, el x0 de una iteración era el x1 de

la iteración anterior, ϕ(x0 ) corresponde al ϕ(x1 ) de la iteración anterior. Si para la aplicación de la
fórmula (9) es necesario calcular ϕ(x1 ), en la función secanteLongi la variable es b y se requiere
longiSimpson(f, a, b1, n) - L
b0 , la aproximación inicial de b, es b0. En la función b1 se puede construir simplemente por b0 + h.
El parámetro eps sirve para detener el proceso iterativo cuando |ψ(b)| ≤eps. El valor maxit indica el
número máximo de iteraciones. El último parámetro es el número de subintervalos que se usará en el
métodode Simpson.
8. Sean k1 k2 k3 k4 los últimos 4 dı́gitos de su documento de identidad. Escriba la función y = fmia(x)
que calcula la función f (x) = (k1 + 1)e(k2 +1)x .
9. Obtenga b para la función anterior con a = k3 y L = k4 + 1.

26
0.15 Solución de ecuaciones, derivación numérica y optimización

Sea f : R → R derivable. Bajo condiciones adecuadas, una manera de hallar un minimizador consiste en
obtener una solución de f ′ (x) = 0. Para esto se puede utilizar el método de la secante o el método de
Newton para f ′ . Para evaluar la derivada, esta se puede aproximar numéricamente.

1. Busque una función que tenga por lo menos una raı́z real y escrı́bala en Scilab en una función llamada
por ejemplo y = f_39( x ).
2. La fórmula del método de la secante para resolver la ecuación g(x) = 0, conocidos xk−1 y xk , es

(xk − xk−1 )g(xk )

xk+1 = xk −
g(xk ) − f (xk−1 )
Para hacer un programa, en cada iteración, solo se requieren los dos últimos x, entonces se utiliza
simplemente

(x1 − x0 )g(x1 )
x2 = x1 −
g(x1 ) − f (x0 )
Obviamente, el x0 de una iteración es el x1 de la anterior y el x1 de una iteración es el x2 de la anterior.
El proceso acaba satisfactoriamente si |g(xk )| ≤ ε. Las terminaciones no deseadas, pero posibles son:
demasiadas iteraciones y un denominador nulo o casi nulo.
Escriba una función [x, info] = secante(g, x0, eps, maxit) que realiza el método de la secante.
El resultado info valdrá
0: si hay un denominador nulo o casi nulo.
1: se obtuvo un solución con la precisión deseada.
2: muchas iteraciones sin obtener un resultado aceptable.
Para hacer únicamente una evaluación de f por iteración el esquema de la función puede ser semejante
a la siguiente porción de código:

x1 = x0 + 0.1
g0 = g(x0)
g1 = g(x1)
for ...

denominador = g1 - g0
if abs(denominador) ..

x2 = ... /denominador
g2 = g(x2)
if abs(g2) <= eps

x0 = x1
x1 = x2
g0 = g1
g1 = g2

3. Verifique su funcionamiento por lo menos con dos funciones.

4. De un ejemplo de una función que cumpla cada una de las siguientes condiciones y escrı́bala en Scilab
• Tiene un minimizador.
• Tiene un maximixador.
• Tiene por lo menos un minimizador y un maximizador.

27
 • No tiene minimizador ni maximizador.
5. Considere las siguientes aproximaciones:

ϕ(x + h) − ϕ(x − h) ϕ(x + h) − 2ϕ(x) + ϕ(x − h)

ϕ′ (x) ≈ , ϕ′′ (x) ≈ ·
2h h2
Sea f : R → R derivable. Si x∗ es un minimizador (un punto de mı́nimo), entonces f ′ (x∗ ) = 0.
Escriba una función [xmin, info] = minSecante(f, x0, eps, maxit), adaptación de la función
secante, para hallar un minimizador de f mediante la obtención de una raı́z de f ′ (x). La derivada
de f se obtendrá numéricamente, es decir, se hará algo semejante a
h = 0.01
h2 = 2*h
g0 = ( f(x0+h) - f(x0-h) )/h2
El resultado info valdrá
0: hay un denominador nulo o casi nulo.
1: se obtuvo un minimizador
2: se obtuvo un maximizador
3: se obuvio un punto crı́tico
4: muchas iteraciones sin obtener un resultado aceptable.

¿Dentro de minSecante, cómo saber si el punto obtenido tal que f ′ (x) ≈ 0 corresponde a un mini-
mizador?
6. Utilice la función minSecante con las funciones de la parte 4.
7. Escriba una función [x, info] = Newton(g, der_g, x0, eps, maxit) que implementa el método
de Newton para resolver f (x) = 0. En la función g está definida g y en der_g está definida g ′ .
Recuérdese que al fórmula del método de Newton para resolver g(x) = 0 es

g(xk )
xk+1 = xk −
g ′ (xk )
8. Escriba una función [xmin, info] = minNewton(f, x0, eps, maxit) que busca un minimizador de
f buscando un raı́z de f ′ (x) = 0, mediante el método de Newton. Las derivadas f ′ y f ′′ se aproximarán
numéricamente.
9. Utilice la función minNewton con las funciones de la parte 4.

28
0.16 Control de paso: método de Fehlberg

Un método muy popular es el de Runge-Kutta-Fehlberg, construido a partir de un método de RK de

orden 5 (el método A) y de un método de RK de orden 6 (el método B). Una de sus ventajas está dada
por el siguiente hecho: los valores K1 , K2 , K3 , K4 y K5 son los mismos para los dos métodos.
En la descripción del algoritmo usaremos la siguiente notación de Matlab y de Scilab: u = [u; t]

MÉTODO RUNGE-KUTTA-FEHLBERG
datos: x0 , y0 , b, h0 , ε, hmin
x = x0 , y = y0 , h = h0
X = [x0 ], Y = [y0 ]
mientras x < b
h = min{h, b − x}
hbien = 0
mientras hbien = 0
calcular yA , yB
e = |yA − yB |/h
si e ≤ ε
x = x + h, y = yB
bienh = 1
X = [X; x], Y = [Y ; y]
fin-si
C0 = 0.84(ε/e)1/4
C = max{C0 , 0.1}, C = min{C, 4}
h = Ch
si h < hmin ent parar
fin-mientras
fin-mientras

Recuérdese que siempre K1 = hf (xi , yi ). La siguiente tabla tiene los coeficientes para el cálculo de los
Kj , donde
Kj = hf (xi + αj h, yi + βj1 K1 + βj2 K2 + · · · + βj,j−1 Kj−1 )

j αj βj1 βj2 ...

1 1
2
4 4
3 3 9
3
8 32 32
12 1932 7200 7296
4 −
13 2197 2197 2197
439 3680 845
5 1 −8 −
216 513 4104
1 8 3544 1859 11
6 − 2 − −
2 27 2565 4104 40

25 1408 2197 1
yA = yi + K1 + 0K2 + K3 + K4 − K5
216 2565 4104 5
16 6656 28561 9 2
yB = yi + K1 + 0K2 + K3 + K4 − K5 + K6
135 12825 56430 50 55

29
Los errores locales son respectivamente O(h5 ) y O(h6 ). La aproximación del error está dada por

|yA − yB |
|error| ≈ e = . (10)
h
El coeficiente para la modificación del valor de h está dado por:
 ε 1/4
C0 = 0.84 , C = min{C0 , 4}, C = max{C, 0.1}. (11)
e
Las fórmulas anteriores buscan que C no sea muy grande ni muy pequeño. Más especı́ficamente, C debe
estar en el intervalo [0.1, 4].
La salida no deseada del algoritmo anterior se produce cuando h se vuelve demasiado pequeño. Esto se
produce en problemas muy difı́ciles cuando, para mantener el posible error dentro de lo establecido, ha
sido necesario disminuir mucho el valor de h, por debajo del lı́mite deseado.
Ejemplo: Resolver, por el método RKF con control de paso, la ecuación diferencial

y ′ = 2x2 − 4x + y
y(1) = 0.7182818

en el intervalo [1, 3], con h0 = 0.5 y ε = 10−6 .

yA = −0.01834063
yB = −0.01830704
e = 0.00006717

h = 0.5 no sirve.

C0 = 0.29341805
C = 0.29341805
h = 0.14670902
yA = 0.51793321
yB = 0.51793329
e = 0.00000057

h = 0.14670902 sirve.

x = 1.14670902
y = 0.51793329
C0 = 0.96535578
C = 0.96535578
h = 0.14162640
yA = 0.30712817
yB = 0.30712821
e = 0.00000029

h = 0.14162640 sirve.

x = 1.28833543
y = 0.30712821

30
 x h ỹ(x)
1.0000000 0.1467090 0.7182818
1.1467090 0.1416264 0.5179333
1.2883354 0.1622270 0.3071282
1.4505624 0.1686867 0.0572501
1.6192491 0.1333497 -0.1946380
1.7525988 0.1329359 -0.3736279
1.8855347 0.1191306 -0.5206051
2.0046653 0.1092950 -0.6137572
2.1139603 0.1024064 -0.6566848
2.2163666 0.0971218 -0.6506243
2.3134884 0.0928111 -0.5948276
2.4062996 0.0891591 -0.4877186
2.4954587 0.0859853 -0.3273334
2.5814440 0.0831757 -0.1114979
2.6646196 0.0806534 0.1620898
2.7452730 0.0783639 0.4958158
2.8236369 0.0762674 0.8921268
2.8999043 0.0743333 1.3535162
2.9742376 0.0257624 1.8825153
3.0000000 2.0855366

1. Escriba una función [X, Y, info] = RKF(f, x0, y0, b, h0, eps, hmin) que implementa el método
RKF anterior.
2. Considere la ecuación diferencial y ′ = −2y/(1 + x2 ) + cos(5x) con la condición inicial y(a) = d, donde
a es el penúltimo dı́gito de su número de cédula y d el último dı́gito. Sea b = a + 2, h0 = 1, ε = 10−4 ,
hmin = 0.001. Encuentre una solución aproximada por el método anterior en el intervalo [a, b].
3. Muestre, para los dos primeros pasos, como en el ejemplo, los valores intermedios: yA , yB , e, C0 , C,
nuevo h, ...
4. Sea X = a:0.1:b . Obtenga la solución aproximada por Scilab en los valores de X.
5. Grafique la solución de Scilab.
6. En la misma gráfica, muestre la solución obtenida con su función RKF.

31
0.17 Método del disparo
Se utiliza para ecuaciones diferenciales de segundo orden con condiciones de frontera.
Considere la ecuación diferencial
y ′′ − 3y ′ + 2y = 4x + e3x
y(a) = ya
y(b) = yb
donde k1 k2 k3 son los 3 últimos dı́gitos de su documento de identidad, a = k1 /10, b = a + 1, ya = k2 y
yb = k3 .
Suponga temporalmente, para facilitar la exposición, que las condiciones iniciales son
y(1) = 3
y(2) = 4 .

1. Obtenga la solución exacta (analı́tica). Puede hacerlo “a mano” o utilizar un software para matemáticas
como Maxima, Derive, Mathematica, Maple, ... o acudir a un amigo(a) que sepa ecuaciones diferen-
ciales lineales de segundo orden con coeficientes constantes no homogenea.
2. Escriba esta ecuación en la forma y ′′ = f (x, y, y ′ ).
3. Mediante un cambio de variable convierta la ecuación en un sistema de dos ecuaciones diferenciales

u′1 = ϕ1 (x, u1 , u2 )
u′2 = ϕ2 (x, u1 , u2 )
u1 (a) = ya
u1 (b) = yb ,
con u1 = y y u2 = y ′ .
4. Convierta este sistema en uno con condiciones iniciales suponiendo algún valor para u2 (a). Por ejemplo
suponga que u2 (1) = 0.5, es decir se tendrı́a
u′1 = ϕ1 (x, u1 , u2 )
u′2 = ϕ2 (x, u1 , u2 )
u1 (1) = 3
u2 (1) = 0.5 .

El anterior sistema se puede escribir

u′ = ϕ(x, u) (12)
 
3
u(1) =
0.5

donde
 
u
u = u(x) = 1
u2
 ′
u
u′ = 1′
u2
 
ϕ1 (x, u1 , u2 )
ϕ(x, u) =
ϕ2 (x, u1 , u2 )

32
 5. Escriba una función der_u = f14(x, u) tal que para x, un real, y u, un vector columna 2 × 1,
implementa la función ϕ cuyo resultado es también un vector columna 2 × 1. Puede ser análoga a

function der_u = f14(x, u)

der_u = zeros(2,1)
der_u(1) = x + u(1) + u(2)
der_u(2) = x*u(1)*u(2)
endfunction

6. Aplique RK4 con h = 0.1 para resolver (12) en el intervalo [a, b]. El resultado será una matriz U de
dos filas en la que las columnas corresponden a las aproximaciones ũ(1), ũ(1.1), ũ(1.2), ..., ũ(2).
7. Si la suposición u2 (1) = 0.5 hubiera sido correcta, entonces se tendrı́a ũ1 (2) = 4. Haga otras suposi-
ciones para u2 (1) (ensayo y error) hasta que ũ1 (2) ≈ 4.
8. Sea t = u2 (a) y g = ũ1 (b) el valor obtenido después de haber aplicado RK4. Es claro que g depende
de t, es decir, g = g(t). Se desea que g(t) = u1 (b) = yb . O sea, se desea resolver

G(t) = 0

donde

G(t) = g(t) − yb .

Adapte el método de la secante para resolver G(t) = 0.

9. Ya resuelta la anterior ecuación, la última matriz U es (aproximadamente) correcta, o sea, en la
primera fila estarán los valores aproximados de u1 = y.
Haga una tabla con los valores x, y(x) exacto (obtenido según la solución en el paso 1.) y ũ1 .
10. Grafique la solución exacta y la solución obtenida numéricamente.

33