Facultad de Ingenierı́a de Minas, Geologı́a Y Civil

Escuela Profesional de Ingenierı́a Civil

MÉTODOS NUMÉRICOS

TRABAJO 1
Resolución de ejercicios propuestos del texto: ”Métodos
numéricos con Matlab” de A. Cordero , J. L. Hueso , E.
Martinez y J. R. Terregrosa

Docente:
Ing. Cristian Castro Pérez

Integrantes:
16162144 AYALA HUAMANGALLI, Nilson Eduardo
16162113 BERROCAL BAUTISTA, Jhordy Edson
16162303 CONDE CONDORI, Eddy Anderson
16162513 FLORES CURIÑAUPA, Jhersson Edgar
16162141 GOMEZ MORENO, Jhon Wilver
16162146 LUZA CASTRO, Erick Javier

Ayacucho - Perú
2019

1
”Métodos numéricos con Matlab” de A. Cordero , J. L.
Hueso , E. Martinez y J. R. Terregrosa

2
 Índice general

CAP 1 VECTORES Y FUNCIONES EN MATLAB 7

PREGUNTA 1.1 7

PREGUNTA 1.2 9

PREGUNTA 1.3 11

PREGUNTA 1.4 13

PREGUNTA 1.5 15

PREGUNTA 1.6 16

PREGUNTA 1.7 17

CAP 2 NÚMEROS COMPLEJOS Y POLINOMIOS 19

PREGUNTA 2.1 19

PREGUNTA 2.2 20

PREGUNTA 2.3 20

PREGUNTA 2.4 21

PREGUNTA 2.5 22

PREGUNTA 2.6 22

PREGUNTA 2.7 23

PREGUNTA 2.8 24

PREGUNTA 2.9 25

PREGUNTA 2.10 26

3
CAP 3 REPRESENTACIÓN GRÁFICA 27

PREGUNTA 3.1 27

PREGUNTA 3.2 28

PREGUNTA 3.3 29

PREGUNTA 3.4 31

PREGUNTA 3.5 34

PREGUNTA 3.6 37

CAP 4 LA ECUACIÓN F(x) = 0 40

PREGUNTA 4.1 40

PREGUNTA 4.2 41

PREGUNTA 4.3 43

PREGUNTA 4.4 45

PREGUNTA 4.5 47

PREGUNTA 4.6 50

PREGUNTA 4.7 54

CAP 5 APLICACIONES DE LOS SISTEMAS LINEALES 57

PREGUNTA 5.1 57

PREGUNTA 5.2 58

CAP 6 VALORES PROPIOS Y VALORES SINGULARES 62

PREGUNTA 6.1 62

PREGUNTA 6.2 64

4
CAP 7 INTERPOLACIÓN POLINÓMICA 66

PREGUNTA 7.1 66

PREGUNTA 7.2 67

PREGUNTA 7.3 69

PREGUNTA 7.4 70

PREGUNTA 7.5 71

PREGUNTA 7.6 72

PREGUNTA 7.7 73

PREGUNTA 7.8 75

CAP 8 TRAZADORES CÚBICOS 76

PREGUNTA 8.1 76

PREGUNTA 8.2 77

PREGUNTA 8.3 79

PREGUNTA 8.4 82

PREGUNTA 8.5 83

CAP 9 MODELOS LINEALES 87

PREGUNTA 9.1 87

PREGUNTA 9.2 92

PREGUNTA 9.3 96

PREGUNTA 9.4 98

PREGUNTA 9.5 99

CAP 10 POLINOMIOS ORTOGONALES 100

PREGUNTA 10.1 100

5
 PREGUNTA 10.2 103

PREGUNTA 10.3 106

PREGUNTA 10.4 109

CAP 11 APROXIMACIÓN TRIGONOMÉTRICA 110

PREGUNTA 11.1 110

PREGUNTA 11.2 111

PREGUNTA 11.3 112

CAP 12 ECUACIONES DIFERENCIALES DE PRIMER ORDEN114

PREGUNTA 12.1 114
Metodo de Euler — 114 • Metodo de Euler Modificado — 115

PREGUNTA 12.2 118

6
 VECTORES Y FUNCIONES EN MATLAB
CAPÍTULO 1
PREGUNTA 1.1

Dados dos números a < b, obtén un vector x cuyas componentes sean

los extremos de los intervalos obtenidos al dividir [a, b] en n partes
iguales. El vector x tiene n + 1 componentes a las que llamaremos
a = x0 < x1 < x2 < x3 < · · · < xn = b . Elige adecuadamente a, by n para obte-
ner representaciones gráficas de los valores en x de las siguientes funciones:

a) f (x) = sin(x) e) f (x) = cosh(x) i) f (x) = 1/x

b) f (x) = cos(x) f) f (x) = tanh(x) j) f (x) = 1/x2

c) f (x) = tan(x) g) f (x) = ln(x) k) f (x) = 1/(1 + x2 )

d) f (x) = sinh(x) h) f (x) = exp(x) l) f (x) = x/(1 + x2 )

SOLUCIÓN

Los n + 1 elementos del vector x se generan con la siguiente expresión:

xi = a + i(b − a)/n

E implementando este procedimiento en Matlab con el código siguiente. Las fun-

ciones se definen el la lı́nea 5.
1 a = −10;
2 b =10;
3 n=100;
4 x=a : ( b−a ) / n : b ;
5 f = @( x ) sin ( x ) ;
6 plot ( x , f ( x ) , ’ linewidth ’ , 1 . 5 )

7
Las gráficas generadas por cada función pedida se muestra a continuación.

8
 PREGUNTA 1.2

Estudia de forma aproximada los siguientes lı́mites de funciones sin resol-

verlos. Toma sucesiones de números que convergen al valor donde se evalúa
el limite, evalua estas sucesiones en la función para obtener su carácter con-
vergente o divergente. y en el caso de convergencia obtener una aproxima-
ción al lı́mite.
ex − 1 2x
a) lı́mx→0 e) lı́mx→∞ 5
tan(2x) x
x − arctan(x) sin(x) − 1
b) lı́mx→0 f) lı́mx→ π2
x3 2x − π
x cos(x) − sin(x) ex − sin(x) − 1
c) lı́mx→0 g) lı́mx→0
x3 sin2 (x)
log(sin(3x)) 2x − x 2
d) lı́mx→0 h) lı́mx→2
log(sin(2x)) x−2
Los lı́mites se calcularán tabulando desde valores cercanos al valor de x y con-
siderando intervalos pertinentes para cada caso. Se muestran las gráficas de las
funciones dadas para poder apreciar el comportamiento de la función.

9
a) = 0,5
b) = 0,3333
c) = −0,3333
d) = 1,0000
e) =∞
f) =0
g) = 0,5
h) = −1,22741

10
 PREGUNTA 1.3

Define las siguientes matrices en Matlab:

 1 1 1 
 2 −1 1   2 3 4
   
  2 4 8 
1 2
A =  3 2 1  B =  4 3 2  C =  1 
     
3 5
5 1 4 1 5 −5 − 16
   
1 1
 

y efectúa las siguientes operaciones:

a) Comprueba la no conmutatividad del producto de matrices calculando

A∗B y B∗A

b) Calcula C3 , utilizando las potencias de matrices elemento a elemento.

c) Define la matriz identidad de orden 3 con la orden correspondiente y

con ella calcula la inversa de A utilizando la división izquierda de ma-
trices.

d) calcula (A ∗ B)−1 de dos formas diferentes usando la teorı́a de matrices.

SOLUCIÓN

Resolviendo el problema directamente en Matlab con el código siguiente

1 A= [ 2 , − 1 , 1 ; 3 , 2 , 1 ; 5 , 1 , 4 ]
2 B=[3 ,2 ,4;4 ,3 ,2;1 ,5 , −5]
3 C= [ 1 / 2 , 1 / 4 , 1 / 8 ; 1 / 3 , 2 / 5 , 1 ; − 1 / 6 , 1 , 1 ]
4

5 AB=A*B
6 BA=B*A
7

8 C3=(C*C) *C
9

10 I=eye ( 3 )
11 Ainv=A\ I
12

13 ABinv1=inv (A*B )
14 ABinv2 = inv ( B ) * inv (A)

Se obtienen las siguientes respuestas en la ventana de comandos.

11
A =
2 -1 1
3 2 1
5 1 4

B =
3 2 4
4 3 2
1 5 -5

C =
0.5000 0.2500 0.1250
0.3333 0.4000 1.0000
-0.1667 1.0000 1.0000

AB =
3 6 1
18 17 11
23 33 2

BA =
32 5 21
27 4 15
-8 4 -14

C3 =
0.2000 0.6556 0.8266
0.2408 1.9723 2.7017
0.1646 2.5433 3.3479

I =
1 0 0
0 1 0
0 0 1

Ainv =
0.5000 0.3571 -0.2143
-0.5000 0.2143 0.0714
-0.5000 -0.5000 0.5000

ABinv1 =
-0.6351 0.0405 0.0946
0.4189 -0.0328 -0.0290
0.3919 0.0753 -0.1100

ABinv2 =
-0.6351 0.0405 0.0946
0.4189 -0.0328 -0.0290
0.3919 0.0753 -0.1100

12
 PREGUNTA 1.4

Utilizando las otras opciones de la instrucción diff de Matlab que puedes

consultar con Help, calcula las siguientes derivadas y derivadas parciales:

d 1
 
a)
dt 1 + t 2
d3  3 3 
b) t sin (t)
dt 3
d2 1
 
c) tan
dt 2 t
∂
d) (x sin(y) − y cos(x))
∂x
∂2
e) (x sin(y) − y cos(x))
∂x∂y

SOLUCIÓN

Usando el siguiente código de Matlab con las instrucciones indicadas, se obtienen

los siguientes resultados:
1 syms t x y
2 FA= ( 1 / ( 1 + t ˆ 2 ) ) ;
3 FB= t ˆ 3 * sin ( t ) ˆ 3 ;
4 FC= tan ( 1 / t ) ;
5 FD= x * sin ( x * y ) − y * cos ( x ) ;
6 % a)
7 FAt= d i f f ( FA , t , 1 )
8 % b)
9 FBt3= d i f f ( FB , t , 3 )
10 % c)
11 FCt2= d i f f ( FC , t , 2 )
12 % d)
13 FDx= d i f f (FD , x , 1 )
14 % e)
15 FDy= d i f f (FD , y , 1 ) ;
16 FDyx= d i f f ( FDy , x , 1 )

13
Resultados imprimidos:

FAt =
-(2*t)/(tˆ2 + 1)ˆ2

FBt3 =
6*sin(t)ˆ3 + 6*tˆ3*cos(t)ˆ3 -27*tˆ2*sin(t)ˆ3 +54*tˆ2*cos(t)ˆ2*sin(t) -21*tˆ3*cos(

FCt2 =
(2*(tan(1/t)ˆ2 + 1))/tˆ3 + (2*tan(1/t)*(tan(1/t)ˆ2 + 1))/tˆ4

FDx =
sin(x*y) + y*sin(x) + x*y*cos(x*y)

FDyx =
sin(x) + 2*x*cos(x*y) -xˆ2*y*sin(x*y)

14
 PREGUNTA 1.5

Los primeros términos no nulos de la serie de MacLaurin para la función

x3 x5
arco tangente son x − + . Calcula el error absoluto y relativo en las
3 5
siguientes aproximaciones de π usando el polinomio en vez de la función
arco tangente.

1 1
 
a) 4 arctan + arctan
2 3
1 1
b) 16 arctan − 4 arctan
5 239
Si p∗ es una aproximación de p se definen error absoluto como |p−p∗ | y error
|p − p∗ |
relativo como si p , 0
p

SOLUCIÓN

El código usado para resolver es el siguiente:

1 McL = @(x) x-xˆ3/3+xˆ5/5;
2 piA =4 * ( McL (1/2)+McL (1/3))
3 piB =16 * McL (1/5) -4 * McL (1/239)
4
5 % error absoluto
6 ErrorAbsolutoPiA = abs(pi -piA)
7 ErrorAbsolutoPiB = abs(pi -piB)
8 %error relativo
9 ErrorRelativoPiA = abs(pi -piA)/pi
10 ErrorRelativoPiB = abs(pi -piB)/pi

Y los resultados imprimidos, los siguientes:

piA = 3.1456
piB = 3.1416

ErrorAbsolutoPiA = 0.0040
ErrorAbsolutoPiB = 2.8376e-05

ErrorRelativoPiA = 0.0013
ErrorRelativoPiB = 9.0323e-06

15
 PREGUNTA 1.6

Calcula la norma euclidea del vector (1, −1, 2, 3, −10, 6) Hazlo de tres formas
diferentes usando los siguientes pseudocódigos que definen dos algoritmos
para determinar la norma euclidea de un vector x = (x1 , x2 , · · · , xn ):
a)

Nombre: normal
Entrada: el vector x
Salida : norma de x
Paso 1: Determinar la longitud de x, que designaremos por n
Paso 2: Hacer suma = 0
Paso 3: Para i = 1, 2, ..., n hacer suma = suma + xi2
Paso 4: Hacer norma = suma( 1/2)

b)
Nombre: norma2
Paso 1: Elevar al cuadrado elemento a elemento el vector x y nombrarlo x2
utilizando la instrucción de MATLAB correspondiente.
Paso 2: Sumar los elementos de x2 utilizando la instrucción de MATLAB
correspondiente.
Paso 3: Hacer norma igual a la raiz cuadrada del último valor obtenido en
el Paso 2 utilizando la instrucción de MATLAB correspondiente.

c)
Utiliza también la instrucción de MATLAB que calcula directamente la nor-
ma euclı́dea.

SOLUCIÓN

Los procedimientos indicados se implemantan en el siguiente código de Matlab:

1 x=[1,-1,2,3,-10,6];
2 % Primer metodo
3 n=length(x); suma =0;
4 for i=1:n
5 suma = suma +x(i)ˆ2;
6 end
7 normaA = suma ˆ(1/2)
8
9 % Segundo metodo
10 normaB=sqrt(sum(x.ˆ2))
11
12 % Tercer metodo
13 normaC = norm(x)

normaA = 12.2882 normaB = 12.2882 normaC = 12.2882

16
 PREGUNTA 1.7

El número π puede aproximarse teóricamente como la suma se la serie

(−1)n
4 ∞
P
n=0 aunque veremos que en la prática este método no es el mejor,
2n + 1
hemos visto ya en el Ejercicio 5 una forma más efectiva de cáculo de π.

a) Calcula las aproximaciones de π para n = 100, 1000, 10000 utilizando

la definición de vectores y la instrucción de MATLAB sum, utilizando el
valor de π de MATLAB calcular la exactitud de las aproximaciones.

b) Crea un archivo de funeión que tenga argumento de entrada n y

que devuelva la aproximación de A para ese número de sumandos.

c) Fijándote en los datos obtenidos en el apartado a) da una estima-

ción del n necesario para aproximar π con una precisión del orden de
10−8

SOLUCIÓN

La implementación de este ejercicio en Matlab es el siguiente código:

1 function[piA]= AproxPi(n)
2 a=0;
3 suma =0;
4 while a<n
5 suma = suma + ((-1)ˆ(a))/(2 * (a)+1);
6 a=a+1;
7 end
8 PiAprox =4 * suma
9 error=pi -4 * suma
10 end

Los resultados imprimidos son:

>> AproxPi(10)
PiAprox = 3.0418
error = 0.0998

>> AproxPi(100)
PiAprox =3.1316
error = 0.0100

>> AproxPi(1000)
PiAprox = 3.1406
error = 1.0000e-03

>> AproxPi(10000)
PiAprox = 3.1415
error = 1.0000e-04

17
>> AproxPi(100000)
PiAprox =3.1416
error =1.0000e-05

>> AproxPi(1000000)
PiAprox =3.1416
error =1.0000e-06

>> AproxPi(10000000)
PiAprox = 3.1416
error = 1.0000e-07

>> AproxPi(100000000)
PiAprox =3.1416
error =1.0000e-08

Donde se puede estimar que con n = 100000000 iteraciones se puede estimar que
el error absoluto es 10e − 8

18
NÚMEROS COMPLEJOS Y POLINOMIOS
CAPÍTULO 2
PREGUNTA 2.1

Deduce las siguientes expresiones a partir de la fórmula de Euler.

eix + e−ix eix − e−ix

cos x = sin x =
2 2i

SOLUCIÓN

De la fórmula de Euler :

eix = cos(x) + i sin(x)

Evaluando para x y −x se obtiene:

eix = cos(x) + i sin(x) (1)

e−ix = cos(−x) + i sin(−x)

e−ix = cos(x) − i sin(x) (2)

Sumando y restando (1) y (2) se obtiene:

eix + e−ix = 2 cos(x)

eix − e−ix = 2i sin(x)

despejando sin x y cos x queda:

eix + e−ix
cos(x) =
2
eix − e−ix
sin(x) =
2i

19
 PREGUNTA 2.2

Deduce de la expresión eai ebi = e(a+b)i las fórmulas del seno y coseno del
ángulo suma a + b.

SOLUCIÓN

A partir de la fórmula de Euler : eix = cos(x)+i sin(x) , para x = a y x = b se obtiene:

eai ebi = (cos a + i sin a) (cos b + i sin b)

eai ebi = cos a cos b + i sin b cos a + i sin a cos b − sin a sin b
eai ebi = (cos a cos b − sin a sin b) + i (sin b cos a + sin a cos b)
eai ebi = (cos(a + b)) + i (sin(a + b))
eai ebi = cos(a + b) + i sin(a + b)
eai ebi = e(a+b)i

PREGUNTA 2.3
√
Prueba que la expresión A cos t+B sin t puede ponerse como A2 + B2 sin(t+
φ) para determinado valor de φ

SOLUCIÓN

Sea φ tal que:

A
sin φ = √
A2 + B2
B
cos φ = √
A2 + B2
√
Si a la expresión A cos t + B sin tla multiplicamos y dividimos por A2 + B2 se ob-
tiene lo siguiente:
√ A cos t + B sin t
!
A2 + B2 √
A2 + B2
√ A cos t B sin t
!
2
= A +B √ 2 +√
A2 + B2 A2 + B2
√ A B
!
= A2 + B2 √ cos t + √ sin t
A2 + B2 A2 + B2
√
= A2 + B2 (sin φ cos t + cos φ sin t)
√
∴ A cos t + B sin t = A2 + B2 sin(t + φ)

20
 PREGUNTA 2.4

Divide el polinomio x3 − 2x − 5 por x − 2

a) manualmente, utilizando el algoritmo de división.

b) por la regla de ruffini

c) utilizando la orden deconv

Evalúa el polinomio en x = 2

a) por sustitución directa

b) usando polyval

SOLUCIÓN
Haciendo el algoritmo de la división se obtiene x2 + 2x + 2 como cociente y −1 de
residuo

Y haciendo el proceso de Ruffini, se obtiene lo mismo.

1 0 -2 -5
2 2 4 4
1 2 2 -1
Usando deconv en Matlab se obtiene lo siguiente:
u = [1 0 -2 -5];
v = [1 -2];
[q,r] = deconv(u,v)

q = 1 2 2
r = 0 0 0 -1
Evaluando el polinomio obtenido en x = 2 se obtiene:
p(x) = x2 + 2x + 2
p(2) = 22 + 2 · 2 + 2 p(2) = 10
usando polyval:
>> p=[1 2 2];
>> polyval(p,2)

ans = 10

21
 PREGUNTA 2.5

halla (x − 1)6 a partir de sus raı́ces.

Usando Matlab para resolver este ejercicio, se define al vector r que cpntiene
las raı́ces del polinomio del enunciado. Usando poly para haalr un polinomio a
partir del vector que contenga sus raı́ces, se obtiene:
>> r = [1 1 1 1 1 1];
>> p=poly(r)

p =

1 -6 15 -20 15 -6 1

Donde p muestra los coeficientes del polinomio del enunciado en su versión ex-
tendida:

p(x) = 1x6 − 6x5 + 15x4 − 20x3 + 15x2 − 6x + 1

PREGUNTA 2.6

Halla los coeficientes del polinomio (2x − 3)5 y represéntalo gráficamente

en el intervalo [1, 2].

Para hallar los coeficientes del polinomio enunciado, se procede con el mismo
procedimiento del ejercicio anterior, donde las raı́ces son 3/2 con multiplicidad 5
>> r = [3/2 3/2 3/2 3/2 3/2];
p=poly(r)

p =
1.0000 -7.5000 22.5000 -33.7500 25.3125 -7.5938

22
 PREGUNTA 2.7

Las cifras de un número entero en notación decimal son los coeficientes

de un polinomio cuyo valor en 10 es el número dado. Para pasar de base
2 a lase 10 basta evaluar el correspondiente polinomio en 2 . Por ejemplo
10112 = 23 22 0·2+1 = 13, Recı́procamente, las cifras de la expresión en base
2 de dicho número y sus cocientes. Expresa el número 53 en base 2.

Usando la función rem de Matlab, que imprime el residuo de una división de

enteros y fix que da la parte entera de una división, hasta obtener como residuo
cero:
1 n=53;
2 base =2;
3 i=1;
4 while n>0
5 c(i)=rem(n,base);
6 n=fix(n/base);
7 i=i+1;
8 end
9 i=i -1;
10 disp(c(i: -1:1))

Cuyo resultado es 1 1 0 1 0 1

23
 PREGUNTA 2.8

La función ex se aproxima en un entorno del origen mediante el polinomio

de taylor:

x2 x3 xn
pn (x) = 1 + x + + + ··· +
2 3! n!
Representa la función exponencial y sus sucesivas aproximaciones en un
intervalo de la forma [−a, a].

SOLUCIÓN

La función exponencial se representará mediante el siguiente código de matlab:

1 hold off
2 syms x
3 dom=linspace (-2,2,100);
4 fun = @(x) exp(x) ;
5 a=10;
6 Z=zeros(1,a+1);
7
8 n=0;
9 while n<=a
10 Zp=zeros(1,n+1);
11 P = diff(fun ,x,n)/factorial(n);
12 f=inline(P);
13 Z(1,n+1)=f(0);
14 Zq=polyval(fliplr(Z),dom);
15 plot(dom ,Zq)
16 hold on
17 n=n+1;
18 end
19 plot(dom ,fun(dom),'k','linewidth ' ,1.1)

24
 PREGUNTA 2.9

La función log(1 + x) se aproxima en un entorno del origen mediante el

polinomio de taylor:

x2 x3 xn
pn (x) = x − + − · · · + (−1)n+1
2 3 n
Representa la función log(1 + x) y sus sucesivas aproximaciones en el inter-
valo < −1, 1 > . Observa la diferencia con el problema anterior.

SOLUCIÓN

El procedimiento seguido es el mismo que el ejercicio anterior.

1 hold off
2 syms x
3 dom=linspace (-1,1,100);
4 fun = @(x) log(x+1) ;
5 a=10;
6 Z=zeros(1,a+1);
7
8 n=0;
9 while n<=a
10 Zp=zeros(1,n+1);
11 P = diff(fun ,x,n)/factorial(n);
12 f=inline(P);
13 Z(1,n+1)=f(0);
14 Zq=polyval(fliplr(Z),dom);
15 plot(dom ,Zq)
16 hold on
17 n=n+1;
18 end
19 plot(dom ,fun(dom),'k','linewidth ' ,1.1)

25
 PREGUNTA 2.10

Halla la descomposición en factores irreductibles reales del polinomio

x5 + x4 + x3 + x2 + x + 1

SOLUCIÓN

Para hallar la descomposición en factores irreductinles de este polinomio, prime-

ro se encontrarán las raı́ces, y a partir de estas se construirá el polinomio factori-
zado. Se usa la siguiente instrucción en Matlab:

>> p=[1 1 1 1 1 1];

>> roots(p)

ans =

0.5000 + 0.8660i
0.5000 - 0.8660i
-1.0000 + 0.0000i
-0.5000 + 0.8660i
-0.5000 - 0.8660i

Donde el vector respuesta ans contiene las raı́ces del polinomio enunciado.
Reconstruyendo el polinomio, se tiene:

p(x) =
(x − (0,5000 + 0,8660i)) (x − (0,5000 − 0,8660i)) (x − (−1,0000 + 0,0000i)) (x − (−0,5000 + 0,8660i)) (x

26
 REPRESENTACIÓN GRÁFICA
CAPÍTULO 3
PREGUNTA 3.1

Las ecuaciones parametricas de una elipse son:

x = a cos t
y = b sin t

donde a y b son los semiejes. Representa la elipse de ecuaciones cartesianas,

x2 y2
9 + 16 = 1, pasándola previamente a ecuaciones paramétricas y obtén en la
misma gráfica la hipérbola xy = 1 para obtener una aproximacion de los
puntos de corte de ambas cónicas.

SOLUCIÓN
Para graficar la curva conica se procede a parametrizar la ecuacion:
Donde:
cos t = x–α
a
y–β
sin t = b
x = 3 cos t
y = 4 sin t

Listing 3.1: Scripts de las representaciones de la elipse y hiperbola.

1 clear all
2 t=0: pi /100:2 * pi;
3 x=3 * cos(t);
4 y=4 * sin(t);
5 s1=plot(x,y)
6 axis equal
7 hold on
8 syms x y
9 f=char ((x * y) - 1);
10 s2=ezplot(f);
11 legend ( [ s1 ,s2 ] , 'elipse ' , 'hiperbola ' )
12 title('ELIPSE vs. HIPERBOLA ')
13 hold off

27
 ELIPSE vs. HIPÉRBOLA
6
elipse
hipérbola
4

y
-2

-4

-6
-6 -4 -2 0 2 4 6
x

Figura 3.1: Graficas de la elipse y la hipérbola juntas.

PREGUNTA 3.2

La funcion y(t) = 1−e−bt , donde t es el tiempo y b¿0,describe muchos proce-

sos de ingenierı́a, como la temperatura de un objeto que esta siendo calenta-
do, etc. Investiga el efecto del parámetro b en la funcion y(t), representando
y para varios valores de b en la misma gráfica.

SOLUCIÓN

-El parametro b es una constante de enfriamiento que depende del tipo de mate-
rial que se calienta.
-Este puede variar de 0 a 1, con valores muy pequeños.
-La temperatura hacemos variar de 0° a 1000°
-Para representar la gráfica, usamos los siguientes codigos.
Listing 3.2: Scripts de las representaciones de la elipse y hiperbola.
1 x=0:1000;
2 syms x y
3 for i=0:0.01:1
4 f=1-exp(-i * x)-y;
5 ezplot (f);
6 hold on;
7 axis equal
8 end

28
 1 - exp(-(7 x)/50) - y

y
-2

-4

-6

-6 -4 -2 0 2 4 6
x

Figura 3.2: Representacion de las asintotas con varios valores de b.

PREGUNTA 3.3

El potencial eléctrico creado en un punto P por una carga q es:

q
V =
4πε0 r

Donde r es la distancia de la carga al punto y 4πε1 0 r = 9x109 . El campo

→
−
eléctrico E es menos el gradiente del potencial. Trabajando en el plano
XY.
→
−
E (x, y) = ∇V (x, y)
Considera un dipolo formado por dos cargas iguales de signo opuesto y va-
lor q = 2x10−6 C, situadas en los puntos (1.5,-1.5) y (-1.5,1.5). El potencial
en un punto creado por dos cargas es la suma de los potenciales correspon-
dientes a cada una de ellas.
a) Representa en 3 dimensiones el potencial creado por el dipolo en una
región del plano que lo contenga.
b) Dibuja en el plano las curvas equipotenciales.

SOLUCIÓN

a) Representaremos el campo en el cuadrado [-5,5]x[-5,5] utilizando la orden

surf.
Listing 3.3: Scripts para graficar el campo en el cuadrado.
1 clear all
2
3 q=2e -6;
4 k=9e9;
5 a=1.5;
6 b= -1.5;
7 x= -5:1/3:5;
8 y=x;
9 [X,Y]= meshgrid(x,y);
10 rp=sqrt ((X-a).ˆ2+(Y-b).ˆ2);

29
11 rn=sqrt ((X+a).ˆ2+(Y+b).ˆ2);
12 V=q * k * (1./rp -1./ rn)
13 surf(X,Y,V)
14 pause

Obtenemos asi la figura:

105
1

0.5

-0.5

-1
5
5
0
0

-5 -5

Figura 3.3: Potencial electrico creado por un dipolo.

b) Debido a la gran variación de los valores del potencial, las curvas de nivel a
intervalos regulares se acumulan en torno a las cargas. Fijamos manualmente los
niveles a representar mediante contour.
1 clear all
2
3 N =10.ˆ(2.5:.5:4.5) ;
4 P=[- fliplr(N) 0 N];
5 contour(X,Y,V,P)

-1

-2

-3

-4

-5
-5 -4 -3 -2 -1 0 1 2 3 4 5

Figura 3.4: Campo y curvas equipotenciales.

30
 PREGUNTA 3.4

a) Representa la superficie engendrada por una circunferencia al girar en

torno a un eje situado en su mismo plano.
b) Representa la superficie generada por una recta que gira en torno a un
eje perpendicular a la misma, al tiempo que se desplaza en la dirección de
dicho eje.
c) Representa la superficie generada por una recta que gira en torno a otra
que se cruza con ella en el espacio.

SOLUCIÓN

a) Sea C una circunferencia situada en el plano XY con centro en (a,0,0) y radio

b, que gira en torno al eje Z. Las coordenadas de un punto(r,z) de la misma se
expresan en función del ángulo t como:

r = a + b cos t
z = b sin t

Con t ∈[0,2π]. Al girar en torno a Z, el punto (r,0,z) recorre las posiciones (x,y,z),
siendo.

x = r cos s
y = r sin s

Siendo s ∈[0,2π] el ángulo de giro. Traducimos a MATLAB esta parametrización

y representamos la superficie tórica.
1 clear all
2 a=2;
3 b=1;
4 s=0: pi /10:2 * pi;
5 t=0: pi /10:2 * pi;
6 [T,S] = meshgrid (t,s);
7 R=a+b * cos(T);
8 Z=b * sin(T);
9 X=R. * cos(S);
10 Y=R. * sin(S);
11 mesh(X,Y,Z)
12 axis equal
13 pause
14 surf(X,Y,Z)
15 axis equal

En la figura, se ha reducido el intervalo de variación de los parámetros con el fin

de mostrar la superficie tórica como deformación de un rectángulo.
b) Consideremos una recta situada sobre el eje OX que se desplaza en la direc-
ción el eje OZ al tiempo que gira en torno al mismo. Tomamos como parámetros

31
 1

-1
3
2
3
1 2
0 1
-1 0
-1
-2
-2
-3 -3

Figura 3.5: Generación de la superficie tórica.

la distancia del eje de giro y la altura t = z de un punto genérico P (x, y, z) de la su-

perficie. Suponemos que el ángulo girado es igual a t. Entonces, las coordenadas
de P son:
x = s cos t
y = s sin t
z=t

Para representar gráficamente la superficie, limitamos los valores s al intervalo

[-2,2] y los de t al intervalo [0,2π].
1 clear all
2 s= -2:0.2:2;
3 t=0: pi /20:2 * pi;
4 [T,S] = meshgrid(t,s);
5 X=S. * cos(T);
6 Y=-S. * sin(T);
7 Z=T;
8 mesh(X,Y,Z)
9 axis equal
c) Suponemos que la generatriz gira en torno al eje OZ, pasando inicialmente
por el punto (d,0,0) con dirección (0,cosα, senα). Un punto de esta recta puede
expresarse en función de un parametro s como:
xs,0 = d
ys,0 = s cos α
zs,0 = s sin α

Al girar en torno al eje OZ, la altura del punto P no varı́a, mientras que la abscisa
y la ordenada se transforman mediante la matriz del giro en el plano XY .
! ! !
xs,t cos t − sin t xs,0
=
ys,t sin t cos t ys,0

32
 6

0
2
2
0
0
-2 -2

Figura 3.6: Helicoide como superficie reglada.

Donde t es el angulo de giro. Por tanto, las coordenadas de un punto genérico de

la superficie son:

xs,t = d cos t − s cos α sin t

ys,t = d sin t + s cos α cos t
zs,t = s sin α

Eligiendo intervalos adecuados para los parámetros, resulta el siguiente código:

1 clear all
2 s= linspace (-1,1,10);
3 t= linspace (0,2 * pi ,25);
4 a=pi /6; d= 1
5 [T,S]= meshgrid(t,s);
6 X= d * cos(T)-S * cos(a). * sin(T);
7 Y= d * sin(T)+S * cos(a). * cos(T);
8 Z=S * sin(a);
9 surf(X,Y,Z)

La superficie obtenida es un hiperboloide de una hoja. Este ejercicio prueba que

se trata de una superficie reglada, es decir, generada mediante rectas.

33
 0.5

-0.5
2
1 2
0 1
0
-1
-1
-2 -2

Figura 3.7: Hiperboloide de una hoja como superficie reglada.

PREGUNTA 3.5

Una circunferencia rueda sin deslizarse sobre una recta.

a) Dibuja la trayectoria (x(t), y(t)) de un punto de la circunferencia mien-
tras avanza hasta dar una vuelta completa, deduciendo las ecuaciones pa-
ramétricas del movimiento de dicho punto
b) Estima la longitud de esta trayectoria midiendo la poligonal utilizada
para representarla gráficamente.
c) Suponiendo que gira a velocidad constante, obtén la expresión analı́tica
de las componentes de la velocidad (x0 (t), y 0 (t)) y dibuja el vector velocidad
sobre varios puntos de la trayectoria.

SOLUCIÓN

a) Consideramos una circunferencia como centro en el punto (0,1) y radio 1 que

gira avanzando sin deslizarse sobre el eje OX. Como el radio es 1, el desplaza-
miento t coincide con el ángulo girado. Salvo que el giro es en sentido negativo.
Las coordenadas del punto de la circunferencia que inicialmente estaba en el ori-
gen se pueden expresar en función de t como:

x = t − sin t
y = 1 − cos t

-Cuando t varı́a de 0 a 2π, la circunferencia da una vuelta completa, recorriendo

una distancia de 2π.
-Representamos gráficamente las ecuaciones paramétricas de la trayectoria del
punto P (x, y), denominada cicloide.

34
 1 clear all
2 t=linspace (0,2 * pi);
3 x=t-sin(t);
4 y=1-cos(t);
5 plot(x,y,[0 2 * pi],[0 0],'k')
6 axis equal
7 hold on
8 z=0: pi /30:2 * pi;
9 a=cos(z)+2;b=sin(z)+1;
10 plot(a,b)
11 hold off

2.5

1.5

0.5

-0.5

-1

0 1 2 3 4 5 6

Figura 3.8: Generación de la cicloide.

b) La longitud de la poligonal {(xk , yk )}n+1

k=1 utilizada para representar gráficamente
la cicloide es: Efectuamos los cálculos con MATLAB.

clear all
t=linspace(0,2*pi);
x=t-sin(t);
y=1-cos(t);
lx=diff(x);
ly=diff(y);
L=sum(sqrt(lx.ˆ2+ly.ˆ2))

L =

7.9997

-Se obtiene L = 79997, lo que constituye una buena aproximación del resultado
análitico, L = 8.

35
 c) Tomando t como el tiempo, las componentes de la velocidad del punto P son:

x0 = 1 − cos t
y 0 = sin t

Eligimos algunos puntos de la trayectoria y dibujamos en ellos los vectores velo-

cidad.
1 clear all
2 t=linspace (0,2 * pi);
3 x=t-sin(t);
4 y=1-cos(t);
5 u=linspace (0,2 * pi ,11);
6 sx=u-sin(u);
7 sy=1-cos(u);
8 dx=1-cos(u);
9 dy=sin(u);
10 plot(x,y,'r:',sx ,sy ,'.' ,[0 2 * pi],[0 0],'k')
11 axis equal
12 hold on
13 quiver (sx ,sy ,dx ,dy ,0.8)
14 hold off

El último argumento de quiver escala los vectores para obtener una representa-
ción más clara.

3.5

2.5

1.5

0.5

-0.5

-1

-1.5
0 1 2 3 4 5 6

Figura 3.9: Velocidad del punto que genera la cicloide.

36
 PREGUNTA 3.6

Muchos fenómenos aleatorios, como los errores de medidas experimenta-

les, se comportan según la llamada distribución de probabilidad normal. La
función densidad de esta distribución, en el caso de media 0 y desviación
tı́pica 1, N (0, 1), es:
1 2
f (x) = √ e−x /2
2π
Si X es una variable aleatoria,
R x la probabilidad de que X sea menor que un
valor x dado es: p(X ≤ x) = −∞ f (t)dt = F(x) Esta fórmula define la llamada
funcion de distribución F(x). La función randn de MATLAB genera valores
aleatorios que siguen la distribución N (0, 1).
a) Representa gráficamente la función de densidad de la distribución
N (0, 1) en el intervalo [−4, 4].
b) Genera 1000 resultados aleatorios con randn agrúpalos en intervalos y
representa las frecuencias relativas de cada subintervalo en un diagrama de
barras. Compara este diagrama con la gráfica de la función densidad.
c) Calcula las frecuencias relativas acumuladas de los subintervalos ante-
riores. Compara su gráfica con de la función de distribuci
 ón, F(x), que en
1
MATLAB puede obtenerse como: F(x) = 2 1 + erf( √x )
2
d) Simula el movimiento Browniano, construyendo dos vectores aleatorios
x e y, del mismo tamaño, e interpretando sus coordenadas como despla-
zamientos de una partı́cula en el plano XY . La posición de la partı́cula en
cada instante se obtiene sumando el desplazamiento (xk , yk ) a la posición
en el instante anterior (pk , qk ). Dibuja la traza de la partı́cula.

SOLUCIÓN

a) Las siguientes instrucciones de MATLAB producen la gráfica de la función

densidad de probabilidad de la distribucion normal.
1 t= -4:0.1:4;
2 y= 1/ sqrt (2 * pi) * exp(-t.ˆ2/2);
3 plot(t,y)

b) Generamos un vector aleatorio con randn. La orden previa sirva para obtener
siempre el mismo resultado .aleatorio”. Comprobamos que la media y la desvia-
ción tı́pica se aproximana los valores de la distribución teórica. Esta aproxima-
ción mejora tomando más valores en la muestra.
1 randn('state ' ,0)
2 x=randn (1 ,1000);
3 mean(x)
4 ans =
5 -0.0431
6 std(x)
7 ans =
8 0.9435

37
 -La orden hist de MATLAB agrupa los datos en intervalos y cuenta las fre-
cuencias. Tomamos intervalos centrados en los puntos ck = −4 + 0,4k, para k =
0, 1, 2, ..., 20. Representaremos las frecuencias e un diagrama de barras. Para com-
pararlo con la función densidad de probabilidad, el área de cada barra debe co-
rresponder a la probabilidad del intervalo correspondiente.
-La altura de la barra multiplicada por la longitud del subintervalo será la fre-
cuencia relativa del mismo.
-La gráfica pedida se obtiene pues mediante.
1 h=0.4;
2 c=-4:h:4;
3 frec= hist(x,c);
4 frel = frec /1000;
5 bar(c,frel/h)
6 hold on
7 plot(t,y)
8 hold off

0.45

0.4

0.35

0.3

0.25

0.2

0.15

0.1

0.05

0
-4 -3 -2 -1 0 1 2 3 4

Figura 3.10: Frecuencias de la muestra y función de densidad de probabilidad

N(0,1).

Observamos en la figura una buena aproximación global, aunque no perfecta al

trabajar con un número finito de valores.
c) Calculamos las frecuencias relativas acumuladas con cumsum y las representa-
mos junto con las probabilidades dadas por la función de distribución. Debido
a que la orden hist centra los intervalos en los puntos ck , cada intervalo va de
ck − h/2 a ck + h/2, siendo h = 0, 4 la longitud del subintervalo. La suma de las
frecuencias hasta el intervalo k-ésimo, corresponde al valor d ela función de dis-
tribucion F en ck +h/2. Por este motivo, desplazamos a la derecha en h/2 la gráfica
de las frecuancias acumuladas.
1 z=(1+ erf(t/sqrt (2)))/2;
2 plot (t,z)
3 hold on
4 frelac = cumsum(frel);
5 bar(c+h/2, frelac)
6 hold off

38
 1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
-4 -3 -2 -1 0 1 2 3 4

Figura 3.11: Frecuencias relativas acumuladas y función de distribución N(0,1).

En la figura muestra el ajuste casi perfecto de las frecuencias relativas acumula-

das y la función de distribución.
d) Simulamos el movimiento Browniano según las instrucciones del enunciado,
obteniendo la figura siguiente. Utilizando comet en lugar de plot, puedes ver
una animacion del movimiento.
1 x= randn (1 ,1000); y= randn (1 ,1000);
2 p= cumsum(x); q= cumsum(y);
3 comet(p,q)
4 axis equal

25

20

15

10

-5

-10

-15

-20

-20 -10 0 10 20 30 40

Figura 3.12: Movimiento Browniano.

39
 LA ECUACIÓN F(x) = 0
CAPÍTULO 4
PREGUNTA 4.1

Determinar las raices del polinomio p(x) = 2x3 − 3x − 1 mediante el metodo

de la biseccion mostrando el proceso paso a paso.

SOLUCIÓN

-Representamos gráficamente el polinomio en el intervalo [-2,2]:

1 clear all
2 x=linspace (-2,2);
3 p=[2 0 -3 -1];
4 y=polyval(p,x);
5 plot (x,y)
6 grid on

-Gráficamente parece que x = −1 es una raı́z del polinomio p(x), en efecto p(−1) =
0. Observando la figura 1 partiremos del intervalo [a1 , b1 ] = [0, 2] para hallar la
raı́z positiva. Procediendo de forma análoga, se obtendrá la otra raı́z negativa.
Verificamos que el polinomio tiene distinto signo en los extremos del intervalo.

10

-5

-10

-15
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

Figura 4.1: p(x) = 2x3 − 3x − 1

p(0) = −1 < 0 y p(2) = 9 > 0

40
Calculamos el punto medio del intervalo y evaluamos en éste polinomio.

c1 = (a1 + b1 )/2 = 1 y p(1) = −2 < 0

A continuación,comprobamos que p(b1 )∗p(c1 ) < 0 por lo que el intervalo [c1 , b1 ] =

[1, 2] contiene a la raı́z, nos centramos en él haciendo a2 = c1 y b2 = b1 , realizamos
una nueva iteración en [a2 , b2 ] = [1, 2]:

c2 = (a2 + b2 )/2 = 1,5, p(c2 ) = 1,25

Ahora p(a2 )∗p(c2 ) < 0 entonces b3 = c2 y a3 = a2 nos quedamos con el subintervalo

izquierdo ya que la raı́z está en [a3 , b3 ] = [1, 1,5] cuyo punto medio.

c3 = (a3 + b3 )/2 = 1,25, p(1,25) = −0,843750

por lo que haremos a4 = c3 y b4 = b3 , nos centraremos en [a4 , b4 ] = [1,25, 1,5] y con-

tinuaremos el proceso hasta que la longitud del intervalo [ak , bk ] sea tan pequeña
como queramos, en este caso tras 25 iteraciones se alcanzo la solución 1.366025,
para tolerancia de 10−7 .
En la figura 2 mostramos los primeros iterados. La ventaja del algoritmo de bi-

16

14

12

10

0
X: 1.366
-2 Y: 0.001117

-4
-0.5 0 0.5 1 1.5 2 2.5

Figura 4.2: Metodo de biseccion.

sección es que siempre converge, si partimos de un intervalo que contiene a la

raı́z. Sin embargo, la convergencia es lenta y puede darse el caso de que el punto
medio calculado en un paso esté más alejado de la raı́z que el iterado del paso
anterior.
PREGUNTA 4.2

2
Halla la menor raı́z positiva de la ecuación e−x −cos x = 0 con una tolerancia
de 10−5 , utilizando el método de regula falsi y mostrando los pasos en las
primeras iteraciones.

SOLUCIÓN

En primer lugar, representamos gráficamente la función para elegir el intervalo

de partida.

41
1 clear all
2 x= -3:0.1:3;
3 y=exp (-x.ˆ2) -cos (x);
4 plot(x,y), grid on

Observando la gráfica de la función (Figura 3) partimos del intervalo [a1 , b1 ] =

[1, 2] para hallar la raı́z.
Evaluamos la función en los extremos:
2
f (a1 ) = e−a1 − cos a1 = −0,172423
2
f (b1 ) = e−b1 − cos b1 = 0,434462

42
 1

0.8

0.6

0.4

0.2

-0.2
-3 -2 -1 0 1 2 3

2
Figura 4.3: e−x − cos x = 0 .

Calculamos el punto c1 donde la recta que une los puntos (a1 , f (a1 )) y (b1 , f (b1 ))
corta al eje x:

c1 = (a1 ∗ f (b1 ) − b1 ∗ f (a1 ))/(f (b1 ) − f (a1 )) = 1,284111

f (c1 ) = −0,090521
A continuación,comprobamos en cuál de los subintervalos [a1 , c1 ] y [c1 , b1 ] se en-
cuentra la raı́z y nos centramos en él. Como f (b1 ) ∗ f (c1 ) < 0 entonces a2 = c1
y b2 = b1 , es decir, la raı́z se encuentra en [a2 , b2 ] = [1,28411, 2], realizamos una
nueva iteración.

c2 = (a2 ∗ f (b2 ) − b2 ∗ f (a2 ))/(f (b2 ) − f (a2 )) = 1,407549

f (c2 ) = −0,024619
De nuevo f (b2 ) ∗ f (c2 ) < 0 y por tanto, hacemos a3 = c2 y b3 = b2 ,nos quedamos
con el subintervalo derecho [a3 , b3 ] = [1,284111, 1,407549] y la siguiente iteración
será:
c3 = (a3 ∗ f (b3 ) − b3 ∗ f (a3 ))/(f (b3 ) − f (a3 )) = 1,439320
f (c3 ) = −0,005119
Por lo que haremos a4 = c3 y b4 = b3 , nos centraremos en [a4 , b4 ] = [1,407549, 1,439320]
y continuaremos el proceso hasta que la longitud del intervalo [ak , bk ] sea tan
pequeña como queramos, en este caso tras 8 iteraciones,se alcanza la solución
1.447414 para tolerancia e 10−5 .
Si exigimos una tolerancia de 10−7 , la solución 1.447414 se alcanza en 11 itera-
ciones.
PREGUNTA 4.3

Utiliza el método de Newton para aproximar la raı́z quinta de 2 con tres

cifras decimales exactas, ilustrando el método paso a paso.

SOLUCIÓN
√
Notemos que 5 2 es la solución de la ecuación f (x) = x5 − 2 = 0, cuya derivada
f 0 (x) = 5x4 utilizamos en el método de Newton.

43
 Tomando una estimación inicial x1 = 1 y aplicando el método de Newton calcu-
lamos las siguientes iteraciones con

x15 − 2 −1 6
x2 = x1 − = 1− = = 1,2
5x14 5 5

x25 − 2 1,22 − 2
x3 = x2 − = 1,2 − = 1,152901
5x24 5 ∗ 1,24
Siguiendo este proceso con dos iteraciones más tenemos x5 = 1,148698, valor que
podemos comparar con el obtenido directamente al hacer 21/5 y observar que se
tienen más de tres cifras exactas.
b) Utilizamos el fichero newton.m para sistematizar los cálculos. En primer lugar,
creamos la función raiz5.m con dos parámetros de salida, la función f (x) de la
que se va a hallar una raı́z y su derivada f 0 (x)
1 function [sol ,x,incr ,k]= newton(fun ,x1 ,tol ,maxiter)
2 k = 1;
3 x = x1;
4 incr = tol +1;
5 while incr > tol & k <= maxiter
6 [y,dy] = feval (fun ,x(k));
7 delta = -y/dy;
8 x(k+1) = x(k) + delta;
9 incr = abs(delta);
10 k = k+1;
11 end
12 if incr > tol
13 sol = [ ];
14 disp('Diverge o insuficientes iteraciones ')
15 else
16 sol = x(end);
17 end

Creamos la función raiz5.m.

1 function [y,dy] = raiz5(x)
2 y = x.ˆ5 -2;
3 dy=5 * x.ˆ4;

[sol,x,incr,k] = newton(’raiz5’,1,1e-4,31)
sol =
1.1487
x =
1.0000 1.2000 1.1529 1.1487 1.1487
incr = error

3.0530e-05
k =número de iteraciones
5

44
 PREGUNTA 4.4

La raı́z de la ecuación tan(x) − c = 0 es relativamente difı́cil de encontrar

cuando c es una cantidad grande. Trata de averiguar la razón de éste com-
portamiento a través de los siguientes pasos:
a) Dibuja la funcion f (x) = tan(x)−c para c = 5 y c = 10 y estima un intervalo
en que se encuentre la raı́z buscada.

SOLUCIÓN

1 % Grafica de funciones en Matlab

2 clc ,clear
3 x = linspace (0,2 * pi);
4 y = tan(x) -5;
5 z = tan(x) -10;
6 plot(x,y)
7 grid on
8 hold on
9 plot(x,z)

El intervalo estimado es de 4rad a 4,7rad.

60

40

20

-20

-40

-60

-80
0 1 2 3 4 5 6 7

Figura 4.4: Gráficos de las funciones ”tan(x)-c”.

b) Aproxima, para cada uno de los valores de c especificados en el apartado a), la

solución a la ecuación no lineal mediante los métodos de bisección y regula falsi
empleando el intervalo estimado en el apartado a), con una tolerancia de 10−4 .
Por el método de bisección se usa el siguiente código:
1 clear all
2 format short;
3 a=input('Introduzca el valor de a: ');
4 b=input('Introduzca el valor de b: ');
5 cont=input('Introduzca el numero de iteraciones cont: ');
6 fun=input('Introduzcal a funcion f(x)=','s');
7 f=inline(fun);
8 for k=1: cont
9 c=(a+b)/2;

45
10 e=abs((b-a)/2);
11 A(k,:)=[k a b c f(c) e];
12 if f(a) * f(c) <0
13 b=c;
14 else
15 a=c;
16 end
17 end
18 fprintf('\n \tk \ta \tb \tc \tf(c) \terror \n')
19 disp(A)
20 fprintf('Solucion :\n c= %8.5f\n',c)
21 fprintf('f(c)= %8.5f\n',f(c))
22 fprintf('error = %8.5f\n',e)

Solución:Para tan(x) − 5
x= 4.51508
f(x)= 0.00227
error= 0.00019
Solución:Para tan(x) − 10
x= 4.61287
f(x)= 0.01534
error= 0.00019

1 function calcular_Callback (hObject , eventdata , handles)

2 F1=str2double(get(handles.contf1 ,'String '));
3 F2=str2double(get(handles.contf2 ,'String '));
4 Fx=get([Link] ,'string ');
5 x1=str2double(get(handles.X1 ,'string '));
6 x2=str2double(get(handles.X2 ,'string '));
7 tol=str2double(get([Link] ,'String '));
8 prod=F1 * F2;
9 if prod <=0
10 e=10;
11 i=1;
12
13 while e>tol;
14 if i>50
15 set([Link] ,'string ','paso la iteracion
50');
16 break
17 end
18 A(i,1)=x1;
19 A(i,2)=x2;
20 m=x1 -(F2 * (x2 -x1))/(F1 -F2);
21 A(i ,5)=m;
22 x=m;
23 Fm=eval(Fx);
24 A(i ,6)=Fm;
25 x=x1;
26 F1=eval(Fx);

46
27 A(i ,3)=F1;
28 x=x2;
29 F2=eval(Fx);
30 A(i ,4)=F2;
31 e=abs(m-x1);
32 if F1 * Fm <0
33 x2=m;
34 else
35 x1=m;
36 end
37 i=i+1;
38 end
39 format long
40 xff=m;
41 set([Link] ,'data ',A);
42 set([Link] ,'string ',xff);
43 set([Link] ,'string ',e);
44 else
45 set([Link] ,'string ','NO CONVERGE ');
46 end

Solución:Para tan(x) − 5
x= 4.51499
f(x)= -0.000089
error= 0.0000029
Solución:Para tan(x) − 10
x= 4.61269
f(x)= -0.003062
error= 0.00000512
c) Analiza razonadamente la velocidad de convergencia de ambos métodos para
el intervalo estimado en el apartado a).
-Por el método de bisección la velocidad de convergencia es más rapida con res-
pecto al método de regula falsi , por el contrario el error del método de bisección
es mayor con respecto al método de regula falsi.

PREGUNTA 4.5

Se construye una caja sin tapadera a partir de una hoja metálica reactan-
gular que mide 10 por 16 centı́metros. Queremos averiguar cuál debe ser
el lado de los cuadrados que hay que recortar en cada esquina para que el
volumen de la caja sea 100 centı́metros cúbicos.

SOLUCIÓN

a) Plantea la ecuación que nos proporcionará el volumen de la caja en función del

lado L de los cuadrados a recortar.

47
 V ol = (16 − 2L)(10 − 2L) ∗ L
V ol = 4L3 − 52L2 + 160L
V ol = 100cm3
4L3 − 52L2 + 160L − 100 = 0

b) Representa gráficamente la ecuación a resolver y obtén un intervalo inicial en

el que esté el valor de L que estamos buscando.

1 % Grafica de funciones en Matlab

2 clc ,clear
3 x = linspace (0 ,10);
4 y = (4 * (x.ˆ3)) -(52 * (x.ˆ2))+(160 * x) -100;
5 plot(x,y)
6 grid on

300

250

200

150

100

50

-50

-100

-150

-200
0 1 2 3 4 5 6 7 8 9 10

Figura 4.5: Gráfico de la función del volumen.

El intervalo que tomaremos será de 3 a 4cm.

c) Resuelve dicha ecuación mediante el método de la bisección, empleando una
tolerancia de 10−9 y tomando como intervalo inicial el obtenido en el apartado b).

1 % m t o d o de la b i s e c c i n
2 function calcular_Callback (hObject , eventdata , handles)
3 F1=str2double(get(handles.contf1 ,'String '));
4 F2=str2double(get(handles.contf2 ,'String '));
5 Fx=get([Link] ,'string ');
6 x1=str2double(get(handles.X1 ,'string '));
7 x2=str2double(get(handles.X2 ,'string '));
8 tol=str2double(get([Link] ,'String '));
9 prod=F1 * F2;
10 if prod <=0
11 e=10;
12 i=1;

48
13
14 while e>tol;
15 if i>50
16 set([Link] ,'string ','paso la iteracion
50');
17 break
18 end
19 A(i,1)=x1;
20 A(i,2)=x2;
21 m=(x1+x2)/2;
22 A(i ,5)=m;
23 x=m;
24 Fm=eval(Fx);
25 A(i ,6)=Fm;
26 x=x1;
27 F1=eval(Fx);
28 A(i ,3)=F1;
29 x=x2;
30 F2=eval(Fx);
31 A(i ,4)=F2;
32 e=abs(m-x1);
33 if F1 * Fm <0
34 x2=m;
35 elseif F1 * Fm >0
36 x1=m;
37 elseif F1 * Fm==0
38 xff=m;
39 e=0;
40 break
41 end
42 i=i+1;
43 end
44 format long
45 xff=m;
46 set([Link] ,'data ',A);
47 set([Link] ,'string ',xff);
48 set([Link] ,'string ',e);
49 if e==0;
50 set([Link] ,'string ','EXACTO PUTIKIZ !!!
');
51 end
52 else
53 set([Link] ,'string ','NO CONVERGE ');
54 end
En 30 iteraciones se obtiene la raiz con una tolerancia de 10−9 .
Solución:
L=3.40175cm
f(L)=-0.0000743
Error=9,31 ∗ 10−10
d) ¿Mejora el método de regula falsi el resultado obtenido por bisección?.

49
 No. Pero nesecita 69 iteraciones para llegar a ese mismo resultado.
Solución:
L=3.40175cm
f(L)=-0.0000743
Error=2,68 ∗ 10−10
e) Obtén el lado de dichos cuadrados empleando el método de punto fijo, con
estimacion inicial el punto medio del intervalo utilizado en apartados anteriores.

1 % Metodo de punto fijo

2 clc ,clear
3 function PuntoFijo(x0 , es , imax , gx)
4 xr = x0;
5 iter = 0;
6 g = inline(gx);
7 do = 0;
8 while (do == 0)
9 xrold = xr;
10 xr = g(xrold);
11 iter = iter + 1;
12 if (xr ˜= 0)
13 ea = abs((xr - xrold)/xr) * 100;
14 end
15 if ((ea < es) || (iter >= imax))
16 break;
17 end
18 end
19 disp('Resultado ')
20 xr

Solución:
L=3.5344cm
f(L)=-7.472471
Error=−3,899 ∗ 10−2

PREGUNTA 4.6

Consideremos el siguiente sistema no lineal:

x2 + x − y 2 = 1
y − sen(x2 ) = 0

SOLUCIÓN

a) Determina gráficamente todas las soluciones del sistema.

b) Transforma el sistema de ecuaciones en una única ecuacion no lineal y aproxi-
ma con precisión de 10−6 sus soluciones mediante el método de punto fijo. Indica
las aproximaciones iniciales utilizadas y las iteraciones necesarias para cada una

50
 3.5

2.5

1.5

0.5
X: 0.7299
0 X: -1.682 Y: 0.5079
Y: 0.3078
-0.5

-1
-4 -3 -2 -1 0 1 2 3 4

Figura 4.6: Gráficos del sistema de ecuaciones no lineales .

de la raı́ces.

x2 + x − sin2 (x2 ) = 1
Las aproximaciones iniciales son:
a=0.7
b=-1.6
Para la primera raı́z se nesecitó solo 5 iteraciones.
1 >> PuntoFijo (0.7, 0.000001 , 5,'-x.ˆ2+1+( sin(x.ˆ2)).ˆ2')
2 Resultado
3
4 xr =
5
6 0.7260

Para la otra raiz se necesitó 20 iteraciones.

1 >> PuntoFijo (-1.6, 0.000001 , 20,'-x.ˆ2+1+( sin(x.ˆ2)).ˆ2')
2 Resultado
3
4 xr =
5
6 -1.6701

c) Repite el apartado b) empleando el método de Newton y comprueba la conver-

gencia cuadrática del método.
1 function evaluar_Callback (hObject , eventdata , handles)
2 a=get([Link] ,'string ');
3 syms x ;
4 dx=diff(a,x);
5 set([Link] ,'string ',char(dx));
6 x0=str2double(get([Link] ,'String '));
7 f=get([Link] ,'string ');
8 Df=get([Link] ,'string ');
9 T=str2double(get([Link] ,'String '));
10 h=1;
11 n=0;

51
12 while h>T;
13
14 if n>20 ;
15 break
16 else
17 b=n+1;
18 A(b,1)=x0;
19 x=A(b,1);
20 A(b,2)=eval(f);
21 A(b,3)=eval(Df);
22
23
24 xa=A(b,1) -(A(b,2)/A(b,3));
25 A(b,4)=xa -A(b,1);
26 x0=xa;
27 h=abs(A(b,4));
28 n=n+1;
29 end
30 format long
31 xff=xa;
32 set([Link] ,'data ',A);
33 set([Link] ,'string ',xff);
34 set([Link] ,'string ',h);
35
36
37 end

Solución:

f (x) = x2 + x − sin2 (x2 )

f 0 (x) = 2x + 1 − 2x sin(2x2 )
x0 = 0,7
tol = 0,000001
xr = 0,725951
error = 1,92 ∗ 10−7

Halla la menor raı́z positiva de la ecuación.

2
e−x − cos(x) = 0

con una tolerancia de 10−6 , utilizando los algoritmos de bisección .

1 % m t o d o de la b i s e c c i n
2 function calcular_Callback (hObject , eventdata , handles)
3 F1=str2double(get(handles.contf1 ,'String '));
4 F2=str2double(get(handles.contf2 ,'String '));
5 Fx=get([Link] ,'string ');
6 x1=str2double(get(handles.X1 ,'string '));
7 x2=str2double(get(handles.X2 ,'string '));
8 tol=str2double(get([Link] ,'String '));

52
 9 prod=F1 * F2;
10 if prod <=0
11 e=10;
12 i=1;
13
14 while e>tol;
15 if i>50
16 set([Link] ,'string ','paso la iteracion
50');
17 break
18 end
19 A(i,1)=x1;
20 A(i,2)=x2;
21 m=(x1+x2)/2;
22 A(i ,5)=m;
23 x=m;
24 Fm=eval(Fx);
25 A(i ,6)=Fm;
26 x=x1;
27 F1=eval(Fx);
28 A(i ,3)=F1;
29 x=x2;
30 F2=eval(Fx);
31 A(i ,4)=F2;
32 e=abs(m-x1);
33 if F1 * Fm <0
34 x2=m;
35 elseif F1 * Fm >0
36 x1=m;
37 elseif F1 * Fm==0
38 xff=m;
39 e=0;
40 break
41 end
42 i=i+1;
43 end
44 format long
45 xff=m;
46 set([Link] ,'data ',A);
47 set([Link] ,'string ',xff);
48 set([Link] ,'string ',e);
49 if e==0;
50 set([Link] ,'string ','EXACTO PUTIKIZ !!!
');
51 end
52 else
53 set([Link] ,'string ','NO CONVERGE ');
54 end

Solución:

53
 x1 = 1 (4.1)
x2 = 2 (4.2)
tol = 0,000001 (4.3)
xr = 1,44741 (4.4)
error = 9,54 ∗ 10−7 (4.5)

La menor raiz positiva es 1.44741 sin considerar al cero como raiz positiva.
Se necesitaron 20 iteraciones para llegar al resultado.

PREGUNTA 4.7

Consideremos el siguiente sistema no lineal:

ex ey + x cos(y) = 0
x+y −1 = 0

SOLUCIÓN

a) Analiza gráficamente las raices del sistema para x ∈ [−6, 6].

1 clear;
2 clc;
3 % Grafica de funciones en Matlab
4 x = linspace (-2 * pi ,2 * pi);
5 y = 1-x;
6 plot(x,y)
7 grid on
8 hold on
9 syms x
10 syms y
11 ezplot(exp(x) * exp(y)+x * cos(y))
12 grid on

x+y-1

X: -4.379
4 Y: 5.379

2
y

X: 3.491
Y: -2.491
-2

-4
X: 5.141
Y: -4.141
-6
-6 -4 -2 0 2 4 6
x

Figura 4.7: Gráficos del sistema de ecuaciones no lineales .

54
 b) Transforma el sistema de ecuaciones en una única ecuación no lineal y aproxi-
ma las raı́ces encontradas en el apartado anterior mediante el método de Newton
con una tolerancia de 10−5 . ¿Cuántas iteraciones son necesarias en cada caso?. Es-
tima la velocidad de convergencia en una de la raı́ces.
-Transformando en una sola ecuación no lineal.

y = 1−x
Reemplazando
f (x) = ex e1−x + x ∗ cos(1 − x) = 0
La − derivada − de − la − f uncion :
f 0 (x) = x ∗ sin(1 − x) + cos(1 − x) = 0

-Se usará el mismo codigo de matlab del fichero newton.m para sistematizar los
cálculos.
-Luego creamos la función raiz1.m para cada raiz con dos parámetros de salida,
la función f (x) de la que se va a hallar una raı́z y su derivada f 0 (x).
1 function [y,dy] = raiz1(x)
2 y = exp(x) * exp(1-x)+x * cos(1-x);
3 dy= x * sin(1-x)+cos(1-x);

-La velocidad de convergencia de las raices son de 0.27 segundos cada una. -Las
iteraciones necesarias para cada caso es de 4,5 y 6 iteraciones según cada caso.
1 [sol ,x,incr ,k] = newton('raiz1 ' ,3,1e-5 ,40)
2 sol =
3 3.4706
4 x =
5 3.0000 3.4675 3.4706 3.4706
6 incr =error
7 2.4369e-06
8 k =iteraciones
9 4
10 [sol ,x,incr ,k] = newton('raiz2 ' ,5,1e-5 ,40)
11 sol =
12 5.1572
13 x =
14 5.0000 5.1757 5.1574 5.1572 5.1572
15 incr =error
16 2.0577e-08
17 k =iteraciones
18 5
19 [sol ,x,incr ,k] = newton('raiz3 ' ,9,1e-5 ,40)
20 sol =
21 9.1554
22 x =
23 9.0000 9.1557 9.1554 9.1554
24 incr =error
25 2.6022e-09

55
26 k =iteraciones
27 4
28 [sol ,x,incr ,k] = newton('raiz4 ' ,11,1e-5 ,40)
29 sol =
30 11.7624
31 x =
32 11.0000 12.2656 11.7700 11.7624 11.7624
11.7624
33 incr =error
34 3.0401e-11
35 k =iteraciones
36 6

56
APLICACIONES DE LOS SISTEMAS LINEALES
CAPÍTULO 5
PREGUNTA 5.1

El diagrama adjunto de la figura representa la discretización del problema

del calor en una placa. se trata de determinar la temperatura en los nodos
interiores de la malla tj; j01, 2, 3, ...,6, conocidas las temperaturas en el bor-
de y suponiendo que hay equilibrio térmico , es decir que las temperaturas
no varı́an.

En el modelo discreto se supone que , en el equilibrio, la temperatura en

cada nodo es la media de las temperaturas en los nodos vecinos.
A) Obtén el sistema lineal correspondiente a los datos de la figura, formu-
lando las condiciones que han de verificar las temperaturas de los nodos
interiores.
B)Itera por el método de Jacobi hasta que la verificación máxima de la tem-
peratura en en nodo sea inferior a 0.01. ¿cuál es la máxima desviación con
respecto a la solución obtenida en el primer apartado?

SOLUCIÓN

57
 PREGUNTA 5.2

Una armadura es una estructura que se usa para soportar en puentes o

edificios. La armadura esta integrada en su totalidad por miembros rectos
dispuestos de modo que forman uno o más triángulos, dando estabilidad a
la estructura. En la figura ?? se muestra una armadura Tı́pica.

En el ejemplo , hay 6 articulaciones y 9 miembros. La fuerza que actúa en

cada articulación, se descompone en sus componentes x e y. Las fuerzas
que actúan sobre la armadura están indicadas a lo largo de cada uno de
los [Link] considera que las fuerzas actúan hacia el centro del miem-
bro,lejos de la articulación. Hay 9 fuerzas Fi , i = 1, 2, ..., 9. Al hacer cero la
suma de todas las fuerzas que actúan horizontal o verticalmente en cada
una de las articulaciones, se obtiene un sistema de ecuaciones lineales 9 ∗ 9;
por ejemplo, en la primera articulación podemos escribir:

sin(45)◦ F1 − F4 − cos(30)◦ F5 = 0
sin(45)◦ F1 + F3 − sin(30)◦ F5 = −1000
A) Escribe El sistema de ecuaciones para determinar las fuerzas indicadas
en la armadura:

SOLUCIÓN

para el cálculo de las fuerzas tomamos encuentra las direcciones de las fuerzas
que nos dan en la figura ??, se muestra el sistema de fuerzas y por condición del
problema todas las fuerzas a compresión esto para el diagrama y el cálculo.
Los diagramas de fuerza de cuerpo libre para cada nodo se muestra en la figura
5.1 las sumas de las fuerzas en ambas direcciones, vertical y horizontal, deben
de ser cero en cada nodo se considera que una fuerza positiva va de izquierda a
derecha y de abajo hacia arriba. las ecuaciones se calculara por nodos:
P
1. Nodo a Fh = 0
◦
P1 + F2 + F1 cos 45 = 0
V
FV = 0
U1 + F1 sin 45◦ = 0
P
2. Nodo b Fh = 0
−F1 sin 45◦ + F4 + F5 cos 30◦ = 0
P
FV = 0
−F3 − F1 cos 45◦ − F5 sin 30◦ = 1000

58
 Figura 5.1: Diagrama de cuerpo libre de la armadura

P
3. Nodo c Fh = 0
−F2 + F6 = 0
P
FV = 0
F3 = 0
P
4. Nodo d Fh = 0
−F4 + F9 sin 45◦ = −500
P
FV = 0
−F7 − F9 cos 45◦ = 0
P
5. Nodo e Fh = 0
−F6 − F5 cos 30◦ + F8 = 0
P
FV = 0
F7 + F5 sin 30◦ = 0
P
6. Nodo f Fh = 0
◦
P 8 − F9 cos 45 = 0
−F
FV = 0
U2 + F9 sin 45◦ = 0
ahora el sistema de ecuaciones es:
V1 + F2 + 0,71F1 = 0
U1 + 0,71F1 = 0
−0,71F1 + F4 + 0,87F5 = 0
−F3 − 0,71F1 − 0,5F5 = 1000
−F2 + F6 = 0
F3 = 0
−F4 + 0,71F9 = −500
−F7 − 0,71F9 = 0
−F6 − 0,87F5 + F8 = 0
F7 + 0,5F5 = 0
−F8 − 0,71F9 = 0
U2 + 0,71F9 = 0

de este sistema de ecuaciones solo despejaremos con respecto a los 9 fuerzas que
se pide por pregunta se muestra en la figura 5.2
B)Resuelva el problema por factorización LU Antes de hallar mostraremos la
ecuación dada por LU.

AX = b Ly = b, Ux = y factorizacion LU (5.1)

59
Figura 5.2: se muestra el sistema de ecuaciones la matriz A, la vector de fuerzas
f) y vector de de coeficientes b

Procederemos a calcular la matriz triangular inferior (L) y la superior (U) en

Matlab [L, U ] = lu(A) se muestra las matrices declaradas dispersas por el orden
que tiene 9*9 y por los ceros que presenta. Calculamos:

Figura 5.3: factorización LU , L triángulo inferior y U triángulo superior declara-

dos dispersos co el comando Matlab sparse

b
y=
L
luego de hallar el vectorb procedemos calcular las fuerzas que sera despejando
x = (F) en la segunda ecuación de la factorización LU.
x = Uy
A continuación se muestra las fuerzas de cada miembro con un método iterativo
directo de factorización LU.

60
Figura 5.4: se muestra las fuerzas calculadas con el método directo de FACTORI-
ZACIÓN LU

61
VALORES PROPIOS Y VALORES SINGULARES
CAPÍTULO 6
PREGUNTA 6.1

Vamos a usar el método de potencias para usar determinar el valor y vector

propio dominante.

 3 0 −5
 
A = 0,2 −1 0 
 
1 1 −2
 

A) aplica el método de las potencias para calcular el valor propio dominan-

te de la matriz A, ası́ como un vector propio asociado.
método de las potencias Es una técnica iterativa para obtener de forma
aproximada valores y vectores propios de una matriz.
definición Si λ es un autovalor de A que, en valor absoluto, es mayor que
cualquier otro autovalor, se llama autovalor dominante y sus autovectores
correspondientes se llaman autovectores dominantes para este método se
toma un vector inicial de la dimensión de la matriz cuadrada.
n o
x0 = n m p

y se genera la siguiente ecuación:

Yk = Axk
1
xk = ( Yk )
ck+1

SOLUCIÓN

este proceso lo aremos en Matlab:

function [ lambda , V]= p o t e n c i a (A, X , t o l , maxiter )

% datos
% A e s una m a t r i z de o r d e n n * n
% X es el vector i n i c i a l de o r d e n n*1
% Tol es la t o l e r a n c i a
%max1 e s e l numero maximo d e i t e r a c i o n e s

62
% resultados
% lambda e s e l a u t o v a l o r dominante
%V e s e l a u t o v e c t o r dominante

lambda =0;
c n t =0;
e r r o r= t o l +1;
i =0;
while e r r o r > t o l & i < maxiter
i = i +1;
Y=A*X ;
% n o r m a l i z a c i o n de y
[m j ]=max( abs ( Y ) ) ;
c1=m;
dc=abs ( lambda−c1 ) ;
Y=(1/ c1 ) * Y ;
%A c t u a l i z a c i o n de X y c r i t e r i o de c o n v e r g e n c i a
dv=norm ( X−Y ) ;
e r r o r=max( dc , dv ) ;
X=Y ;
lambda=c1 ;

end
V=X ;
valor inicial: n o
x0 = 1 1 1
matriz para el cálculo:
 3 0 −5
 
A = 0,2 −1 0 
 
1 1 −2
 

la respuesta se muestra en la figura 6.1 siguiente

Figura 6.1: se da el valor propio dominante de la matriz A, ası́ como un vector

propio asociado ´por el método de potencia

63
 PREGUNTA 6.2

la figura corresponde a un sistema mecánico de masa m1 y m2 . Ligadas a

tres resortes con rigidez constante KK , i = 1, 2, 3. Denotamos por X1 y X2
los desplazamientos de las masas con respecto a su posición de equilibrio.
La Ley de Hooke, que expresa que la fuerza de reacción elástica es propor-
cional al desplazamiento y de sentido contrario a él, F = −Kx, nos da las
ecuaciones del equilibrio para el sistema.

SOLUCIÓN
A) Escribe la ecuación de equilibrio Primero se mostrara un diagrama de cuerpo
libre para las dos masas y con la Ley de Hooke F = −Kx calcularemos las ecua-
ciones donde los desplazamientos tienen sentido contrario a la fuerza (fuerza
restauradora) estas se muestra en la figura 6.4.

Figura 6.2: Diagrama del cuerpo libre

F1 = −K1 X1 − K2 (X1 − X2 )
F2 = K2 (X1 − X2 ) − K1 X1
Donde X1 y (X1 − X2 ) son las deformaciones que afectan a dichas masas.
B) da valores a las constantes elásticas y halla los desplazamientos X1 y X2 que
sufren las masas a sendas fuerzas, f1 y f2 . Estas fuerzas, en el equilibrio, se com-
pensan con las fuerzas de reacción elástica :f1 = F1
Kg Kg Kg
en nuestro caso tomaremos la constante K1 = 10 s2 ,K2 = 20 s2 y K3 = 15 s2 , las
fuerzas F1 = 19,6N y F2 = 29,4N
19,6 = −10X1 − 20(X1 − X2 )

64
 29,4 = 20(X1 − X2 ) − 15X1
Dando la forma de un sistema de ecuaciones lineales:

19,6 = −30X1 − 20X2

29,4 = 5X1 − 20X2

llevamos a la forma matricial:
matriz para el cálculo: " #
19,6
F=
29,4
" #
−30 −20
B=
5 −20
para hallar los desplazamientos tomaremos el método LU cuya ecuación es:

AX = b Ly = b, Ux = y factorizacion LU (6.1)

en nuestro caso al vector columna F = b y la matriz de coeficientes B = A obten-

dremos la matriz triangular inferior y superior A.

Figura 6.3: se observa el cálculo de la matriz triangular inferior y superior en

Matlab

procedemos a calcular y = b/L y luego los desplazamientos X = y/U

Figura 6.4: la figura muestra las fuerzas y los desplazamientos

65
INTERPOLACIÓN POLINÓMICA
CAPÍTULO 7
PREGUNTA 7.1

La tabla adjunta proporciona el numero de viajeros que suben al metro en

una estación a determinada hor a de lamañana.

hora 8 9 10 11 12 13 14
viajeros 41 35 21 9 11 17 32

a) Óbten el estimado de viajeros que suben en esta estación a las 10:30,

mediante interpolación lineal.

b) Realiza la estimación del apartado anterior mediante interpolación

cuádratica.

c) Utiliza el comando polyf it para hallar el polinomio de interpolación de

grado 3 y polyval para obtener la estimación deseada.

SOLUCIÓN

a) Para hallar el polinomio de grado 1,p1 (x) = a1 (x) + a2 , que nos de la estimación
de viajeros para x = 10, 5 tomamos los puntos de interpolación (10, 21)(11, 9),
y calculamos la recta que pasa por ellos:

p1 (10) = (10)a1 + a2 = 21

p1 (11) = (11)a1 + a2 = 9
obteniendo a1 = −12 y a1 = 141, por lo que el número de viajeros estiamdos a
las 10:30 será p1 (10, 5) = −12 ∗ 10, 5 + 141 = 15.
b) Hallaremos el polinomio de interpolacion cuadrático, p2 (x) = a1 x2 + a2 x + a3 , a
partir de los puntos (10,21),(11,9),y (12,11). El sistema de ecuaciones lineales
para este caso resulta:
    
100 10 1 a1  21
121 11 1 a2  =  9 
    
144 12 1 a3 11
    

Lo resolvemos mediante las siguientes instruciones en MATLAB:

66
 
1 x = 8:14;
2 y = [ 41 35 21 9 11 17 3 2 ] ;
3 x2 = x ( 3 : 5 ) ;
4 y2 = y ( 3 : 5 ) ;
5 A = vander ( x2 ) ;
6 p2 = A\ y2 ;
 

La parábola que interpola los puntos es p2 (x) = 7x2 − 159x + 911. Obtenemos
la estimación requerida con polyval (p2,10.5), resulta 13.25 viajeros.

c) Para hallar polinomios de interpolacion de mayor grado, vamos incorporan-

do los puntos por su proximidad a 10.5, donde se requiere la estimación.
Podrı́amos seguir el proceso anterior, construir la matriz de van der Monde
y resolver el ssitema asociado. Otra alternativa es usar el comando polyf it.
 
1 x3 = x ( 2 : 5 ) ;
2 y3 = y ( 2 : 5 ) ;
3 p3 = p o l y f i t ( x3 , y3 , 3 )
4
5 p3 =
6
7 0.0003
8 −0.0350
9 1.6167
10 16.0000
 

Con lo que obtenemos p3 (x) = 2x3 − 59x2 + 565x − 1729, y el número estimado
de viajeros es 14.

PREGUNTA 7.2

La solubilidad del cloruroámonico en el agua toma, a distintas temperatu-

ras, los valores dados en la siguiente tabla:

Temperatura 10 20 30 40 50 60
Solubilidad 33 37 42 46 51 55

a) Elige los 4 puntos centrales para hallar el plinomio de interpolación de

grado 3, obteniendo el mismo número de condicion del sistema de ecua-
ciones lineales asociado.

b) Óbten ahora el mismo polinomio de interpolación realizando un ndes-

plazamiento del origen al punto medio de las abscisas de inter4polación.
Comprueba la mejora del condicionamiento del sistema.

c) Representa el polinomio hallando junto a los puntos de interpolación y

utilizando para obtener una estimación de la solubilidad a las tempera-
turas 25,35, y 45.

SOLUCIÓN

67
a) Obtenemos al polinomio pedido ejecutando las siguientes instrucciones:
 
1 x = ([Link]) ’;
2 y = [33 37 42 46 51 5 5 ] ’ ;
3 x3 = x ( 2 : 5 ) ;
4 y3 = y ( 2 : 5 ) ;
5 A = vander ( x3 ) ;
6 ncon = cond (A)
7
8 ncon =
9
10 3.9175 e+06
11 p = A\ y3
12
13 p =
14
15 0.0003
16 −0.0350
17 1.6167
18 16.0000
 

El polinomio de grado 3 resulta:

p(x) = 0, 000333x3 − 0, 035x2 + 1, 61666x − 16

b) Veamos que el número de condición mejora, al utilizar otra expresión del po-
linomio, realizando un desplazamiento del origen.
 
1 xm = mean ( x3 ) ;
2 xd = x3 − xm ;
3 Ad = vander ( xd ) ;
4 ncon = cond (Ad)
5 pd = Ad\ y3
 

La expresión del polinomio desplazado será:

pd(x) = 0, 000333x − 353 + 0, 391667x − 35 − 44

El número de condición se ha reducido en 3 ordenes de magnitud, 3,822944x103 .

Notemos que, para evaluar el polinomio desplazado utilizando polyval, he-
mos de desplazar también los puntos donde se realiza la evaluación:
 
1 xg = l i n s p a c e ( x ( 1 ) , x ( end ) ) ;
2 yg = p o l y v a l ( pd , xg−xm ) ;
3 p l o t ( x , y , ’ * ’ , xg , yg )
 

68
 En la Figura 7.1 se ha representado el polinomio junto a los puntos de in-
terpolación. Teniendo en cuenta el desplazamiento de los datos, hallamos las
estimaciones de la solubilidad en las temperaturas perdidas, 25,35, y 45.
 
1 t = p o l y v a l ( pd , [ 2 5 , 3 5 , 4 5 ] −xm)
2 t =
3 39 ,7500 44 ,0000 48 ,2500
 

PREGUNTA 7.3

Las diferentes contracciones de un resorte dependiendo de las cargas apli-

cades vienen dadas en la tabla:

Carga(Kp) 5 10 15 20 25
Contracción(mm) 49 105 172 253 352

a) Escribe el polinomio de Lagrange de grado 2 para estimar la contracción

producida por una carga de 13 kp.

b) Edita una funcion en MATLAB para evaluar el polinomio de intyerpola-

ción producida de Lagrange en un punto dado a tomado como paráme-
tro de entrada, junto a los puntos de interpolación y ejeútala para obte-
ner el valor del polinomio anterior en x = 13

SOLUCIÓN

69
a) Elegimos los tres nodos más cercanos a x = 13, que son 10,15 y 20. El polino-
mio de Lagrange será de la forma:

(x − 10) (x − 20) (x − 10) (x − 20) (x − 10) (x − 20)

p2 (x) = 172× +105× +253×
(15 − 10) (15 − 20) (10 − 15) (10 − 20) (20 − 15) (20 − 10)

b) Utilizaremos como parámetros de entrada los puntos de interpolación y el

valor deonde se va a realizar la estimación:
 
1 f u n t i o n z =l a g r a n g e ( x , y , a )
2 m = length ( x ) ;
3 z = zeros ( size ( a ) ) ;
4 f o r i = 1 :m
5 Li = 1 ;
6 f o r j = 1 :m
7 if i = j
8 Li = Li * ( a−x ( j ) ) / ( x ( i ) − x ( j ) ) ;
9 end
10 end
11 z = z + y ( i ) * Li ;
12 end
 

Obtenemos la contracción del resorte para una masa de 13kp con la siguiente
llamada:
 
1 x = [10 15 2 0 ] ;
2 y = [105 172 2 5 3 ] ;
3 v = lagrange ( x , y , 1 3 )
 

PREGUNTA 7.4

a) Resuelve el sistema de ecuaciones lineales a partir del cual se obtiene

el polinomio de interpolación de Newton de grafo 3, que se ajusta a los
datos del 7.1

b) Óbten el polinomio de Newton de grado 3 contruyendo la tabla de di-

ferncias divididas, paso a paso.

c) Escribe el polinomio en su forma aislada para evaluarlo en x = 10, 5.

a) El poliniomio de Newton grado 3 para los punjtos más próximos a 10,5 y con-
siderando el orden de proximidad (10,21),(11,9),(12,11) y (9,35) tienen la si-
guiente forma: P3 (x) = c4 + c3 (x − 10) + 7(x − 10)(x − 11) + c1 (x − 10)(x − 11)(x − 12)
inponiendo las condiciones de interpolación obtendremos el siguiente sistema
de ecuacuiones lineales:

P3 (10)=21=c4
P3 (11)= 9 =c4 +c3
P3 (12)=11=c4 +c3 +2c2
P3 (9) =35=c4 −c3 +2c2 −6c1

70
 Al resolver obtenemos el polinomio de Newton: P3 (x) = 21 − 12(x − 10) + 7(x −
10)(x − 11) + 2(x − 10)(x − 11)(x − 12)
b) Construimos el polinomio sin tener que resolverlo un sitema lineal, partiendo
de los púntos de interpólación, construimos la tabla de diferencias divididas:

x y
10 21
9−12
11 9 11−10 = -12
11−9 2−(−12)
12 11 12−11 = 2 12−10 = 7
35−11 −8−2 5−7
9 35 9−12 = −8 9−11 = 5 9−10 =2

c) Expresamos el polinomio en su forma anidada:

P3 (x) = c4 + c3 (x − x1 ) + c2 (x − x1 )(x − x2 ) + c1 (x − x1 )(x − x2 )(x − x3 )
= c4 + (x − x1 )(c3 + (x − x2 )(c2 + (x − x3 )c1 ))
= ((c1 (x − x3 ) + c2 )(x − x2 ) + c3 )(x − x1 ) + c4

Evaluamos de dentro hacia afuera:

P3 (10, 5) = ((2(10, 5 − 12) + 7)(10, 5 − 11) − 12)(10, 5 − 10) + 12 = 14
El valor estimado de viajeros a los 10:30 horas,14, coincide con el obtenido a
esta misma hora con el polinomio normal de interpolación de grado 3 calcu-
lado en el apartado c) del problema 3.1.

PREGUNTA 7.5

Consideramos la función de Runge en el intervalo [−5, 5]

1
y=
1 + x2
Vamos a comprobar que cuando una función tiene caractierı́sticas muy di-
ferentes a las de los polinomios, la interpolaciónpolómica con nodos equies-
paciados puede resultar inadecuada.

a) Óbten el polinomio de interpolación de grado 4 que aproxima a la fun-

cion Range tomando nodos equiespaciados en el intervalo [−5, 5].

SOLUCIÓN
a) En primer lugar codificamos la función de Runge:
 
1 function y = f ( x )
2 y = 1./(1+ x ˆ2)
 

Calculamos el polinomio de interpolación de Newton de grado 4, tomando 5

nodos equidistantes en [−5, 5] y utilizando la tabla de difencias divididas:
 
1 x = linspace ( −5 ,5 ,5);
2 y = favel ( f , x ) ;
3 p = ta b l a d d ( x , y ) ;
 

71
 Los polinomios de Newton obtenidos son:

P0 (x) = 0, 038462 (7.1)

P1 (x) = P0 (x) + 0, 039788(x + 5) (7.2)
P2 (x) = P1 (x) + 0, 061008(x + 5)(x + 25)x (7.3)
P3 (x) = P2 (x) − 0, 026525(x + 5)(x + 2, 5)x(x − 2, 5) (7.4)

PREGUNTA 7.6

De un muelle con un extremo fijo colgamos distintas masas obteniendo l

siguente tabla:

Masa(kg) 5 10 15
Longitud(cm) 395 620 795

SOLUCIÓN

Para este ejecicio utilizamos la interpolación cuadrática para estimar la longitud

del muelle para una masa de 7 kg.
Imponiendo las condiciones de interpolación (5,395) , (10,620) y (15,795) teme-
nos un sistema con la siguiente expresión:
    
 25 5 1 a1  395
100 10 1 a2  = 620
    
225 15 1 a3 795
    

Creamos vectores de los puntos a interpolar:

 
1 x = [ 5 ; 10 ; 1 5 ] ;
2 y = [395 ; 620 ; 7 9 5 ] ;
3 A = [ vander ( x ) ] ;
4 p = A\ y
5 p =
6 −1.0000
7 60.0000
8 120.0000
 
Por lo que tenemos el polinomio de interpolación:

P2 (x) = −x2 + 60x + 120

Comprobamos que cumple las condiciones de interpolación:

 
1 polyval (p , x )
2
3 ans =
4 395.0000
5 620.0000
6 795.0000
 

72
 En m = 7kg:
 
1 v = polyval (p , 7 )
2 v =
3 491.0000
 
Representamos el poliniomio junto a los puntos de interpolación:
 
1 xn = 3 : . 1 : 2 0 ;
2 yn = p o l y v a l ( p , xn ) ;
3 p l o t ( x , y , ’ * ’ , xn , yn )
4 x l a b e l ( ’ Masas ’ )
5 y l a b e l ( ’ Longitudes ’ )
6 t i t l e ( ’ Polinomio de grado 2 ’ )
 

PREGUNTA 7.7

Un estudio sobre productividad en una fábrica ha determinado que el pro-

medio de piezas producidas por hora por sus trabajadores está relacionado
con el número de horas transcurridas a partir del comienzo de la jornada.
A partir de la siguienmte tabla haz una estimación del número de piezas
producidas a las 3,5 y 7 horas de trabajo.

Horas de trabajo 0 2 4 6 8
Piezas por hora 0 62 116 74 8

SOLUCIÓN

73
 Desplazamos del origen al valor medio de los nodos c=10.
 
1 x=[5;10;15];
2 y=[395;620;795];
3 A=vander ( x ) ;
4 nc=cond (A)
5
6 nc =
7
8 1.1079 e+03
9
10 p=A\ y
11
12 p =
13
14 −1.0000
15 60.0000
16 120.0000
17
18 c=mean ( x )
19
20 c =
21
22 10
23
24 xd=x−c ;
25 Ad=vander ( xd ) ;
26 ncd=cond (Ad)
27
28 ncd =
29
30 35.4120
31
32 pd=Ad\ y
33
34 pd =
35
36 −1
37 40
38 620
 
La expresión del polinomio sera:

P2 (x) = −x − 102 + 40x − 10 + 620

Como el polinomio de interpolación es único, al resolver el polinomio obtemos:

P2 (x) = −x2 + 60x + 120

74
 PREGUNTA 7.8

Hallar la cota del error al aproximar la función sen(x) con nodos equiespa-
ciados en [0, π] por el polinomio interpolador de tercer grado.

SOLUCIÓN

Según la expresión del error para este caso:

hn+1
   
|f (x) − Pn (x)| ≤ (n+1)!
máxx∈[0,π] f (n+1) (x)

Como la derivada n-ésima de al función está acotada por la unidad, podemos ga-
rantizar :

 
máxx∈[0,π] f (n+1) (x) ≤ 1

y decir que para cuaqluier valor de n el error sestá acotado por:

hn+1
 
|f (x) − Pn (x)| ≤ (n+1)!

1 π 4
 
en particular para n = 3 y h = π/4 será: E ≤ 24 4 = 0,0159

75
TRAZADORES CÚBICOS
CAPÍTULO 8
PREGUNTA 8.1

Escribe el sistema lineal que verifican los valores de la derivada segunda del
spline cúbico en los nodos, suponiendo conocidas éstas en los extremos.

SOLUCIÓN

Si rk es la derivada segunda del spline en el punto xk , se verifican las ecuación:

hk h + hk+1 h
rk + k rk+1 + k+1 rk+2 = ek+1 − ek (8.1)
6 3 6
k = 1, 2, ..., n − 1 (8.2)

donde:
yk+1 − yk
hk = xk+1 + xk y ek =
xk+1 − xk
sabemos que la derivada segunda de los extremos es conocida, pasamos al segun-
do mienbro el téremino en r y de la primera ecucación y el termino en rn + 1 de la
última, con lo que dichas ecuaciones quedan:

hn−1 h − hn h
rn−1 + n−1 rn = en − en−1 − n rn−1
6 3 6
Por tanto la matriz del sistema es:
 h +h h2 
 1 3 2 6 
 h
 2 h2 +h 3 h 3


 6 3 6 
 h3 h4 +h3 h4 
6 3 6
 
M =  
... ... ...


 

hn−2 hn−2 −hn−1 hn−1 
 

6 3 6
 
hn−1 hn−1 +hn 

6 3

76
 el termino independiente:

 e2 − e1 − h61 r1 
 
 

 e3 − e2 

 e4 − e3

g = 
 
.. 

 . 


 en − en−2


en − en−2 − h6n rn+1
 

La matriz es tridiagonal simetrica estrictamente diagonal dominante y dispersa.

El sistema Mr = g no necesita pivotación parcial par su resolución directa y es
apropı́ado mara los metodos iterativos.

PREGUNTA 8.2

Obtengamos estos vectores con MATLAB para la función de Runge en el

intervalo [-1,1]:
1
f (x) =
1 + 25x2

SOLUCIÓN
 
1 n=10; % Tomamos 10 s u b i n t e r v a l o s
2 h=2/n ; % Abscisas
3 x =[ −1:h : 1 ] ; % Ordenadas de Runge
4 y = 1 . / ( 1 + 25* x . ˆ 2 ) ; % y t i e n e n+1 c o m p o n e n t e s
5 a=y ( 1 : n )
6
7 a =
8
9 0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000
0.2000 0.1000 0.0588
10
11 b= d i f f ( y ) . / d i f f ( x )
12
13 b =
14
15 0.1018 0.2059 0.5000 1.5000 2.5000 −2.5000 −1.5000
−0.5000 −0.2059 −0.1018
 
Luego de obtener los coeficientes del spline, su evaluación en un punto x , elegir
los coeficientes ak , bk del polinomio correspondiente, qk (x) y evaluarlo según su
formula.
 
1 s =[b ’ , a ’ ]
2
3 s =
4
5 0.1018 0.0385
6 0.2059 0.0588
7 0.5000 0.1000

77
 8 1.5000 0.2000
9 2.5000 0.5000
10 −2.5000 1.0000
11 −1.5000 0.5000
12 −0.5000 0.2000
13 −0.2059 0.1000
14 −0.1018 0.0588
 
Para evaluar el splline, y los coeficientes, se necesitan los nodos. como ambos
datos: mkpp(comando, crea una estrucctura similar al spline)
 
1 ps=mkpp( x , s )
2 xg = ( − 1 : . 0 0 5 : 1 ) ’ ;
3 yg =1./(1+25* xg . ˆ 2 ) ;
4 ys=ppval ( ps , xg ) ;
5 p l o t ( x , y , ’ * ’ , xg , yg , xg , ys )
 

Comparando el spline lineal con el polinomio de grado 10 que interpola los mis-
mos nodos, observamos que este últimno presenta más oscilaciones, sobretodo en
los extremos del intervalo.

78
 PREGUNTA 8.3

Halla el spline cúbico que interpola la fuinción y = sen(x)+cos(x/2), toman-

do 7 nodos equiespaciados en el intervalo [−2π, 2π], subiendo que el valor
de la derivada segunda en ambos extremos del intervalo es 0.25.

SOLUCIÓN

Determinamos en primer lugar loos nodos y valores a interpolar

 
1 n=6;
2 x= l i n s p a c e ( −2* pi , 2 * pi , n+1)
3 x =
4
5 −6.2832 −4.1888 −2.0944 0 2.0944 4.1888 6.2832
6 y=s i n ( x )+ c o s ( x / 2 )
7 y =
8
9 −1.0000 0.3660 −0.3660 1.0000 1.3660 −1.3660 −1.0000
 
Los valores de la derivada segunda en los extremos son las componentes primera
y última del vector r:
r1 = 0,25; rn1 = 0,25. Obtendremos las restantes componetes r resolviendo el sis-
tema lineal que satisfacn las derivadas segundas.
Al pasar al segundo mienbro los terminos correspondientes de r1 y rn+1 queda un
sistema cuadrado de orden n − 1 con matriz tridiagonal simétrica.
Construimos su matriz por diagonales, operando adecuadamente con los elemen-
tos del vector h.
 
1 h= d i f f ( x ) ;
2 ds=h ( 2 : n − 1 ) / 6 ;
3 dl=h ( 1 : n − 1 ) / 3 ;
4 dr=h ( 2 : n ) / 3 ;
5 dp=dl+dr
6 M=diag ( dp)+ diag ( ds , −1)+ diag ( ds , 1 )
7
8 M =
9
10 1.3963 0.3491 0 0 0
11 0.3491 1.3963 0.3491 0 0
12 0 0.3491 1.3963 0.3491 0
13 0 0 0.3491 1.3963 0.3491
14 0 0 0 0.3491 1.3963
 
los terminos independientes:
 
1 e= d i f f ( y ) . / h
2
3 e =
4
5 0.6522 −0.3495 0.6522 0.1748 −1.3045 0.1748
6
7 g= d i f f ( e )

79
 8
9 g =
10
11 −1.0018 1.0018 −0.4775 −1.4792 1.4792
 
Si las derivadas no son nulas, hay que ajustar el primer y el último elemnto de g:
[-1.001757600061189 1.001757600061189 -0.477464829275686 -1.479222429336875
1.479222429336875]
 
1 g ( 1 ) = g (1) − h ( 1 ) * r1 /6
2 g=
3
4 −1.001757600061189 1.001757600061189 −0.477464829275686 −1.479222429336875 1
5 g ( n−1)=g ( n−1) −h ( n ) * rn1 /6
6 g=
7
8 −1.0890 1.0018 −0.4775 −1.4792 1.3920
 
Resolviendo el sistema lineal para obtener los valores desconocidos de las segun-
das derivadas, rk .
 
1 r=M\g ’ ;
2 r =[ r1 ; r ; rn1 ] ;
3
4 r =
5
6 0.2500
7 −1.0415
8 1.0460
9 −0.2727
10 −1.3232
11 1.3277
12 0.2500
 
Una vez obtenidas las deriavdas segundas, aplicamos la formula de los coeficien-
tes de los polinomios interpolantes:
 
1 r1=r ( 1 : n ) ;
2 r r =r ( 2 : n + 1 ) ;
3 a=y ( 1 : n ) ’ ;
4 b=e ’ − ( r r +2* r1 ) . * h ’ / 6 ;
5 c=r1 /2
6 d= d i f f ( r ) . / h ’ / 6 ;
 
Componemos la matriz del spline y representarlo graficamente.
 
1 s =[d ; c ; b ; a ] ;
2 ps=mkpp( x , s ) ;
3 xg= l i n s p a c e ( −2* pi , 2 * pi ) ;
4 ys=ppval ( ps , xg ) ;
5 yg=s i n ( xg )+ c o s ( xg / 2 ) ;
6 p l o t ( x , y , ’ * ’ , xg , yg , xg , ys , ’ −. ’ )
7
8 s =
9

80
10 −0.1028
11 0.1661
12 −0.1049
13 −0.0836
14 0.2110
15 −0.0858
16 0.1250
17 −0.5207
18 0.5230
19 −0.1363
20 −0.6616
21 0.6639
22 0.8412
23 0.0124
24 0.0172
25 0.8270
26 −0.8442
27 −0.8394
28 −1.0000
29 0.3660
30 −0.3660
31 1.0000
32 1.3660
33 −1.3660
 

En la gráfica observamos la aproxiamcion de la función mediante spline calcula-

do conociendo las derivadas segundas en los extremos.

81
 PREGUNTA 8.4

Las gráficas muestran distintas aproximaciones de la función de Runge.

1
y=
1 + 25x2
obtenidas interpolando sus valores en 9 nodods equidistantes en [0, 1]

SOLUCIÓN

Se obtiene con la función spline de MATLAB.

 
1 x = −1:.5:1; % Nodos
2 y =1./(1+25* x . ˆ 2 ) ; % Valores
3 xg = − 1 . 1 : 0 . 0 5 : 1 . 1 ; % A b s c i s a s para e l g r a f i c o
4 yg =1./(1+25* xg . ˆ 2 ) ; % Ordenadas
5 ys= s p l i n e ( x , y , xg ) ;
6 p l o t ( x , y , ’ * ’ , xg , yg , xg , ys )
 

Spline No-Nodo

La interpolación segmentaria lineal de esta función, a la que corresponde la se-

gunda gráfica:

82
 Spline lineal

PREGUNTA 8.5

Diseñe una curva ”suave”mediante interpolación degmentaria de manera

que sirva de guı́a para que el brazo robótico, partiendo de un punto deter-
minad, se mueva a una posición concreta para coger un objeto, lo deje en
otra pososicı́on y vuelva al punto de partida. Las tablas siguinetes dan el
codigo de mlas actuaciones del brazo y las posiciones qeu debe alacanzar.
Acci ón Significado
0 posici ón de partida
1 posici ón intermedia
2 posici ón de coger objeto
3 posici ón de soltar objeto
x y acci ón
0 0 0
2 4 1
6 4 1
7 6 2
12 7 1
15 1 3
8 −1 1
4 −2 1
0 0 0

SOLUCIÓN

83
 En priemr lugar reepresentamos gráficamente dichos puntos conectándolos me-
diante rectas, lo que equivale a la interpolación polinomica segmentaria lineal.[Link]

Obsevamos los mmovimentos bruscos, entonces es necesario una trayectoria que

permita al brazo moverse consuavidas de una posición a otra sin poiner en peligro
a los objetos que transporta.
Para diseñar el spline cúbico divideremos la trayectoria en tres trozos: desde la
posición de partida hasta coger el objeto, desde ese ı́nstante hasta que soltamos
el objeto y, por último, desde que se suelta el obejto hasta que se vuelve a la
posición original. Esta división es lógica puesto que el brazo ha de detenerse al
final de cada uan de estas trayectorias.
 
1 a c c i o n =[0 1 1 2 1 3 1 1 0 ] ;
2 c o g e r=f i n d ( a c c i o n ==2);
3 s o l t a r =f i n d ( a c c i o n ==3);
4 f i n =l e n g t h ( x ) ;
5 x1=x ( 1 : c o g e r ) , y1=y ( 1 : c o g e r )
6 x1 =
7
8 0 2 6 7
9
10
11 y1 =
12
13 0 4 4 6
14 x2=x ( c o g e r : s o l t a r ) , y2=y ( c o g e r : s o l t a r )
15
16 x2 =
17

84
18 7 12 15
19
20
21 y2 =
22
23 6 7 1
24 x3=x ( s o l t a r : f i n ) , y3=y ( s o l t a r : f i n ) ;
25
26 x3 =
27
28 15 8 4 0
29
30 y3 =
31
32 1 −1 −2 0
 
Para represntar los splines; obtenemos suficientes valores de x en cada tramo.
Calculamos la diostancia minioma entre los mpuntyos conicidos en toda la tra-
yuectoria. Ası́, el incremento en la variable x para definir el spline será la decima
parte de esta distancia mı́nima.
 
1 i n c r =min ( abs ( x ( 2 : f i n ) − x ( 1 : f i n − 1 ) ) ) / 1 0 ;
 
De este modo, al menos hábra 10 puntos de representación entre cada par de
puntos de la tabla anterior, distribuimos del siguente modo:
 
1 t 1= l i n s p a c e ( x ( 1 ) , x ( c o g e r ) ) ;
2 t 2= l i n s p a c e ( x ( c o g e r ) , x ( s o l t a r ) ) ;
3 t 3= l i n s p a c e ( x ( s o l t a r ) , x ( f i n ) ) ;
 
Calculamos los spline no-noddo mediante el comando de MATLAB:
 
1 ys1= s p l i n e ( x1 , y1 , t 1 ) ;
2 ys2= s p l i n e ( x2 , y2 , t 2 ) ;
3 ys3= s p l i n e ( x3 , y3 , t 3 ) ;
 
Y dibujamos la trayectoria resultante:
 
1 p l o t ( [ t 1 t 2 t 3 ] , [ ys1 ys2 ys3 ] , x , y , ’ o ’ )
2 xlabel ( ’ x ’ ) , ylabel ( ’y ’ ) , grid
3 t i t l e ( ’ R e c o r r i d o d e l brazo ’ ) ;
 

85
86
MODELOS LINEALES
CAPÍTULO 9
PREGUNTA 9.1

El paro registrado en febrero de los últimos años, según datos de la INEM

se resume en la siguiente tabla:

año número de parados

1991 2362373
1992 2337472
1993 2471228
1994 2774579
1995 2575873
1996 2427015
1997 2262721
1998 2067830
1999 1783940
2000 1659820
2001 1598920

a) Halla la recta, previa formación de las ecuaciones normales.

b) Ajusta un polinomios mı́nimo-cuadrático de grado 3 , Compara el

número de condición de las ecuaciones normales según se efectúe o
no el centrado de las abcisas.

c) Halla el error cuadrático y el ı́ndice de determinación del polino-

mio.

d) Estima el número de parados en febrero de 2002, de seguir esa

tendencia.

e) Obtén el mismo polinomio con POLIFIT, de tres maneras: a partir

de los datos iniciales, centrando las abscisas y centrando y escalando.
Conpara la magnitud relativa de los coeficientes. Genera gráficos con
los polinomios obtenidos de cada manera.

87
 SOLUCIÓN

a) El los coeficientes de la recta ajustada y = ax + b se obtienen con el siguiente

procedimiento:

Sean:
" P 2 P #
 1991 2362373
   
x x
 1992
   N= P P 0




 2337472 
 x x
 1993 2471228
    " P #
   xy
 1994  
2774579
 C= P
 




 y
 1995 2575873
    " #
   a
x =  1996 y =  2427015
   P=
 1997
  b




 2262721 

 1998 2067830
   
  
   
 1999   1783940 
   
 2000   1659820 
   
2001 1598920

Las ecuaciones normales son las siguientes:

NP = C
" P 2 P #" # " P #
a a a xy
P P 0 = P
a a b y

La solución de este sistema es:

P = N−1 C

El código en Matlab que solucione este sistema y la ecuación de la recta ajustada

se muetra a continuación:
Listing 9.1: Código en Matlab
1 x =[1991 , 1992 , 1993 , 1994 , 1995 , 1996 , 1997 , 1998 , 1999 , 2000
2 y=[2362373 ; 2337472 ; 2471228 ; 2774579 ; 2575873 ; 2427015 ; 2262721
3 plot ( a , b , ’ o ’ )
4 hold on
5

6 l =length ( a ) ;
7 N=[sum( a . ˆ 2 ) , sum( a ) ; sum( a ) , sum( a . ˆ 0 ) ] ;
8 C=[sum( a . * b ) ; sum( b ) ] ;
9 P=N\C ;
10 Q=P ’ ;
11 k=polyval (Q, x ) ;
12 plot ( x , k )

88
El resultado en la consola de P = N\C es P = 1.0e+08 *[ -0.0009 1.8941]
Por lo que la ecuación de la recta ajustada es la siguiente:
 
y = −0,0009e + 08x + 1,8941e + 08
 

La gráfica de la recta es la siguiente:

b) Para encontrar un polinomio de ajuste de grado 3 procedemos de manera simi-

lar al procedimiento anterior, pero en este caso, se usará un método más general,
donde:

N = AT A

Donde la matriz A para un ajuste de grado 3 con n = 11 tiene la siguiente forma:

x13 x12 x11 x10 

 

x23 x22 x21 x20 
 

x33 x32 x31 x30 
 
A = 
 .. .. .. .. 

 . . . . 
xn xn2
3 xn xn0
1 

Luego, la matriz N es:

x13 x12 x11 x10   P 6

 
 3 P 5 P 4 P 3
x23 x33 xn3
  
 x1 ···   x23 x22 x21 x20   P xi5 x
P i4
x
P i3
x
P i2
 
 2
 x x22 x32 ··· xn2 x31 x30  =  P xi4 x x x
  
x33 x32
 
N =  11 P i3 P i2 P i1
  
x21 x31 xn1
  
 x1 ···   .. .. .. ..   P xi x
P i2
x
P i1
x
P i0

 0
x1 x20 x30 ··· xn0
 
 . . . .  xi3 xi xi xi

xn xn2
3 1
xn xn 0
 

Luego, los coeficientes se obtienen al resolver el siguiente sistema:

P = N−1 C

89
   y1   P 3 
 
 3
 x1 x23 x33
··· xn3
  y   P xi yi 
 2
 x x22 x32
···   2   x2 y 
xn2
Donde C = AT y =  11   y3  =  P i i 
 x1 x21 x31
··· xn1  ..   P xi1 yi 
   
 0   .  
x1 x20 x30
··· xn0 xi0 yi

 
y1 1
Resolviendo el sistema mencionado en Matlab con el siguiente código se obtienen
 a 
 
 b 
los coeficientes P =   correspondientes al polinomio de ajuste:
 
 c 
 
d

y = ax3 + bx2 + cx + d

1 x =[1991 , 1992 , 1993 , 1994 , 1995 , 1996 , 1997 , 1998 , 1999 , 2000
2 y=[2362373 ; 2337472 ; 2471228 ; 2774579 ; 2575873 ; 2427015 ; 2262721
3 plot ( a , b , ’ o ’ )
4 hold on
5

6 l =length ( a ) ;
7 N=[sum( a . ˆ 2 ) , sum( a ) ; sum( a ) , sum( a . ˆ 0 ) ] ;
8 C=[sum( a . * b ) ; sum( b ) ] ;
9 P=N\C ;
10 Q=P ’ ;
11 k=polyval (Q, x ) ;
12 plot ( x , k )

La matriz de coeficiente que resulta es la siguiente:

1,0e + 11 ∗ [0,0000 − 0,00000,0075 − 5,2564]
Esta corresponde al siguiente polinomio:
 
P (x) = 1011 0,0000x3 − 0,0000x2 + 0,0075x − 5,2564

90
c)

(yi − P (xi ))2 , donde P (xi ) es el polinomio ob-

P
El error cuadrático está dado por:
[Link]:

Errorcuadratico = 1,5214e + 11

El ı́ndice de determinación de calcula de la siguiente manera:

(p(x ) − ȳ)2
P
I= P i
(yi − ȳ)2
Haciendo los cálculos en Matlab con la siguiente orden: indet=sum((k-mean(y)).ˆ2)/sum((y-m
se obtiene:

I = 0,8961

d)

Para estimar el número de parados para el 2002 , solo hace falta evaluar x = 2002
en el polonomio P (x). Esto usando el comando polyval de Matlab con el cual
se estima 1,1175e06 parados.
e)

Usando el comando polyfit(x,y,3) de Matlab se obtiene la matriz de coefi-

cientes que corresponden al polinomio:

P (x) = 1,0e13 (0,0000x3 − 0,0000x2 + 0,0039x − 2,6221)

91
 PREGUNTA 9.2

La densidad del aire (en g/m3) varı́a con la altura h según la tabla siguiente:

h 7 10 21 27 34 39 43 47 51 55 59 59 61
p 556 369 191 75 26.2 9.9 4.4 2.3 1.4 0.8 0.5 0.33 0.25

a) Halla el ajuste polinómico de menor grado con indice de determi-

nación mayor que 0.99. Cual es el error cuadrático?

b) Dibuja los datos y el polinomio ajustado. Explica por qué este ajuste
no resulta adecuado para alturas grandes.

c) Ajusta los datos a un modelo exponencial. Compara el error y el

comportamiento del ajuste con el caso anterior

SOLUCIÓN

Hallando polinomios sucesivos con el procedimiento realizado en el problema

anterior y hallando el ı́ndice de determinación denotado por:

(p(x ) − ȳ)2
P
I= P i
(yi − ȳ)2

El código en Matlab usado para hallarlos es similar al anterior:

1 x =[7 , 10 , 21 , 27 , 34 , 39 , 43 , 47 , 51 , 55 , 59 , 59 , 6 1 ] ;
2 y =[556 ; 369 ; 191 ; 75 ; 2 6 . 2 ; 9 . 9 ; 4 . 4 ; 2 . 3 ; 1 . 4 ; 0 . 8 ; 0 . 5 ; 0 . 3
3 plot ( x ’ , y , ’ o ’ )
4 hold on
5 l =length ( x ) ;
6 A=[ x ’ . ˆ 5 , x ’ . ˆ 4 , x ’ . ˆ 3 , x ’ . ˆ 2 , x ’ . ˆ 1 , x ’ . ˆ 0 ] ;
7 N=A’ *A ;
8 C=A’ * y ;
9 P=N\C ;
10 Q=P ;
11 k=polyval (Q, x ) ;
12 plot ( x , k )
13 disp ( P )
14

15 i n d e t=sum ( ( k−mean ( y ) ) . ˆ 2 ) / sum ( ( y−mean ( y ) ) . ˆ 2 ) ;

16 disp ( i n d e t )

El número de columnas de la matriz A en la lı́nes 6 del código variará según el

grado del polinomio a ajustar.

92
 Polinomio de grado 1 Polinomio de grado 2
I = 0.7329 I = 0.9704

Polinomio de grado 3 Polinomio de grado 4

I = 0.9872 I = 0.9879

Polinomio de grado 5 Polinomio de grado 6

I = 0.9915 I = 0.9968
El polinomio de 5 y 6to grado tiemen un ı́ndice de determinación mayor que 0.99.

93
b)

La siguiente gráfica muestra polinomios de 5,6, 7 y 8vo grado, estos no tienden

a seguir la tendencia de los puntos, por lo que no es apropiado para realizar
pronósticos ya que solo se ajustan a los pares de datos.

c)
Para realizar un ajuste exponencial se hace la siguiente transformación para ha-
llar una función de la forma y = kebx :

y = kebx
ln y = ln(kebx )
ln y = ln k + bx, ln k = a, ln y = z
z = a + bx

El problema se reduce a ajustar un polinomio de grado 2 (recta) de la forma

z = a + bx , para hacerlo solo hace falta cambiar la matriz y por otra z = ln y,
para ello, se procede a calcularlo con el siguiente código de Matlab:

94
 1 x =[7 , 10 , 21 , 27 , 34 , 39 , 43 , 47 , 51 , 55 , 59 , 59 , 6 1 ] ;
2 y =[556 ; 369 ; 191 ; 75 ; 2 6 . 2 ; 9 . 9 ; 4 . 4 ; 2 . 3 ; 1 . 4 ; 0 . 8 ; 0 . 5 ; 0 . 3
3 l y=log ( y ) ;
4 p l o t ( x ’ , y , ’ o ’ , ’ LineWidth ’ , 2 )
5

6 hold on
7

8 A=[ x ’ . ˆ 1 , x ’ . ˆ 0 ] ;
9 N=A’ *A ;
10 C=A’ * l y ;
11 P=N\C ;
12

13 d=l i n s p a c e ( 0 , 7 0 ) ;
14 r=polyval ( P , d ) ;
15 exd=exp ( P ( 1 , 1 ) . * d+P ( 2 , 1 ) ) ;
16 p l o t ( d , exd )
17 %p l o t ( d , r , ’ LineWidth ’ , 2)
18 disp ( P )

los coeficientes de la recta z = a + bx son a = −0,1462 y b = 7,8142, pero como

ln y = z, la ecuación queda ası́:

z = a + bx
z = −0,1462 + 7,8142x
ln y = −0,1462 + 7,8142x

y = e−0,1462+7,8142x
 

La gráfica de esta última función y de los pares de puntos se muestra a continua-

ción. Se puede observar que el ajuste es adecuado ya que sigue la tendencia de
los puntos.

95
 PREGUNTA 9.3

La tabla siguiente muestra los valores de la temperatura ambiente medida

a distintas horas del dı́a

Hora 6 8 10 12 14 16 18 20
Grados 7 9 12 18 21 19 15 10

Normaliza las abscisas y ajusta un polinomio minimo-cuadrático de cuarto

grado por los siguientes procedimientos:

a) Planteando un sistema lineal sobredeterminado y resolviendo di-

rectamente las ecuaciones normales.

b) Aplicando previamente la factorización QR según el algoritmo de

Gram- Scmidt a la matriz del sistema.

c) Aplicando el algoritmo qr de MATLAB para obtener la factoriza-

ción.

d) Resolviendo el sistema sobredeterminado mediante la contraba-

rra,.

e) Calculando directamente el polinomio mediante polyfit

SOLUCIÓN

a)

Resolveremos direcamente el sistema siguiente:

NP = C

conde la matriz de coeficientes es:

P = N−1 C

Además,

N = AT A

Donde
x13 x12 x11 x10 
 

x23 x22 x21 x20 
 

x33 x32 x31 x30 
 
A = 
 .. .. .. .. 

 . . . . 
xn xn2
3 xn xn0
1 

96
 Para hallar la matriz de coeficientes, se usará el siguiente código de Matlab:
1 h=[6 ,8 ,10 ,12 ,14 ,16 ,18 ,20] ’;
2 t =[7 ,9 ,12 ,18 ,21 ,19 ,15 ,10] ’;
3 plot ( h , t , ’ o ’ , ’ linewidth ’ , 1 . 5 )
4 hold on
5 A=[h . ˆ 4 , h . ˆ 3 , h . ˆ 2 , h . ˆ 1 , h . ˆ 0 ] ;
6 N=A’ *A ;
7 C=A’ * t ;
8 P=inv (N) *C ;
9 Q=P ’ ;
10 y f=polyval (Q, h ) ;
11

12 d=l i n s p a c e ( 3 , 2 1 , 1 0 0 )
13 k=polyval (Q, d ) ;
14 plot (d , k , ’ linewidth ’ ,1)
15 disp ( P )
16

17 P2=polyval ( P , h ) ;
18 e r r =sum ( ( k− t ) . ˆ 2 ) ;
19 disp ( e r r )
20 %i n d i c e de d e t e r m i n a c i o n
21 i n d e t=sum ( ( yf −mean ( t ) ) . ˆ 2 ) / sum ( ( t −mean ( t ) ) . ˆ 2 ) ;
22 disp ( i n d e t )

La gráfica generada se muestra a continuació[Link] ı́ndice de determinación es

0.9871

b)

c)

97
d)

PREGUNTA 9.4

La tabla adjunta muestra valores de una señal sinusoidal

y = K sin(ωt + φ)

de la que se conoce la frecuencia ω = 1 y se quiere estimar la amplitud K y

la fase φ. El cambio

P = K cos φ, Q = K sin φ

convierte el problema en uno lineal:

y = P sin(ωt) + Q cos(ωt)

t 0 1.5 3 4.5 6
y -0.0713 -2.1641 0.8403 2.2658 -0.5083

Estima los valores desconocidos, halla el eror cuadrático y representa los

datos y el ajuste.

98
PREGUNTA 9.5

Año Abonados
1988 11689
1989 29783
1990 54712
1991 108451
1992 180296
1993 257250
1994 411930
1995 928955

Ajusta el polinomio de grado 7 a los datos de la tbla anterior. Modifi-

ca previamente el origen de tiempos para mejorar el comportamiento
numérico del problema. Representa gráficamente el polinomio obte-
nido.

Elige el menor grado de polinomio mı́nimo cuadrático que aproxime

estos datos con la ocndición de que el ı́ndice de determinación sea
mayor que 99 %

Compara las predicciones efectuadas mediante los polinomios ante-

riores con los datos de la tabla siguiente

Año miles de abonados

1996 2988
1997 4330
1998 7051
1999 14884
2000 24266
2001 29656

Con toda la información de que dispones, calcula un informe para

una operadora de telefonı́a móvil estimando la evolución del merca-
do en 2002 y 2003. Comprueba tus predicciones con los datos reales
publicados.

99
POLINOMIOS ORTOGONALES
CAPÍTULO 10
PREGUNTA 10.1

Aproxima la función f (x) = ex en el intervalo [0, 3] mediante un polinomio

de quinto grado. Halla el error cuadrático.

SOLUCIÓN

Se sabe que si {φn (x)} es un conjunto infinito de funciones ortogonales sobre un

intervalo [a, b], entonces una función f (x) definida sobre este intervalo, entonces
es posible determinar un cpnjunto de de coeficientes cn , n = 0, 1, 2, · · · para el que:

f (x) = c0 φ0 + c1 φ1 + c2 φ2 + · · ·

Si se multiplica esta expresiónon por φm y luego se integra se tiene:

Z b Z b Z b Z b
f (x)φm (x) dx = c0 φ0 φm (x) dx + c1 φ1 φm (x) dx + c2 φ2 φm (x) dx + · · ·
a a a a

Donde φn yφm pertenecen al mismo conjunto ortogonal. Entonces cada integral

del lado derecho es 0 , pues si dos funciones f (x) y g(x) son tortogonales, enton-
ces:
Zb


f (x), h(x) = f (x)g(x) dx
Rb a
0 = a f (x)g(x) dx

excepto cuando m = n, entonces se tiene:

Z b Z b Z b
f (x)φn (x) dx = cn φn (x)φm (x) dx = cn φn2 dx
a a a

Donde a0 φ0 +a1 φ1 +a2 φ2 +a3 φ3 +a4 φ4 +a5 φ5 Es el polinomio de grado 5 a ajustar.

Si hacemos que:

f (x) = ex = a0 φ0 + a1 φ1 + a2 φ2 + a3 φ3 + a4 φ4 + a5 φ5 =
a0 x0 + a1 x1 + a2 x2 + a3 x3 + a4 x4 + a5 x5 + q(x)

100
O sea, el conjunto ortogonal esel siguiente:
n o
{φ0 , φ1 , φ2 , φ3 , φ4 , φ5 } = x0 , x1 , x2 , x3 , x4 , x5

Para aproximar una función de 5to grado a la función ex , se usarán polinomios

ortogonales, donde al minimizar el error entre ambas curvas nos generará el po-
linomio pedido.
Z3
2
E = |ex − (a0 φ0 + a1 φ1 + a2 φ2 + a3 φ3 + a4 φ4 + a5 φ5 )|2 dx
0

Donde q(x) es una función ortogonal a a0 x0 + a1 x1 + a2 x2 + a3 x3 + a4 x4 + a5 x5 , y de

la defición de producto escalar y la propiedad de polinomios ortogonales:

 Rb
f (x), h(x) = a f (x)g(x) dx
Rb
0 = a f (x)g(x) dx

Si se multiplica por cada φi y se hace que q(x) sea ortogonal a φ0 , φ1 , φ2 , φ3 , φ4

, φ5 y φ6 . Entonces se obtiene el siguiente sistema:








f (x), φ0 = a0 φ0 , φ0 + a1 φ1 , φ0 + a2 φ2 , φ0 + a3 φ3 , φ0 + a4 φ4 , φ0 + a5 φ5 , φ0








f (x), φ1 = a0 φ0 , φ1 + a1 φ1 , φ1 + a2 φ2 , φ1 + a3 φ3 , φ1 + a4 φ4 , φ1 + a5 φ5 , φ1








f (x), φ2 = a0 φ0 , φ2 + a1 φ1 , φ2 + a2 φ2 , φ2 + a3 φ3 , φ2 + a4 φ4 , φ2 + a5 φ5 , φ2








f (x), φ3 = a0 φ0 , φ3 + a1 φ1 , φ3 + a2 φ2 , φ3 + a3 φ3 , φ3 + a4 φ4 , φ3 + a5 φ5 , φ3








f (x), φ4 = a0 φ0 , φ4 + a1 φ1 , φ4 + a2 φ2 , φ4 + a3 φ3 , φ4 + a4 φ4 , φ4 + a5 φ5 , φ4








f (x), φ5 = a0 φ0 , φ5 + a1 φ1 , φ5 + a2 φ2 , φ5 + a3 φ3 , φ5 + a4 φ4 , φ5 + a5 φ5 , φ5

Que esequivalente a:

D E D E D E D E
hex , 1i = a0 h1, 1i + a1 hx, 1i + a2 x2 , 1 + a3 x3 , 1 + a4 x4 , 1 + a5 x5 , 1
D E D E D E D E
hex , xi = a0 h1, xi + a1 hx, xi + a2 x2 , x + a3 x3 , x + a4 x4 , x + a5 x5 , x
D E D E D E D E D E D E D E
ex , x2 = a0 1, x2 + a1 x, x2 + a2 x2 , x2 + a3 x3 , x2 + a4 x4 , x2 + a5 x5 , x2
D E D E D E D E D E D E D E
ex , x3 = a0 1, x3 + a1 x, x3 + a2 x2 , x3 + a3 x3 , x3 + a4 x4 , x3 + a5 x5 , x3
D E D E D E D E D E D E D E
ex , x4 = a0 1, x4 + a1 x, x4 + a2 x2 , x4 + a3 x3 , x4 + a4 x4 , x4 + a5 x5 , x4
D E D E D E D E D E D E D E
ex , x5 = a0 1, x5 + a1 x, x5 + a2 x2 , x5 + a3 x3 , x5 + a4 x4 , x5 + a5 x5 , x5

El sistema anterior puede reducirse a una matriz N de coeficientes y otra A de

términos independientes

Na = A

Los coeficientes buscados se obtienen con:

a = N−1 A

101
Donde:
D E Zb
i−1 j−1
Nij = x , x = xi−1 xj−1 dx
a
D E Zb
Ai = f (x), xi−1 = f (x)xi−1 dx
a

Al calcular cada Nij y Bi se ontienen la siguientes matrices:

9 81 243 243
 
 3 9 
2 4 5 2
 
 
 9 81 243 243 2187 

 2 9 
4 5 2 7 
81 243 243 2187 6561
 
 9
 

N =  81
 4 5 2 7 8 
243 243 2187 6561 
2187
 
 4 5 2 7 8
 

 243 243 2187 6561 59049 
 2187 
 5 2 7 8 10
 
 243 2187 6561 59049 177147


 2187 
2 7 8 10 11

e3 − 1
 
 
 2e3 + 1
 

 5e3 − 2
 

A = 
 12e3 + 6



 33e3 − 24
 

87e3 + 120
 

Como a = N−1 A , reemplazando estas últimas matrices y resolviendo se obtiene:

a0 0,9953
   
   

 a1  
  1,064 

a2 0,297
   
a =   = 
   

 a3   0,419 
  −9,859 × 10−2
  
a4
 

  
a5 4,07 × 10−2


Donde el polinomio ajustado es el siguiente:

P5 (x) = 0,9953 + 1,064x + 0,297x2 + 0,419x3 + −9,859 × 10−2 x4 + 4,07 × 10−2 x5

102
 PREGUNTA 10.2
Z 1  2
Halla a, b, c y d para que la integral ex − (a + bx + cx2 + dx3 ) dx sea
−1
mı́nima

SOLUCIÓN

Hallar los coeficientes a, b, c y d que hace que la integral mostrada sea mı́nima im-
plica hallar la aproximación de un polinomio de tercer grado mı́nimo cuadrático
que se ajuste a la función f (x) = ex en [−1, 1].
Del ejemplo anterior se tiene la siguiente matriz:

Nc = A → c = N−1 A

103
Donde:
D E Zb
i−1 j−1
Nij = x , x = xi−1 xj−1 dx
a
D E Zb
Ai = f (x), xi−1 = f (x)xi−1 dx
a

Integrando cada elemento de la matriz se obtiene:

 2,0000 −0,0000 0,6667 0,0000 

 
 −0,0000 0,6667 0 0,4000 
N = 

 0,6667 0 0,4000 −0,0000 

 
0,0000 0,4000 −0,0000 0,2857

 2,3504 
 
 0,7358 
A = 
 
 0,8789 

 
0,4495

 a 
 
 b 
Entonces, la matriz c =   es:
 
 c 
 
d

 a   0,9963
   

 b   0,9980 
c =   = 
   
 c   0,5367


   
d 0,1761

Por lo tanto, el polinomio que hace mı́nima a la integral del enunciado es:

0,9963 + 0,9980x + 0,5367x2 + 0,1761x3

104
 El programa usado para resolver este ejercicio es el siguiente:
1 N=z e r os ( 4 , 4 ) ;
2 f o r a =1:4
3 f o r b=1:4
4 orN = @( x ) x . ˆ ( a − 1 ) . * x . ˆ ( b − 1 ) ;
5 n i= i n t e g r a l ( orN , − 1 , 1 ) ;
6 N( a , b)= n i ;
7 end
8 end
9 disp (N)
10

11 A=z e r os ( 4 , 1 ) ;
12 f = @( x ) exp ( x ) ;
13 f o r k=1:4
14 C = @( x ) f ( x ) . * x . ˆ ( k − 1 ) ;
15 a i = i n t e g r a l (C, − 1 , 1 ) ;
16 A( k , 1 ) = a i
17 end
18 disp (A)
19 Pol=N\A
20 Dom=l i n s p a c e ( − 3 , 3 , 1 0 0 ) ;
21 Ran=polyval ( flipud ( Pol ) ,Dom ) ;
22 p l o t (Dom, f (Dom) , ’ y ’ ,Dom, Ran , ’ k ’ , ’ l i n e w i d t h ’ , 2 )

Cuya gráfica generada es la siguiente:

105
 PREGUNTA 10.3

Aproxima f (x) = |x| en [−1, 1] por un polinomio de grado 4, utilizando po-

linomios de Legendre.

SOLUCIÓN

Generando los 5 primeros polinomios de Legendre por la fórmula de Rodrigues:

1 dn h 2 n
i
Pn (x) = (x − 1)
2n n! dxn
Se tiene:

P0 (x) = 1
P1 (x) = x
3 1
P2 (x) = x2 −
2 2
5 3 3
P3 (x) = x − x
2 2
35 4 15 2 3
P4 (x) = x − x +
8 4 8
Estos polinomios son ortogonales respecto al pero w(x) = 1

Z1 

 0 ,i , j
Pi (x)Pj (x) dx = 
 2
 Cj = 2j + 1 , i = j

−1 

Aproximamos la función f (x) = |x| en el intervalo [−1, 1] a la combinación lineal

P (x) de polinomios de Legendre

P (x) = a0 P0 (x) + a1 P1 (x)a2 P2 (x) + a3 P3 (x) + a4 P4 (x)

Los coeficientes ak se calculan del siguiente modo

Z1
1
ak = f (x)Pk (x) dx
Ck −1
2
Ck = , k = 0, 1, 2, 3, · · · , n
2j + 1

Integrando, calculamos los coeficientes a0 , a1 , a2 , a3 ya4

106
 Z 1
1 1
a0 = (1) |x| dx =
2 −1 2
Z 1
3
a1 = (x) |x| dx = 0
2 −1
Z 1
5 3 2 1 5

a2 = x − |x| dx =
2 −1 2 2 8
Z 1
7 5 3 3

a3 = x − |x| dx = 0
2 −1 2 2
Z 1
9 35 4 15 2 3 3

a4 = x − x + |x| dx = −
2 −1 8 4 8 16

Por lo que el polinomio ajustado P (x) es el siguiente:

1 5 3
P (x) = P0 (x) + P2 (x) − P4 (x)
2 8 16
Graficando la función en Matlab usando el siguiente código se obtiene la siguien-
te figura:
1

2 x=l i n s p a c e ( − 1 . 5 , 1 . 5 , 1 0 0 ) ;
3 P0 = @( x ) 1 ;
4 P2 = @( x ) ( 3 / 2 ) . * x . ˆ 2 − 1 / 2 ;
5 P4 = @( x ) 1 / 8 . * ( 3 5 . * x . ˆ 4 − 3 0 . * x . ˆ 2 + 3 ) ;
6 f = @( x ) abs ( x ) ;
7 PA = @( x ) ( 1 / 2 ) . * P0 ( x ) + ( 5 / 8 ) . * P2 ( x ) − ( 3 / 1 6 ) . * P4 ( x ) ;
8 plot ( x , f ( x ) , ’k ’ , ’ linewidth ’ ,2)
9 hold on
10 p l o t ( x , PA( x ) , ’ l i n e w i d t h ’ , 2 )

107
108
 PREGUNTA 10.4

Verifica la fórmula de Rodrigues para los primeros polinomios de Legendre:

1 dn 2
Pn (x) = (x − 1)n
2n n! dxn

SOLUCIÓN

Evaluando directamente la fórmula de Rodrigues, se obtienen los diez primeros

polinomios de Legendre:

P0 (x) = 1
P1 (x) = x
3 1
P2 (x) = x2 −
2 2
5 3 3
P3 (x) = x − x
2 2
35 4 15 2 3
P4 (x) = x − x +
8 4 8
63 5 −35 3 15
P5 (x) = ·x + ·x + ·x
8 4 8
231 6 −315 4 105 2 −5
P6 (x) = ·x + ·x + ·x +
16 16 16 16
429 7 −693 5 315 3 −35
P7 (x) = ·x + ·x + ·x + ·x
16 16 16 16
6435 8 −3003 6 3465 4 −315 2 35
P8 (x) = ·x + ·x + ·x + ·x +
128 32 64 32 128
12155 9 −6435 7 9009 5 −1155 3 315
P9 (x) = ·x + ·x + ·x + ·x + ·x
128 32 64 32 128
46189 10 −109395 8 45045 6 −15015 4 3465 2 63
P1 0(x) = ·x + ·x + ·x + ·x + ·x −
256 256 128 128 256 256

109
APROXIMACIÓN TRIGONOMÉTRICA
CAPÍTULO 11
PREGUNTA 11.1

Dada la ecuación diferencial resuelva el P.V.I:

solucion

grafica

110
 PREGUNTA 11.2

integrar la ecuación diferencial de primer orden que describe la carga de

un condensador.

dq Vε q
= −
dt R RC
solucion

grafica

111
 PREGUNTA 11.3

Comprobar que la ED

(2x + y)dx + (x − 3y)dy = 0

es exacta. Resolver dicha ecuacion y representar las soluciones

solucion

grafica

112
113
ECUACIONES DIFERENCIALES DE PRIMER ORDEN
CAPÍTULO 12
PREGUNTA 12.1

Un proyectil de masa m=0.11 kg se lanza verticalmente hacia arriba con

una velocidad inicial de 8m/s siendo frenado por la accion de la gravedad
y por la resistencia del aire. Se satisface la ecuacion diferencial

mv 0 = −mg − k.v. |v|

donde g = 9,8m/s2 y k = 0,002. Hallar la velocidad al cabo de 2 segundos

usando el metodo Euler, Heun, Euler modificado y Runge-Kutta

12.1.1 Metodo de Euler

12.1.2 Metodo de Euler Modificado

114
Metodo de Euler Modificado con n=20

115
grafica

116
Metodo de Euler Modificado con n=10

117
 PREGUNTA 12.2

Un liquido de baja viscosidad fluye de un tanque conico invertido por un

orificio circular a razon de
√
dx 2
p x
= −0,6πr 2g
dt x
donde r=3 mm es el radio del orificio, x es la altura de la superficie del
liquido desde el vertice del cono y A(x) es el area transversal del deposito
a una altura x. Supongamos que la altura inicial del liquido es de 8 m y el
volumen inicial de 60πm3 . Hallar el nivel del liquido al cabo de una hora,
mediante

118
119
ANEXOS

120