Métodos Numéricos en Los Lenguajes Matlab y Microsoft Visual C
Métodos Numéricos en Los Lenguajes Matlab y Microsoft Visual C
Métodos Numéricos en Los Lenguajes Matlab y Microsoft Visual C
EN LOS LENGUAJES
MATLAB Y MICROSOFT
VISUAL C#.NET
Prologo
Este libro ha sido diseado para los estudiantes y profesionales, especialmente de las reas de la Ingeniera, Economa, Ciencias Administrativas,
Ciencias Exactas, entre otras que tengan conocimientos bsicos de Matlab
y/o slidos en Microsoft Visual C#.NET.
Para el caso de Microsoft Visual C#.NET., se ha puesto nfasis en la reutilizacin de cdigo con el empleo de una librera comn de clases.
No pretende cubrir temas nuevos que vienen con Microsoft Visual Studio .NET como:
Libreras .NET (que trabajan con Framework)
ADO.NET en Visual C # .NET.
Aplicaciones Windows Form.NET
Bibliotecas de controles.NET
Servicio de Windows.NET
Servicio de Web.NET
Debug and Release.NET
Pues se considera que es materia de otra obra con temas avanzados.
Esta obra, se encuentra formada por ocho captulos, y dos apndices,
estructurados de la siguiente manera:
En el captulo I, se consideran temas relacionados con teora y propagacin de errores, iniciando con definiciones bsicas, se pasa a ilustrar con
ejemplos para luego complementar con programas en Matlab y Microsoft
Visual C#.NET mostrando los resultados de la ejecucin de los programas.
En el captulo II, se describen mtodos grficos y numricos para calcular races reales de funciones trascendentes por los mtodos de la Biseccin,
Aproximaciones Sucesivas, Aproximaciones Sucesivas Modificado, Rgula
- Falsi, Secante, Newton - Raphson y Muller. En la siguiente seccin se presenta el Mtodo de Horner para calcular races reales de funciones. En la
ltima seccin de este captulo se presenta el Mtodo de Newton - Bairstow
para calcular races reales y complejas de funciones polinmicas de la forma
F(x) = 0. Todos los mtodos se describen con un anlisis matemtico, algoritmo en seudocdigo, ejemplos ilustrativos e implementacin en Matlab y
Microsoft Visual C#.NET mostrando los resultados de la ejecucin de los
programas.
En el captulo III, se desarrollan mtodos para resolver sistemas de ecuaciones lineales y simultneas por Gauss, Gauss Jordan, Gauss Seidel y
Crout. En la segunda seccin se propone un mtodo para calcular la inversa
de una matriz a travs de la aplicacin de la resolucin de una secuencia
de sistemas de ecuaciones lineales y simultneas por Crout. En la tercera
seccin se presenta una variante del Mtodo de Crout para resolver sistemas simtricos de ecuaciones lineales y simultneas. En la cuarta seccin
se presenta otra variante del Mtodo de Crout para resolver sistemas simtricos en banda de ecuaciones lineales y simultneas. En la seccin final de
este captulo se presenta otra variante del Mtodo de Crout para resolver
sistemas simtricos Sky - Line de ecuaciones lineales y simultneas. Todos
los mtodos se describen con un anlisis matemtico, algoritmo en seudocdigo, ejemplos ilustrativos e implementacin en Matlab y Microsoft Visual
C#.NET mostrando los resultados de la ejecucin de los programas.
En el captulo IV, se describen algoritmos para interpolar por los mtodos de Interpolacin Lineal, Polinomios de Lagrange y Newton mediante
diferencias divididas. Todos los mtodos se describen con un anlisis matemtico, algoritmo en seudocdigo, ejemplos ilustrativos e implementacin
en Matlab y Microsoft Visual C#.NET mostrando los resultados de la ejecucin de los programas.
En el captulo V, se aplican los algoritmos descritos en los captulos anteriores para ajustar una muestra estadstica a una funcin polinmica mediante los mtodos de los Mnimos Cuadrados; Polinomios de Lagrange y
Newton, mostrando el diagrama de dispersin y el grfico de la familia de
polinomios procesados. Los mtodos se describen con un anlisis matemtico, algoritmo en seudocdigo, ejemplos ilustrativos e implementacin en
Matlab y Microsoft Visual C#.NET mostrando los resultados de la ejecucin
de los programas.
En el captulo VI, se describen mtodos para calcular derivadas de primer orden y de orden superior en un punto dado, especificando la longitud
1
NMEROS APROXIMADOS Y ERRORES
CAPITULO I
NUMEROS APROXIMADOS Y ERRORES
Generalidades
Nmero Aproximado .- Se dice que un nmero es aproximado cuando difiere
ligeramente del nmero exacto.
Sean: R = nmero exacto y
r = nmero aproximado de R
Entonces se presentan dos casos:
a) cuando r > R, r es aproximado por exceso, y
b) cuando r < R, r es aproximado por defecto.
Error de un nmero aproximado.-Es la diferencia que existe entre el nmero
exacto y el nmero aproximado, es decir:
Error de r = Er = R - r
(1.1).
Error absoluto de un nmero aproximado.-Es igual al valor absoluto del error de
un nmero aproximado.
Matemticamente se define as: Error absoluto de r = EAr = R - r
(1.2).
Existen dos casos en los nmeros exactos:
a) R es conocido.
b) R es desconocido.
Para el caso (a), el error absoluto de un nmero aproximado, se lo puede determinar
aplicando la frmula (1.2).
En la prctica se presenta el caso (b), por consiguiente debe introducirse el concepto
de un error superior estimado, denominado:
Cota de error absoluto.-Que se define como un nmero cualquiera Cr no menor que
el error absoluto de un nmero aproximado,
Entonces Cr EAr
(1.3).
De (1.2) se puede deducir que: R = r Er
(1.4).
Ejemplos 1.1:
a) si se conoce que A = 2 tiene una variacin de 1.41 A 1.42 entonces
Ea = 1.42 - 1.41 = 0.01
b) si se sabe que A = 2 tiene una variacin de 1.411 A 1.414 entonces
Ea = 1.414 - 1.411 = 0.003
Por consiguiente, se puede tener un nmero infinito de cotas de error absoluto, se
debe seleccionar entonces el menor de los valores, segn las circunstancias.
Cuando se escribe un nmero aproximado, procedente de una medida, es comn dar
su cota de error absoluto, por ejemplo si la longitud de un segmento es l = 235 cm
con un error de 0.5 cm, se escribe L = 235 cm 0.5 cm. En este caso la cota de error
absoluto es Cl = 0.5 cm y la magnitud exacta de la longitud del segmento cae dentro
del margen 234.5 cm L 235.5 cm
La cota de error absoluto no indica las caractersticas de la calidad de una medida;
suponiendo que al medir las longitudes de dos segmentos se tiene L1 = 80.8 cm
0.2 cm y L2 = 5.2 cm
0.2 cm, independientemente de que los dos errores coincidan, la primera es mejor
que la segunda medida, por consiguiente, debe introducirse el criterio de:
Error relativo.-Que es la relacin entre el error absoluto y el mdulo (valor absoluto)
del nmero exacto: r = EAr / R (1.5).
Aplicando el mismo criterio de la cota de error absoluto, se define la:
Cota de Error relativo.-De un nmero aproximado como cualquier CRr no menor
que el error relativo de dicho nmero. Por consiguiente r CRr (1.6), en la
prctica se utiliza frecuentemente:
Er = rCRr (1.7).
Ejemplo:
Conocida la cota de error relativo, determinar los lmites del nmero exacto.
Solucin:
De (1.4) R = r Er, reemplazando en esta ecuacin 1.7 se obtiene: R = r rCRr, por
consiguiente R = r(1 CRr) (1.8) o tambin r(1 - CRr) R r(1 + CRr) (1.9).
Ejemplo 1.2:
Dado el nmero aproximado a = 42.53 y la cota de error relativo
CRa = 2 por mil, averiguar entre que lmites se encuentra el nmero exacto.
Solucin:
Aplicando (1.9) y reemplazando valores:
42.53(1 - 0.002) A 42.53(1 + 0.002); 42.53*0.998 A 42.53*1.002; entonces:
42.445 A 42.615
10
Digitos significativos.
Cualquier nmero positivo a puede representarse en forma decimal con una cantidad
finita o infinita de dgitos:
(1.10)
a = m10m + m-110m-1 + m-210m-2 + ... m-n+110m-n+1 + ...
donde los i son los dgitos del sistema numrico decimal,
m 0, es el dgito ms significativo, m es la potencia ms elevada de diez del
nmero a y n el nmero de dgitos.
Ejemplo 1.4:
14325.98 = 1*104 + 4*103 + 3*102 + 2*101 + 5*100 + 9*10-1 + 8*10-2
Se denomina dgito significativo de un nmero aproximado a cualquier dgito
diferente de cero o cualquier cero situado entre dgitos significativos; son adems
aquellos ceros que se utilizan para indicar posiciones a la derecha de un nmero o
de otro dgito significativo. Todos aquellos ceros que estn a la izquierda del nmero
y que solamente se indican para la posicin del punto decimal no se deben
considerar como dgitos significativos.
Ejemplos 1.5:
645.8004 tiene 7 dgitos significativos.
0.0000038 tiene 2 dgitos significativos.
78.4000 tiene 6 dgitos significativos.
78.4 tiene 3 dgitos significativos.
11
12
TERMINO
SENO
0.392699082
0.392699082
-0.010093189
0.382605893
0.000077825
0.382683718
-0.000000286
0.382683432
TERMINO
COSENO
-0.077106284
0.922893716
0.000990897
0.923884612
-0.000005094
0.923879519
13
c) Desarrollar una aplicacin Matlab y otra en Visual C#, para que parametrice el
ejemplo anterior, ingresando el ngulo en radianes x; y el nmero de trminos; la
aplicacin debe desplegar como resultados: sen(x); cos(x); tan(x), EAtan(x)
Solucin:
% ErroresTangente.m
% Calculo de errores de la tangente usando series de Taylor.
while 1
ang = input('Ingrese el ngulo en grados: ');
% Validacin:
nTer = input('Ingrese el nmero de trminos (salir < 1 o > 12: ');
if (nTer < 1 | nTer > 12) break; end;
% Conversin de grados a radianes:
x = pi * ang / 180;
% Clculo del seno mediante la serie de Taylor:
senx = 0;
aux = 1;
signo = 1;
for i = 1 : nTer
fac = 1;
for j = 2 : aux
fac = fac * j;
end;
senx = senx + signo * x ^ aux / fac;
signo = -signo;
aux = aux + 2;
end;
% Clculo del coseno mediante la serie de Taylor:
cosx = 0;
14
aux = 0;
signo = 1;
for i = 1 : nTer
fac = 1;
for j = 2 : aux
fac = fac * j;
end;
cosx = cosx + signo * x ^ aux / fac;
signo = -signo;
aux = aux + 2;
end;
% Calcular las tangentes aproximada y exacta:
tanx = senx / cosx;
TAN = tan(x);
% Calcular los errores:
erAbs = abs(TAN - tanx);
erRel = erAbs / abs(TAN);
% Desplegar resultados:
disp(sprintf('Angulo (radianes): %f', x));
disp(sprintf('Tangente (exacta): %f', TAN));
disp(sprintf('Seno: %f', senx));
disp(sprintf('Coseno: %f', cosx));
disp(sprintf('Tangente (aproximada): %f', tanx));
disp(sprintf('Error absoluto: %f', erAbs));
disp(sprintf('Error relativo: %f\n', erRel));
end;
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
Este libro no pretende, incorporar un curso de Matlab o Visual C#, sin embargo otras
aplicaciones de estas herramientas de desarrollo, se las describe ms adelante y en
el Apndice A para Matlab; Apndice B para Visual C#.
15
16
Como se puede apreciar, existen nueve cuadros de edicin, mismos que permitirn
ingresar los datos; y presentar los resultados; y dos botones, mismos que tienen la
funcionalidad:
- Procesar: Ejecuta el mtodo para calcular los errores y desplegar los resultados.
Y
- Salir: Finaliza la ejecucin del programa.
Se describe a continuacin el cdigo fuente:
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
// Declaracion de variables
double ang, x, sen = 0, cos = 0, tan, TAN, erAbs, erRel, fac;
int nTer, signo = 1, aux, i, j;
try
{
// Transferencia de captura de datos a variables:
ang = double.Parse(txtAngGrad.Text);
nTer = int.Parse(txtNTer.Text);
// Validacin:
if (nTer < 1 || nTer > 12) return;
// Conversin de grados a radianes:
x = Math.PI * ang / 180;
// Clculo del seno mediante la serie de Taylor:
aux = 1;
for (i = 1; i <= nTer; i ++) {
fac = 1;
for (j = 2; j <= aux; j ++) fac *= j;
sen += signo * Math.Pow(x, aux) / fac;
signo = -signo;
aux += 2;
}
// Clculo del coseno mediante la serie de Taylor:
aux = 0;
signo = 1;
for (i = 1; i <= nTer; i ++) {
fac = 1;
for (j = 2; j <= aux; j++) fac *= j;
cos += signo * Math.Pow(x, aux) / fac;
signo = -signo;
aux += 2;
}
// Calcular las tangentes aproximada y exacta:
tan = sen / cos;
TAN = Math.Tan(x);
// Calcular los errores:
erAbs = Math.Abs(TAN - tan);
erRel = erAbs / Math.Abs(TAN);
// Transferir variables a pantalla:
txtAngRad.Text = Math.Round(x, 9).ToString();
txtTanExacto.Text = Math.Round(TAN, 9).ToString();
17
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Salir:
private void cb_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
programa:
18
19
conociendo que
2,- Sea f(x) = 3x4 + 6x3 - 2x2 5x + 2. Para X = 5.03. Calcular el error absoluto y
relativo de f(x) que se comete al truncar a X a los decimales.
20
3.- La funcin arco seno, se encuentra definida por la siguiente serie de Taylor:
sen -1(x) = x + (x3/3)(1/2) + (x5/5)(1*3)/(2*4) + (x7/7)(1*3*5)/(2*4*6) +
Calcular el error absoluto y relativo cuando se toman los cuatro primeros trminos de
la serie.
21
22
2
CLCULOS DE RAICES REALES
DE ECUACIONES DE LA FORMA F (X) = 0
23
24
CAPITULO II
CALCULO DE RAICES REALES DE ECUACIONES DE LA FORMA F (X) = 0
Generalidades
Este captulo se centra en la localizacin de los ceros de funciones en el plano
cartesiano.
Ejemplos de ecuaciones de este tipo son:
1)
2)
3)
4)
5)
A*X + B = 0
A*X2 + B*X + C = 0
A0*Xn+ A1*Xn-1 + A2*Xn-2 + . . . + An-1*X + A n = 0
X2*log10(2*X) 8 = 0
X5 * cos 2*X + ex *sen X + 3 = 0
Dependiendo del caso, se puede aplicar un mtodo para resolver el problema, tales
mtodos pueden ser exactos o algebraicos, aproximados (grficos o numricos), etc.
As se tiene que para resolver la ecuacin del ejemplo 1, bastar con despejar X de
esta ecuacin, con lo que X = -B/A. En el ejemplo 2, se aplica la frmula para el
clculo de races de ecuaciones cuadrticas, con lo que
X=
B B 2 4 AC
2A
En el caso de los otros ejemplos, resulta difcil o casi imposible encontrar las races
reales, aplicando un mtodo exacto, por lo que se recurren a mtodos aproximados,
los cuales se analizan ampliamente en este captulo.
Mtodos Grficos:
Se revisa a continuacin, un mtodo grfico, el cual consiste en hacer Y = F(X),
seguidamente se elabora una tabla, dando valores en forma ascendente a la variable
independiente X, paralelamente se calcula el valor Y. Este grupo de pares
ordenados, se lo representa en un sistema de coordenadas rectangulares y
finalmente se unen tales puntos, generndose un diagrama de dispersin. La
interseccin entre la funcin F(X) y el eje horizontal, constituyen los ceros de la
funcin o las races reales de la ecuacin F(X) = 0.
Ejemplo 2.2:
25
Com
o puede observarse, la funcin F(X), corta el eje horizontal, en el rango de las
abscisas 3 y 3.5, lo que significa que la raz es aproximadamente 3.2.
Otro mtodo grfico consiste en transformar la ecuacin F(X) = 0 en f1(X) = f2(X) y
haciendo Y1 = f1(X), Y2 = f2(X), se elabora una tabla, asignando un grupo de valores
a X, calculndose Y1 y Y2. Los puntos de esta tabla, se los representa en un
diagrama de dispersin, en el cual se dibujan las funciones f1(X) y f2(X). Las
intersecciones de stas representan las races reales de la ecuacin F(X) = 0.
Ejemplo 2.3:
Dada la ecuacin: X2*log10(2*X) 8 = 0, encontrar una raz real, utilizando el mtodo
grfico alternativo.
26
Solucin:
Se hace 8/X = X*log10(2*X)
(2)
Entonces Y1 = f1(X) = 8/X, Y2 = f2(X) = X*log10(2*X)
Asignando a X los valores del ejemplo anterior, se tiene la tabla:
X
0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00 5.50
y1 = f1(X) 16.00 8.00 5.33 4.00 3.20 2.67 2.29 2.00 1.78 1.60 1.45
y2 = f2(X) 0.00 0.30 0.72 1.20 1.75 2.33 2.96 3.61 4.29 5.00 5.73
Los puntos que representan a las dos funciones, se los lleva al sistema de ejes de
coordenadas en el plano de una manera independiente, obtenindose as el
diagrama:
Mtodos Numricos:
Estos involucran la implementacin de un algoritmo y se separan en dos fases:
Localizacin del cambio de signo de la funcin F(X), y
27
Afinamiento de los valores de las abscisas (XI y XD), entre las que se encuentra
la raz real.
XI = LIX
XD = XI
YD = F (XD)
XI = XD
YI = YD
XD = XD +DX
YD = F (XD)
P = YI * YD
REPETIR DESDE EL PASO (d) MIENTRAS (P > 0) Y (XD < LSX)
SI P = 0 ENTONCES
SI YI = 0 ENTONCES
XI ES RAIZ EXACTA
CASO CONTRARIO
XD ES RAIZ EXACTA
XD = XD + DX
FIN SI
FIN SI
Ejemplo 2.4:
Localizar los cambios de signo de la ecuacin: F(X) = X2*log10(2*X) 8 = 0
28
Solucin:
Se adoptan los valores LIX = 0.5, LSX = 5, DX = 0.5, realizando el seguimiento del
algoritmo, se genera la tabla siguiente:
XI
0.50
0.50
1.00
1.50
2.00
2.50
3.00
3.50
4.00
4.50
5.00
XD
0.50
1.00
1.50
2.00
2.50
3.00
3.50
4.00
4.50
5.00
5.50
YI = F(XI)
-8.00
-7.70
-6.93
-5.59
-3.63
-1.00
2.35
6.45
11.32
17.00
YD = F(XD) P = YI*YD
-8.00
-7.70
(+)
-6.93
(+)
-5.59
(+)
-3.63
(+)
-1.00
(+)
2.35
(-)
6.45
(+)
11.32
(+)
17.00
(+)
23.50
(+)
XD
-4.00
-3.00
-2.00
-1.00
0.00
1.00
2.00
3.00
4.00
5.00
YI = F(XI)
-32.00
-2.00
10.00
10.00
4.00
-2.00
-2.00
10.00
40.00
YD = F(XD) P = YI*YD
-32.00
-2.00
(+)
10.00
(-)
10.00
(+)
4.00
(+)
-2.00
(-)
-2.00
(+)
10.00
(-)
40.00
(+)
94.00
(+)
Por consiguiente, existen races en los intervalos [-3, -2], [0, 1] y [2, 3].
Segunda fase: Afinamiento de los valores que se encuentran cercanos a la raz.
29
30
Teorema 2.1:
Dada
una
funcin continua
F(X)
en
el
intervalo
cerrado [A, B]
en las abscisas.
Existe por lo
menos una raz
real,
si
F(A)*F(B) <= 0.
Si
adems
F(X), mantiene
el signo en ese
intervalo, entonces existe una sola raz (ver figuras 2.3 y 2.5).
Existen varios mtodos para calcular races reales, los que analizaremos en este
captulo son:
Biseccin.
Aproximaciones sucesivas.
Aproximaciones sucesivas modificado.
Rgula - falsi.
Secante y.
Newton - Raphson.
Muller.
Cuyos algoritmos, necesitan las siguientes entradas: F(X), XI, XD, PRECISION
Donde F(X) es la funcin continua; XI y XD son las abscisas izquierda y derecha,
respectivamente, entre las que puede estar por lo menos una raz real; PRECISION
es el margen de error que se produce al efectuar el clculo de una raz.
Seguidamente se revisan cada uno de ellos.
31
Mtodo de la
Biseccin:
Consiste
en
acortar
el
intervalo de las
abscisas,
calculando
su
valor medio Xm,
luego se calcula
F(Xm)
y
se
aproxima a la
abscisa
cuya
ordenada tiene el
mismo signo de
F(Xm). (Ver la figura 2.6).
Algoritmo para el mtodo de la Biseccin
a)
b)
c)
d)
e)
YI = F(XI)
Xm = (XI+XD) / 2
Ym= F(Xm)
P = YI * Ym
SI P > 0 ENTONCES
XI = Xm
YI = Ym
CASO CONTRARIO
XD = Xm
FIN SI
f) REPETIR DESDE EL PASO (b) MIENTRAS |Ym| > PRECISION
g) ESCRIBIR RESULTADOS
Ejemplo 2.6:
Calcular la raz de la ecuacin del ejemplo 2.5, para XI = -3, XD = -2, PRECISION =
1E-6,
usando el algoritmo del mtodo de la biseccin.
32
Solucin:
Para calcular la raz, se elabora la tabla:
ITERACION
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
XI
-3.000000
-3.000000
-3.000000
-3.000000
-2.937500
-2.906250
-2.906250
-2.898438
-2.898438
-2.896484
-2.895508
-2.895508
-2.895264
-2.895142
-2.895142
-2.895111
-2.895111
-2.895111
-2.895107
-2.895107
-2.895107
-2.895107
-2.895107
-2.895107
-2.895107
XD
Xm=(XI+XD)/2 YI = F(XI) Ym= F(Xm) P = YI*Ym
-2.000000
-2.500000
-2.000000 5.875000
(-)
-2.500000
-2.750000
-2.000000 2.453125
(-)
-2.750000
-2.875000
-2.000000 0.361328
(-)
-2.875000
-2.937500
-2.000000 -0.784912
(+)
-2.875000
-2.906250
-0.784912 -0.203278
(+)
-2.875000
-2.890625
-0.203278 0.081142
(-)
-2.890625
-2.898438
-0.203278 -0.060537
(+)
-2.890625
-2.894531
-0.060537 0.010435
(-)
-2.894531
-2.896484
-0.060537 -0.025018
(+)
-2.894531
-2.895508
-0.025018 -0.007283
(+)
-2.894531
-2.895020
-0.007283 0.001578
(-)
-2.895020
-2.895264
-0.007283 -0.002852
(+)
-2.895020
-2.895142
-0.002852 -0.000637
(+)
-2.895020
-2.895081
-0.000637 0.000471
(-)
-2.895081
-2.895111
-0.000637 -0.000083
(+)
-2.895081
-2.895096
-0.000083 0.000194
(-)
-2.895096
-2.895103
-0.000083 0.000056
(-)
-2.895103
-2.895107
-0.000083 -0.000014
(+)
-2.895103
-2.895105
-0.000014 0.000021
(-)
-2.895105
-2.895106
-0.000014 0.000004
(-)
-2.895106
-2.895107
-0.000014 -0.000005
(+)
-2.895106
-2.895107
-0.000005 -0.000001
(+)
-2.895106
-2.895106
-0.000001 0.000001
(-)
-2.895106
-2.895106
-0.000001 0.000000
(-)
-2.895106
-2.895107
-0.000001 0.000000
(+)
33
Conforme a lo que se observa en las figuras 2.7 (a) y (b) existe convergencia para
cuando |fm'(X) | 1; en (c) y (d) no existe convergencia porque |fm'(X) | > 1, por
consiguiente hay que considerar esta situacin en el algoritmo.
Algoritmo para el mtodo de aproximaciones sucesivas:
a) Ym = (XI+XD)/2
b) Xm = Ym
c) Ym = fm(Xm)
d) REPETIR DESDE EL PASO (b) MIENTRAS ( |Ym - Xm| > PRECISION ) Y
( |fm'(Ym)| 1)
e) ESCRIBIR RESULTADOS
Ejemplo 2.7:
Calcular la raz de la ecuacin del ejemplo 2.5, para XI = -3, XD = -2, PRECISION =
1E-6,
usando el algoritmo del mtodo de aproximaciones sucesivas.
34
Solucin:
Debe elaborarse una ecuacin Ym = Xm = fm(Xm), tal que |fm'(-2.78)| 1.
Despejando X, del primer trmino de la ecuacin F(X) = 0, se tiene que:
Ym = Xm = fm(Xm) = (7*X - 4)(1/3), por consiguiente, fm'(Xm) = (7/3)*(7*X - 4)(-2/3).
Para calcular la raz, se elabora la tabla:
ITERACION
0
1
2
3
4
5
6
7
8
9
10
11
Xm
-2.500000
-2.780649
-2.862886
-2.886109
-2.892600
-2.894408
-2.894912
-2.895052
-2.895091
-2.895102
-2.895105
Ym =
fm(Xm)
-2.500000
-2.780649
-2.862886
-2.886109
-2.892600
-2.894408
-2.894912
-2.895052
-2.895091
-2.895102
-2.895105
-2.895106
YP =
fm'(Ym)
|Ym - Xm|
0.284688
0.280125
0.278869
0.278520
0.278424
0.278397
0.278389
0.278387
0.278386
0.278386
0.278386
0.280649
0.082237
0.023223
0.006491
0.001809
0.000504
0.000140
0.000039
0.000011
0.000003
0.000001
35
YI = fm(XI)
YD = fm(XD)
TGT = (YD - YI) / (XD - XI)
BETA = 1 / (1 - TGT)
XI = XD
YI = YD
XD = XD + BETA* (YD - XD)
REPETIR DESDE EL PASO (b) MIENTRAS |XD - XI| > PRECISION
ESCRIBIR RESULTADOS
Ejemplo 2.8:
Calcular la raz de la ecuacin del ejemplo 2.5, para XI = -3, XD = -2, PRECISION =
1E-6,
usando el algoritmo del mtodo de aproximaciones sucesivas modificado.
Solucin:
Despejando X, del segundo trmino de la ecuacin F(X) = 0, se tiene que:
Y = X = fm(X) = (X3 + 4) / 7. Para calcular la raz, se elabora la tabla:
ITERACION
0
1
2
3
4
5
6
XI
-3.000000
-2.000000
-2.833333
-2.935065
-2.893916
-2.895084
-2.895107
XD
-2.000000
-2.833333
-2.935065
-2.893916
-2.895084
-2.895107
-2.895107
YI = fm(XI) YD = fm(XD)
-3.285714
-0.571429
-0.571429
-2.677910
-2.677910
-3.040633
-3.040633
-2.890831
-2.890831
-2.895026
-2.895026
-2.895107
-2.895107
-2.895107
TGT
BETA
|XD - XI|
2.714286
2.527778
3.565486
3.640455
3.590627
3.592104
-0.583333
-0.654545
-0.389790
-0.378723
-0.386007
-0.385787
0.833333
0.101732
0.041149
0.001168
0.000022
0.000000
36
YI = F(XI)
YD = F(XD)
Xm = XI - YI * (XD - XI) / (YD - YI)
Ym= F(Xm)
P = YI * Ym
SI P > 0 ENTONCES
XI = Xm
YI = Ym
CASO CONTRARIO
XD = Xm
YD = Ym
FIN SI
g) REPETIR DESDE EL PASO (c) MIENTRAS |Ym| > PRECISION
h) ESCRIBIR RESULTADOS.
Ejemplo 2.9:
37
Calcular la raz de la ecuacin del ejemplo 2.5, para XI = -3, XD = -2, PRECISION =
1E-6,
usando el algoritmo del mtodo de Rgula - Falsi.
Solucin:
Para calcular la raz, se elabora la tabla:
ITERACION
0
1
2
3
4
5
6
XI
-3.000000
-3.000000
-3.000000
-3.000000
-3.000000
-3.000000
-3.000000
XD
-2.000000
-2.833333
-2.892054
-2.894959
-2.895099
-2.895106
-2.895106
YI = F(XI) YD = F(XD)
-2.000000 10.000000
-2.000000 1.087963
-2.000000 0.055307
-2.000000 0.002681
-2.000000 0.000130
-2.000000 0.000006
-2.000000 0.000000
Xm
-2.833333
-2.892054
-2.894959
-2.895099
-2.895106
-2.895106
Ym = F(Xm) P= YI*Ym
1.087963
0.055307
0.002681
0.000130
0.000006
0.000000
(-)
(-)
(-)
(-)
(-)
(-)
Mtodo de la Secante:
Este mtodo utiliza
la
frmula
de
recurrencia del de
Regula
Falsi,
calculndose valores
consecutivos en las
abscisas de uno de
los
extremos
y
manteniendo fijo el
otro extremo. El
mtodo
tiene
la
restriccin de que la
funcin F(X) en el
intervalo donde se
encuentra la raz debe mantener F(X) el mismo signo, como lo indica la figura 2.10.
YI = F(XI)
YD = F(XD)
XI = XI - YI * (XD - XI) / (YD - YI)
YI = F(XI)
REPETIR DESDE EL PASO (c) MIENTRAS |F(XI)| > PRECISION
ESCRIBIR RESULTADOS.
38
Ejemplo 2.10:
Calcular la raz de la ecuacin del ejemplo 2.5, para XI = -3, XD = -2, PRECISION =
1E-6,
usando el algoritmo del mtodo de la Secante.
Solucin:
Se har variar el punto que se encuentra a la izquierda y se mantendr fijo el punto
de la derecha. Los resultados del proceso se muestran en la siguiente tabla:
ITERACION
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
XI
-3.000000
-2.833333
-2.935065
-2.870721
-2.910541
-2.885557
-2.9011
-2.891378
-2.897438
-2.893653
-2.896015
-2.89454
-2.89546
-2.894886
-2.895244
-2.895021
-2.89516
-2.895073
-2.895127
-2.895093
-2.895115
-2.895101
-2.89511
-2.895105
-2.895108
-2.895106
-2.895107
-2.895106
-2.895107
-2.895106
-2.895107
YI = F(XI)
-2.000000
1.087963
-0.738974
0.437324
-0.282132
0.172490
-0.109059
0.067524
-0.042358
0.026356
-0.016483
0.010275
-0.006418
0.004004
-0.002500
0.001560
-0.000974
0.000608
-0.000379
0.000237
-0.000148
0.000092
-0.000058
0.000036
-0.000022
0.000014
-0.000009
0.000005
-0.000003
0.000002
-0.000001
XD
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
-2.000000
YD
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
10.000000
39
Xm = (XI + XD) / 2
Ym = F(Xm)
YP = F'(Xm)
Xm = Xm - Ym / YP
REPETIR DESDE EL PASO (b) MIENTRAS |Ym| > PRECISION
ESCRIBIR RESULTADOS.
Ejemplo 2.11:
Calcular la raz de la ecuacin del ejemplo 2.5, para XI = -3, XD = -2, PRECISION =
1E-6,
usando el algoritmo del mtodo de Newton - Raphson.
Solucin:
F(X) = X3- 7*X + 4; por consiguiente F'(X) = 3 * X2 - 7; la tabla siguiente muestra los
resultados del proceso de este algoritmo:
ITERACION
0
1
2
3
4
5
Xm
Ym = F(Xm) YP = F'(Xm)
-2.500000
-3
5.875000
11.750000
-2.9
-2.000000
20.000000
-2.895118 -0.089000
18.230000
-2.895107 -0.000207
18.145124
-2.895107 0.000000
18.144925
40
Como se puede observar hasta aqu, este es uno de los mtodos ms simples
porque involucra pocas operaciones, y a la vez el ms eficiente.
Ejemplo 2.12:
Calcular la raz de la ecuacin del ejemplo 2.4, para XI = 3, XD = 3.5, PRECISION =
1E-6,
usando el algoritmo del mtodo de Newton - Raphson.
Solucin:
F(X) = X2*log10(2*X) 8; por consiguiente F'(X) = 2*X*log10(2*X) + 2*X / ln(10); la tabla
siguiente muestra los resultados del proceso de este algoritmo:
ITERACION
0
1
2
3
4
5
6
7
8
9
Xm
Ym = F(Xm) YP = F'(Xm)
3.250000
3.177666
0.586397
8.106851
3.163712
0.109740
7.864297
3.161213
0.019541
7.817668
3.160772
0.003443
7.809322
3.160694
0.000606
7.807850
3.160681
0.000106
7.807591
3.160678
0.000019
7.807545
3.160678
0.000003
7.807537
3.160678
0.000001
7.807536
Mtodo de Muller:
Se basa en el mtodo de la Secante, reemplazando la lnea recta por una parbola
de la forma:
Y = Parabola(x) =
a(x x2)2 + b(x x2) + c .
(Observar la figura 2.12)
En la figura 2.12: En la RAIZ:
Parabola(X) = F(X) (2.17);
Para encontrar los coeficientes a,
b, c, se plantea el siguiente
sistema de ecuaciones:
F(x0) = a(x0 x2)2+b(x0 x2)+c
(2.18)
F(x1) = a(x1 x2)2+b(x1 x2)+c
(2.19)
F(x2) = c (2.20)
41
Resolviendo el sistema:
a = (d1 d0)/(h1 + h0) (2.21);
b = ah1 + d1 (2.22), y
c = F(x2) (2.23).
Donde:
h0 = x1 x0 (8); h1 = x2 x1 (2.24);
d0 =(F(x1) F(x0)) / h0 (2.25);
d1 =(F(x2) F(x1)) / h1 (2.26);
Para encontrar la frmula de recurrencia se puede aplicar la expresin para el clculo
de races de la ecuacin cuadrtica; sin embargo, por problemas de redondeo se
aplica la relacin:
2c
; (2.27 )
x3 = x2
b b 2 4ac
Como se puede apreciar en la frmula de recurrencia, se tienen dos valores en el
denominador de la fraccin, en este caso, para lograr la convergencia, se debe tomar
el mximo en valor absoluto de los dos. Seguidamente se describe el
correspondiente algoritmo.
Algoritmo para el mtodo de Muller:
a)
b)
c)
d)
e)
f)
g)
h)
i)
j)
k)
l)
m)
n)
o)
p)
Ejemplo 2.13:
Calcular la raz de la ecuacin del ejemplo 2.5, para XI = -3, XD = -2, PRECISION =
1E-6,
usando el algoritmo del mtodo de Muller.
42
Solucin:
F(X) = X3- 7*X + 4; las tablas siguientes muestran los resultados del proceso de este
algoritmo:
Iter
x0
x1
-3
-2
2
3
x2
-2.5
-2
-2.5 2.89303534
-2.5 2.89303534 2.89514769
Iter
h0
x3
2.89303534
2.89514769
2.89510651
y0 =
f(x0)
y1 =
f(x1)
-2
10
y2 = f(x2)
y3 = f(x3)
5.875 0.03754400
10
5.875 0.03754400 0.00074717
5.875 0.037544 0.00074717 0.00000003
h1
d0
d1
a
b
c
R
DEN
0.500000 12.000000 8.250000 7.500000 12.000000 5.875000 17.895530 29.895530
0.393035 8.250000 14.852242 7.393035 17.757966 0.037544 17.789199 35.547165
0.002112 14.852242 18.127298 8.288183 18.144806 0.000747 18.144123 36.288929
1 1.000000
2 0.500000
3 0.393035
x0
x1
x2
x3
3.5
3.5
3.25
3.25
3.1607115
5
3.1607115
5
3.1606775
8
Iter
h0
h1
d0
1 0.500000 0.250000 6.698179
2 0.250000 0.089288 7.064215
y0 = f(x0)
0.996638
7
2.352451
0
d1
y1 = f(x1)
y2 = f(x2)
2.352451
0
0.586397
3
0.586397
3 0.0002185
0.000218
5 0.0000000
y3 = f(x3)
DEN
7.064215 1.464141
6.698179 0.586397
6.436718 13.134897
6.565001 1.471353
6.433627 0.000219
6.433527 12.867153
43
44
% FPOL.m
function y = FPOL(x)
y = x ^ 3 - 7 * x + 4;
return
%derivada.m
%Por definicion de la derivada
function y = derivada(x)
h = 1e-12;
y = (FPOL(x + h) - FPOL(x)) / h;
return
% fm.m
function y = fm(x)
signo = 1;
r = 7 * x - 4;
if r < 0
r = -r;
signo = -1;
end;
y = signo * r ^ (1 / 3);
return
%derivadafm.m
%Por definicion de la derivada
function y = derivadafm(x)
h = 1e-12;
y = (fm(x + h) - fm(x)) / h;
return
% Localiza_Raiz.m
function [p, xi, xd] = Localiza_Raiz(nombre_f, li, dx, b)
xd = li;
yd = feval(nombre_f, xd);
while 1
xi = xd;
yi = yd;
xd = xd + dx;
yd = feval(nombre_f, xd);
p = yi * yd;
if (p <= 0) | (xd >= b)
45
break;
end;
end;
if (p == 0)
if (xi == 0)
fprintf('Raiz exacta = %12.5f\n', xi);
elseif (xd == 0)
fprintf('Raiz exacta = %12.5f\n', xd);
xd = xd + dx;
end;
end;
return
% Biseccion.m
function b = Biseccion(nombre_f, xi, xd)
NIter = 100; Precision = 1e-7;
fprintf('Clculo de la raiz entre %12.5f y %12.5f, por el mtodo de la biseccin:\n', xi,
xd);
yi = feval(nombre_f, xi);
for Iter = 1 : NIter
x = (xi + xd) / 2;
y = feval(nombre_f, x);
if (abs(y) <= Precision)
break;
elseif (y * yi > 0)
xi = x; yi = y;
else
xd = x;
end;
end;
y = abs(y);
if (y > Precision)
fprintf('No se alcanza la precisin requerida: %12.5f ', y);
end;
fprintf('Raz =: %12.5f, en %d iteraciones\n', x, Iter);
b = x;
return
% AproxSuc.m
function b = AproxSuc(nombre_f, n_fp, xi, xd)
NIter = 100; Precision = 1e-7;
fprintf('Clculo de la raiz entre %12.5f y %12.5f, por el mtodo de aproximaciones
sucesivas:\n', xi, xd);
x = (xi + xd) / 2;
for Iter = 1 : NIter
yp = feval(n_fp, x);
46
if abs(yp) >= 1
fprintf('No existe convergencia\n');
break;
end;
x0 = x;
x = feval(nombre_f, x);
if (abs(x - x0) <= Precision)
break;
end;
end;
if abs(yp) < 1
x0 = abs(x - x0);
if (x0 > Precision)
fprintf('No se alcanza la precisin requerida: %12.5f ', x0);
end;
fprintf('Raz =: %12.5f, en %d iteraciones\n', x, Iter);
b = x;
else
b = -11111;
end;
return
% AproxMod.m
function b = AproxMod(nombre_f, xi, xd)
NIter = 100; Precision = 1e-7;
fprintf('Raiz entre %12.5f y %12.5f, por el mtodo de aproximaciones sucesivas
modificado:\n', xi, xd);
dif = xd - xi;
yi = feval(nombre_f, xi);
for Iter = 1 : NIter
yd = feval(nombre_f, xd);
tgt = (yd - yi) / dif;
alfa = 1 / (1 - tgt);
xi = xd;
yi = yd;
xd = xd + alfa * (yd - xd);
dif = xd - xi;
if (abs(dif) <= Precision)
break;
end;
end;
dif = abs(dif);
if (dif > Precision)
fprintf('No se alcanza la precisin requerida: %12.5f ', dif);
end;
fprintf('Raz =: %12.5f, en %d iteraciones\n', xd, Iter);
47
b = xd;
return
% RegulaFalsi.m
function x = RegulaFalsi(nombre_f, xi, xd)
NIter = 100; Precision = 1e-7;
fprintf('Raiz entre %12.5f y %12.5f, por el mtodo de Rgula - Falsi:\n', xi, xd);
yi = feval(nombre_f, xi);
yd = feval(nombre_f, xd);
for Iter = 1 : NIter
x = xi - yi * (xd - xi) / (yd - yi);
y = feval(nombre_f, x);
if y * yi > 0
xi = x;
yi = y;
else
xd = x;
yd = y;
end;
y = abs(y);
if (y <= Precision)
break;
end;
end;
if (y > Precision)
fprintf('No se alcanza la precisin requerida: %12.5f ', y);
end;
fprintf('Raz =: %12.5f, en %d iteraciones\n', x, Iter);
return
% Secante.m
function xi = Secante(nombre_f, xi, xd)
NIter = 100; Precision = 1e-7;
fprintf('Raiz entre %12.5f y %12.5f, por el mtodo de la Secante:\n', xi, xd);
yi = feval(nombre_f, xi);
yd = feval(nombre_f, xd);
for Iter = 1 : NIter
xi = xi - yi * (xd - xi) / (yd - yi);
yi = feval(nombre_f, xi);
if (abs(yi) <= Precision)
break;
end;
end;
yi = abs(yi);
if (yi > Precision)
fprintf('No se alcanza la precisin requerida: %12.5f ', yi);
48
end;
fprintf('Raz =: %12.5f, en %d iteraciones\n', xi, Iter);
return
% Newton_Raphson.m
function nr = Newton_Raphson(nombre_f, nombre_fp, xi, xd)
NIter = 100; Precision = 1e-7;
fprintf('Clculo de la raiz entre %12.5f y %12.5f, por el mtodo de Newton Raphson:\n', xi, xd);
x = (xi + xd) / 2;
for Iter = 1 : NIter
x = x - feval(nombre_f, x) / feval(nombre_fp, x);
y = abs(feval(nombre_f, x));
if (y <= Precision)
break;
end;
end;
if (y > Precision)
fprintf('No se alcanza la precisin requerida: %12.5f ', y);
end;
fprintf('Raz =: %12.5f, en %d iteraciones\n', x, Iter);
nr = x;
return
% Muller.m
function x3 = Muller(nombre_f, x0, x1)
NIter = 100; Precision = 1e-7;
fprintf('Raiz entre %12.5f y %12.5f, por el mtodo de Muller:\n', x0, x1);
x2 = (x0 + x1) / 2;
y0 = feval(nombre_f, x0);
y1 = feval(nombre_f, x1);
y2 = feval(nombre_f, x2);
for Iter = 1 : NIter
h0 = x1 - x0;
h1 = x2 - x1;
d0 = (y1 - y0) / h0;
d1 = (y2 - y1) / h1;
a = (d1 - d0) / (h1 + h0);
b = a * h1 + d1;
c = y2;
rc = sqrt(b^2 - 4 * a * c);
den = b - rc;
if abs(b + rc) > abs(b - rc)
den = b + rc;
end;
x3 = x2 - 2 * c / den;
49
y3 = feval(nombre_f, x3);
if (abs(y3) <= Precision)
break;
end;
x0 = x1;
x1 = x2;
x2 = x3;
y0 = y1;
y1 = y2;
y2 = y3;
end;
y3 = abs(y3);
if (y3 > Precision)
fprintf('No se alcanza la precisin requerida: %12.5f ', y3);
end;
fprintf('Raz =: %12.5f, en %d iteraciones\n', x3, Iter);
return
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
50
51
52
Botn Procesar:
private void cb_procesar_Click(object sender, EventArgs e)
{
double liminf, limsup, delta, estado, xi, xd = 0;
try
{
liminf = double.Parse(txtInicial.Text);
limsup = double.Parse(txtFinal.Text);
delta = double.Parse(txtIncremento.Text);
CRaicesReales.a = liminf;
CRaicesReales.b = limsup;
CRaicesReales.dx = delta;
CRaicesReales.maxIter = 100;
CRaicesReales.precision = 1e-8;
CRaicesReales.lb = listBox1;
CRaicesReales.punf = F;
CRaicesReales.punfp = FP;
CRaicesReales.punfm = fm;
CRaicesReales.punfpm = fmp;
xi = liminf;
xd = xi;
do {
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
estado = CRaicesReales.LocCambioSigno();
xi = CRaicesReales.xi;
xd = CRaicesReales.xd;
if (estado < 0)
{
CRaicesReales.Biseccion();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.AproxSucesivas();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.AproxSucesivasMod();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
53
CRaicesReales.NewtonRaphson();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.RegulaFalsi();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.Secante();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.Muller();
}
xi = xd;
}
while ((xd += delta) < limsup);
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void cb_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
54
55
56
3
1
0
0
POSITIVAS
NEGATIIVAS
1
1
2
0
57
X-K
R = F(K)
Q(X)
A1
+ B0 * K
B1
A2
+ B1 * K
B2
...
...
...
AN-1
+ BN-2 * K
BN-1
AN
+ BN-1 * K
R = F(K)
58
F(K) = Q(K) (2.58), esto significa que para calcular F'(K) debe aplicarse otra divisin
sinttica sobre el esquema de HORNER, con lo que se tiene:
A0
B0
C0
A1
+ B0 * K
B1
+ C0 * K
C1
A2
+ B1 * K
B2
+ C1 * K
C2
...
...
...
...
...
AN-1
+ BN-2 * K
BN-1
+ CN-2 * K
F'(K)
AN
+ BN-1 * K
R = F(K)
-5
-5
16
22
48
8
6
11
28
24
78
F(2) = 43
14
39
F'(2) = 102
59
Ejemplo 2.18:
Calcular la raz entre XI = 2 y XD = 3, de la ecuacin del ejemplo 2.5, con una
precisin de 0.0001, empleando los esquemas de HORNER y el mtodo de Newton Raphson.
Solucin:
A este mtodo combinado se lo conoce tambin como de Birge - Vieta. los resultados
del proceso de estos algoritmos, se presentan a continuacin:
Xm = (2 +3) / 2 = 2.5
Iteracin 1:
1
-7
2.5
2.5
6.25
F(X) =-1.875
2.5
2.5
-0.75
12.5
2.125
F'(X) =11.75
-7
2.3191
5.3785
-3.7606
2.3191
2.3191
-1.6215
10.7569
F(X) = 0.2394
4.6383
F'(X) = 9.1354
-7
2.2929
5.2576
-3.9952
2.2929
2.2929
-1.7424
10.5152
F(X) = 0.0048
4.5859
F'(X) = 8.7728
-7
2.2924
5.2551
-4.0000
2.2924
2.2924
-1.7449
10.5102
F(X) = 0.0000
4.5848
F'(X) = 8.7653
2.3191
2.2929
2.2924
60
El esquema de HORNER, se adapta a los algoritmos que calculan races reales, los
cuales emplean la ecuacin F(X) = 0, es decir, Biseccin, Regula - Falsi, Secante,
Newton Raphson y Muller. Se incorporan entonces, estos cinco mtodos e los
programas respectivos. Para cada lenguaje, se elabora un solo programa, desde el
cual se ingresan el grado del polinomio, los coeficientes del polinomio,
almacenndolos en un arreglo, el intervalo en el que se encuentran todas las races
reales de la ecuacin; y el incremento de variacin en las abscisas. El programa
llama a cada funcin o subrutina, de una manera adecuada para localizar el cambio
de signo y luego prueba a cada uno de los mtodos para calcular la raz y la
despliega por pantalla. Se describen a continuacin una lista de funciones y
subrutinas.
1) Una funcin que reciba el valor de X y calcule y devuelva F(X), usando el mtodo
de la divisin sinttica o esquema de HORNER.
2) Una funcin que reciba el valor de X y calcule y devuelva F(X) , usando el
mtodo de la divisin sinttica o esquema de HORNER..
3) Una funcin que reciba lmite inicial (XI), el lmite final (LSX), y el incremento (DX)
y calcule y devuelva el intervalo del cambio de signo XI y XD, y adems el estado
en el nombre de la funcin.
4) Una subrutina que reciba el intervalo del cambio de signo XI y XD y calcule e
imprima la raz en ese intervalo.
5) Una funcin para ingresar desde el teclado el grado y los coeficientes del
polinomio y almacenarlos a estos ltimos en un arreglo.
6) Una subrutina para desplegar sobre la pantalla los coeficientes del polinomio.
Se describe a continuacin el cdigo fuente necesario y los resultados en Matlab:
% raices2.m
% Calculo de raices reales de funciones polinomicas.
% FUNCION PRINCIPAL:
global A N
nf = 'HORNER';
nfp = 'HORNERP';
[N, A] = LeePolinomio;
li = input('Ingrese el lmite izquierdo: ');
ld = input('Ingrese el lmite derecho: ');
dx = input('Ingrese el incremento de variacin: ');
xi = li;
while 1
61
62
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
63
64
if (CHorner.n > 0)
{
// Definir el arreglo de coeficientes de la funcion:
CHorner.a = new double[CHorner.n + 1];
// Leer el arreglo de coeficientes y generar los resultados en el listbox:
CLectEscrMatrices.n = CHorner.n + 1;
CLectEscrMatrices.b = CHorner.a;
CLectEscrMatrices.obent = objes;
CLectEscrMatrices.lb = listBox1;
CLectEscrMatrices.LeeVector();
listBox1.Items.Add("CALCULO DE RAICES DE LA FUNCION POLINOMICA DE GRADO " +
CHorner.n.ToString());
CLectEscrMatrices.EscribeVector(" DE COEFICIENTES");
CRaicesReales.a = liminf;
CRaicesReales.b = limsup;
CRaicesReales.dx = delta;
CRaicesReales.maxIter = 100;
CRaicesReales.precision = 1e-8;
CRaicesReales.lb = listBox1;
CRaicesReales.punf = CHorner.horner;
CRaicesReales.punfp = CHorner.fphorner;
xi = liminf;
xd = xi;
do
{
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
estado = CRaicesReales.LocCambioSigno();
xi = CRaicesReales.xi;
xd = CRaicesReales.xd;
if (estado < 0)
{
CRaicesReales.Biseccion();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.NewtonRaphson();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.RegulaFalsi();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.Secante();
CRaicesReales.xi = xi;
CRaicesReales.xd = xd;
CRaicesReales.Muller();
}
xi = xd;
}
while ((xd += delta) < limsup);
// Cerrar flujo de entrada:
objes.CerrarLectura();
}
}
catch (Exception ex)
{
65
}
}
MessageBox.Show(ex.Message);
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
66
67
68
Teorema 2.4
Si una funcin polinmica tiene grado impar, entonces existe por lo menos una raz
real.
Teorema 2.5
El residuo de la divisin entre una funcin polinmica F(X) y un polinomio de la forma
X2 + p*X + q es BN-1 * X + BN-1 * p + BN, donde BN-1 y BN son los coeficientes de X-1 y
X-2, respectivamente del polinomio divisor (Q(X)) de esa divisin.
Ejemplo 2.19:
Dado el polinomio: F(X) = 3X6 + 2X5 - 4X4 + 2X3 + 3X2 - X + 2 y el divisor
D(X) = X2 + p*X + q = X2 + 2X + 1, determinar el cuociente Q(X) y BN-1 y BN.
Solucin:
Efectuando la divisin de F(X) para de D(X):
3X6
-3X6
+ 2X5
- 6X5
- 4X4
- 3X4
- 4X5
4X5
- 7X4
+ 8X4
X4
- X4
+ 2X3
+ 4X3
+ 6X3
- 2X3
4X3
- 4X3
+ 3X2
- X2
+ 2X2
- 8X2
- 6X2
6X2
-X
- 4X
- 5X
+12X
7X
69
+ 2 X2 + 2X + 1
3X4 - 4X3 + X2 + 4X 6
+6
+8
B0
A1
A2
- B0 * q
- B0 * p - B1 * p
B1
B2
...
...
...
...
AN-2
AN-1
AN
- BN-4 * q - BN-3 * q - BN-2 * q
- BN-3 * p - BN-2 * p - BN-1 * p
BN-2
BN-1
BN
p
q
Ejemplo 2.20:
Dado el polinomio del ejemplo 2.19, mediante el esquema de Newton - Bairstow,
determinar BN-1 y BN.
Solucin:
El esquema que se genera es:
3
2
-4
2
-3
4
-6
8
-2
3
-4
1
4
3
-1
-8
-6
-1
-4
12
BN-1 = 7
2
6
-14
BN = -6
p=2
q=1
Las siguientes frmulas permiten determinar las cotas, para races complejas:
||RAIZ| (|A0 | + | Amx|) / |A0 |
(2.87);
(2.88)
70
aplicar la frmula : X = ( P P 2 4 * Q ) / 2
71
Ejemplo 2.3:
Calcular las races de la ecuacin:
F(X) = 2X5 - 4X4 + 7X3 - 14X2 + 6X - 10 = 0; usar una PRECISION DE 1E-7
Solucin:
Clculo de races reales:
Localizacin del cambio de signo:
2
-4
7
-14
6
-10 1
2
-2
5
-9
-3
2
-2
5
-9
-3
F(1) = -13
2
2
-4
4
0
7
0
7
-14
14
0
6
0
6
-10 2
12
F(2) = 2
Como F(1) = -13 y F(2) = 2, una raz est real entre 1 y 2, procedemos al calculo,
para un
X0 = (1 + 2)/2 = 1.5:
2
-4
7
-14
6
-10
X0=1.5
3
-1.5
8.25
-8.625
-3.938
2
-1
5.5
-5.75
-2.625
F(X0) =-13.94
3
3
12.75
10.5
2
2
8.5
7 F'(X0) =7.875
X1 = X0 - F(X0)/F'(X0) = 1.5 - (-13.938)/ 7.875 = 3.26984127
Como F(X0) > PRECISION, se procede a la iteracin siguiente:
Iteracin 2:
2
-4
6.53968
2 2.53968
6.53968
2 9.07937
7
8.30436
15.3044
29.6881
44.9924
6
-10
-14
404.984
50.0428 117.854
F(X)=394.984
36.0428 123.854
147.118 598.907
183.161 F'(X)=722.76
3.26984
7
3.93986
10.9399
18.7731
29.713
-10
-14
6
133.471
29.7931
43.01
F(X)=123.471
15.7931
49.01
80.9188
263.38
96.7118 F'(X)=312.3
9
72
2.72335
7
1.52771
8.52771
12.3678
20.8955
6
-14
-10
19.8534 13.6272
45.6942
5.85337 19.6272
F(X)=35.6942
48.6469 126.882
54.5003 F'(X)=146.50
2.3281
7
0.35215
7.35215
9.04219
16.3943
Iteracin 6:
2
-4
3.96919
2 -0.0308
3.96919
2 3.93839
7
-0.0611
6.93886
7.81611
14.755
-14
6
-10
15.3253 2.76265
18.2655
1.32535 8.76265 F(X)=8.26549
34.1735 73.9964
35.4989 F'(X)=82.75
9
X = 2.08447 - 8.26549 / 82.759 = 1.98459655
6
-14
-0.4548
13.7708
5.54521
-0.2292
29.2827 57.6595
29.0535 F'(X)=63.204
-10
11.005
F(X)=1.005
2.08447
1.9846
7
-0.1233
6.87674
7.62827
14.505
6
-14
-10
-0.9091
13.5382
10.0224
5.09089
-0.4618
F(X)=0.02241
55.3089
28.556
28.0942 F'(X)=60.399
8
73
1.9687
Iteracin 8:
2
-4
3.93665
2 -0.0634
3.93665
2
3.8733
7
-0.1247
6.87531
7.62391
14.4992
6
-14
-10
-0.9195
13.5328
10
5.08047
-0.4672
F(X)=1.2E-05
28.5392 55.2548
28.072 F'(X)=60.335
1.96832
7
-0.1247
6.87531
7.62391
14.4992
6
-14
-10
-0.9195
13.5328
10
5.08046
-0.4672
F(X)=3.4E-12
55.2548
28.5392
28.072 F'(X)=60.335
1.96832
74
Resolviendo el sistema:
CN-1*P + CN-2*Q = B N
-0.040098*P + 3.938624*Q = 1.077989
CN-2*P + CN-3*Q = B N-1
3.938624*P + 0.208443*Q = -0.153693
P = -0.053478159; Q = 0.273152559
Iteracin 2:
P = P + P = -0.121426706; Q = Q + Q = 1.012096162
2 -0.063351 6.8753052
-0.467167468
5.080462767
-2.024192
-0.181673967
-4.931852776
0.2428534 0.0217964
0.591701321
-0.006938336
2 0.1795027 4.8729093 BN-1=-0.057140113 BN=0.141671655
Aplicando por segunda ocasin el mtodo para el nuevo polinomio:
2
0.1795027
4.8729093
-0.057140113
-2.024192
-0.427464972
p=0.121426706
q=1.01209616
2
p=0.121426706
q=1.01209616
2
0.2428534
0.0512853
0.352137721
2 CN-3=0.4223561 CN-2=2.9000023
CN-1=-0.132467365
Resolviendo el sistema:
CN-1*P + CN-2*Q = B N
-0.132467365*P + 2.9000023*Q = 0.141671655
CN-2*P + CN-3*Q = B N-1
2.9000023*P + 0.4223561*Q = -0.057140113
P = -0.026641078; Q = 0.047635336
Iteracin 3:
P = P + P = -0.148067784; Q = Q + Q = 1.059731498
2 -0.063351 6.8753052
-0.467167468
5.080462767
-2.119463
-0.246689416
-5.076442519
0.2961356 0.0344679
0.709290604
-0.000676119
2 0.2327848 4.7903101 BN-1=-0.00456628 BN=0.003344129
Aplicando por segunda ocasin el mtodo para el nuevo polinomio:
2
0.2327848
4.7903101
-0.00456628
-2.119463
-0.560513606
p=0.148067784
q=1.05973149
8
p=0.148067784
q=1.05973149
8
0.2961356
0.0783161
0.407062502
2 CN-3=0.5289204 CN-2=2.7491632
CN-1=-0.158017383
Resolviendo el sistema:
CN-1*P + CN-2*Q = B N
-0.158017383*P + 2.7491632*Q = 0.003344129
CN-2*P + CN-3*Q = B N-1
2.7491632*P + 0.5289204*Q = -0.00456628
75
P = -0.001874275; Q = 0.001108687
Iteracin 4:
P = P + P = -0.149942059; Q = Q + Q = 1.060840185
2 -0.063351 6.8753052
-0.467167468
5.080462767
-2.12168
-0.250924113
-5.080460302
p=0.14994206
q=1.06084018
5
0.2998841 0.0354663
0.718086181
-8.09751E-07
2 0.2365334 4.7890911 BN-1=-5.40042E-06 BN=1.65473E-06
Aplicando por segunda ocasin el mtodo para el nuevo polinomio:
2
0.2365334
4.7890911
-5.40042E-06
p=-0.14994206
-2.12168
-0.569053236
q=1.06084018
5
0.2998841
0.0804315
0.412017129
2
CNCN-2=2.7478423
CN-1=-0.157041507
=0.5364175
3
Resolviendo el sistema:
CN-1*P + CN-2*Q = B N
-0.157041507*P + 2.7478423*Q = 1.65473E-06
CN-2*P + CN-3*Q = B N-1
2.7478423*P + 0.5364175*Q = -5.40042E-06
P = -2.05991E-06; Q = 4.84469E-07
Iteracin 5:
P = P + P = -0.149944119; Q = Q + Q = 1.06084067
2 -0.063351 6.8753052
-0.467167468
5.080462767
-2.121681
-0.250928599
-5.080462767
0.2998882 0.0354674
0.718096066
-6.6472E-14
0.2365375 4.7890912 BN-1=-4.43312E-13 BN-1=1.98966E-12
p=0.14994412
q=1.06084067
0.2365375
4.7890912
-2.121681
0.0804339
CN-2=2.7478438
-4.43312E-13
-0.569062237
0.412023016
CN-1=-0.157039221
p=-0.14994412
q=1.06084067
0.2998882
2 CN-3=0.5364257
Resolviendo el sistema:
CN-1*P + CN-2*Q = B N
-0.157039221*P + 2.7478438*Q = 1.98966E-12
2.7478438*P + 0.5364257*Q = -4.43312E-13
CN-2*P + CN-3*Q = B N-1
P = -2.99344E-13; Q = 7.06973E-13
P = P + P = -0.149944119; Q = Q + Q = 1.06084067
76
77
% raices3.m
% Calculo de raices reales o complejas de funciones polinomicas.
% FUNCION PRINCIPAL:
global A N
nf = 'HORNER';
[N, A] = LeePolinomio;
% Obtener una copia de A:
B = A;
NB = N;
li = input('Ingrese el lmite izquierdo: ');
ld = input('Ingrese el lmite derecho: ');
dx = input('Ingrese el incremento de variacin: ');
xi = li;
78
while 1
[p, xi, xd] = Localiza_Raiz(nf, xi, dx, ld);
if (p < 0)
Raiz = Muller(nf, xi, xd);
% DEFLACIONAR EL POLINOMIO: */
for i = 2 : N
A(i) = A(i) + A(i - 1) * Raiz;
end;
N = N - 1;
end
xi = xd;
xd = xd + dx;
if (xd > ld) break; end;
end;
if mod(N, 2) == 1
disp('Error: NO SE CALCULARON TODAS LAS RAICES REALES');
return;
end;
if N > 1
while N > 2
NewtonBairstow;
end;
% POLINOMIO RESTANTE, VERIFICAR SI LAS RAICES SON REALES O
COMPLEJAS */
P1 = -A(2) / 2 / A(1); P2 = A(2) * A(2) - 4 * A(1) * A(3);
if (P2 < 0)
disp(sprintf('RAICES COMPLEJAS = %f +/- %f * I\n', P1, sqrt(-P2)/2/A(1)));
else
disp(sprintf('RAICES REALES MULTIPLES O CERCANAS: %f Y %f\n', P1 sqrt(P2)/2/A(1), P1 + sqrt(P2)/2/A(1)));
end;
end;
fprintf('USO DEL METODO MATLAB:');
roots(B)
%NewtonBairstow.m
function NewtonBairstow
global A N;
NIter = 100; Precision = 1e-7;
% CALCULAR LAS COTAS DEL MODULO DE LA RAIZ */
AMax = abs(A(1));
for I = 2 : N + 1
if AMax < abs(A(I))
AMax = abs(A(I));
end;
end;
79
80
disp(sprintf('RAICES REALES MULTIPLES O CERCANAS: %f Y %f\n', (-P1 sqrt(Q1)) / 2, (-P1 + sqrt(Q1)) / 2));
end;
disp(sprintf('EN %d NITERACIONES\n', iter));
% TRANSFERIR LOS COEFICIENTES B(I) A LOS A(I): */
for I = 1 : N - 1
A(I) = B(I);
end;
N = N - 2;
return;
81
Como
se
puede
apreciar, existen tres
cuadros de edicin,
mismas que permitirn
ingresar los datos; un
cuadro de lista, el cual
permitir presentar los
resultados; y cuatro
botones, mismos que
tienen la funcionalidad:
- Lectura Proceso:
Valida las entradas;
abre una ventana
para
registrar
el
nombre del archivo el
cual permite leer los
datos del polinomio
desde el medio de
almacenamiento
permanente; ejecuta
los mtodos para
calcular las races
reales y/o complejas y despliega los resultados.
- Guardar Como : Abre una ventana para registrar el nombre del archivo en el
que se almacena permanentemente los resultados obtenidos en el cuadro lista.
- Limpiar: deja en blanco el cuadro lista; y
- Salir: Finaliza la ejecucin del programa.
Se describe a continuacin el cdigo fuente:
Botn Lectura Proceso:
private void b_lectura_Click(object sender, EventArgs e)
{
double liminf, limsup, delta, estado, xi, xd = 0, Raiz;
int I;
try
{
// Transferir valores de los texbox a memoria:
liminf = double.Parse(txtInicial.Text);
limsup = double.Parse(txtFinal.Text);
delta = double.Parse(txtIncremento.Text);
// Abrir el archivo de entrada a traves del cuadro de dialogo:
CEntSalArchivos objes = new CEntSalArchivos("c:\\");
objes.AbrirArchivoEnt();
// Leer el grado del polinomio:
CHorner.n = Int16.Parse(objes.LeerLinea());
if (CHorner.n > 0)
{
82
83
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
Las
referencias
a
las
clases
CEntSalArchivos,
CRaicesReales,
CRaicesRealesComplejas, CHorner, CLectEscrMatrices y CListBox, provienen de la
84
85
2.-. Para la ecuacin y rango del ejercicio (1), localizar los cambios de signo para un
incremento en las abscisas de uno.
3.-. Calcular una raz real para cualquiera de los cambios de signo del ejercicio 2, por
el mtodo de la biseccin usando una precisin de 1E-1. Indicar en cuantas
iteraciones se produce la convergencia.
86
4.-. Calcular la raz real para cualquiera de los cambios de signo del ejercicio 2, por el
mtodo de aproximaciones sucecivas usando una precisin de 1E-1. Indicar en
cuantas iteraciones se produce la convergencia.
87
5.-. Calcular una raz real para el cambio de signo del ejercicio 4, por el mtodo de
aproximaciones sucesivas modificado usando una precisin de 1E-2. Indicar en
cuantas iteraciones se produce la convergencia.
6.-. Calcular una raz real para el cambio de signo del ejercicio 5, por el mtodo de
regula - falsi usando una precisin de 1E-1. Indicar en cuantas iteraciones se
produce la convergencia.
88
7.-. Calcular una raz real para el cambio de signo del ejercicio 6, por el mtodo de la
secante usando una precisin de 1E-1. Indicar en cuantas iteraciones se produce la
convergencia.
89
8.-. Calcular una raz real para el cambio de signo del ejercicio 7, por el mtodo de
Newton - Raphson usando una precisin de 1E-2. Indicar en cuantas iteraciones se
produce la convergencia.
90
9.-. Calcular una raz real para el cambio de signo del ejercicio 8, por el mtodo de
Muller usando una precisin de 1E-2. Indicar en cuantas iteraciones se produce la
convergencia.
10.-. Indicar una ventaja y una desventaja de los mtodos para calcular races reales
(Biseccin, Secante, Newton Raphson y Muller). A su criterio cual mtodo le parece
el ms eficiente y el menos eficiente; por que?.
.
91
13.- Calcular una raz real para cualquiera de los cambios de signo del ejercicio 12,
por el mtodo de Newton - Raphson - Horner usando una precisin de 1E-2. Indicar
en cuantas iteraciones se produce la convergencia.
92
93
15.- Desarrollar un programa en Matlab y otro Visual C#.NET para ingresar desde la
interface de usuario un intervalo dado por XI y XD y desde archivos, los coeficientes
de los polinomios del problema del numeral 14. A continuacin el programa debe
buscar a) Una discontinuidad, en [XI, XD], en el caso de existir, debe desplegar el
mensaje de error apropiado y finalizar. b) localizar la primera raz real en [XI, XD] y
calcularla con una precisin de 1E-8 por el mtodo de Horner Muller y desplegarla.
c) en el caso de no existir la raz en [XI, XD], desplegar el mensaje apropiado y
finalizar. Ejecutar los programas con un grupo de datos que coincidan con los casos
(a), (b) y (c).
16.- Considerar la siguiente ecuacin F(x) = 3x4 + 4x3 + 2x2 - 6x + 5 = 0.Calcular las
dos primeras races complejas por el mtodo de Newton - Bairstow, empleando una
presicin de 1e-5. Indicar en cuantas iteraciones se produce la convergencia.
94
17.- Ejecutar los programas en Matlab y Visual C#.NET para el clculo de races
complejas, ingresando los datos del polinomio del numeral 16.
95
96
3
ANLISIS MATRICIAL
97
98
CAPITULO III
ANALISIS MATRICIAL
Generalidades
En este captulo, se analizarn los siguientes temas:
Generalidades.
Resolucin de sistemas de ecuaciones lineales y simultneas.
Clculo del determinante de una matriz.
Clculo de la matriz inversa.
Caso particular en la resolucin de sistemas de ecuaciones lineales y
simultneas: sistemas simtricos.
Caso particular en la resolucin de sistemas de ecuaciones lineales y
simultneas: sistemas simtricos en banda.
Caso particular en la resolucin de sistemas de ecuaciones lineales y
simultneas: sistemas simtricos del tipo Sky - Line.
Ecuacin Lineal:
Es una secuencia de trminos cuyas variables o incgnitas estn elevados a la
potencia uno, cero, y no estn afectadas por otras variables o incgnitas.
Ejemplos:
2x + 3y 2z = 37 es una ecuacin lineal; -8x + 2xz = 3 no es una ecuacin lineal.
Sistema de Ecuaciones Lineales y Simultneas:
Es una agrupacin de ecuaciones lineales para las que el nmero de ecuaciones es
igual al nmero de incgnitas
Ejemplo:
2x 3y 2z = 37
-18x 16 y + 12z=18
2y 4z = 3
Es un sistema lineal de tres ecuaciones con tres incgnitas.
99
2x y = 4
Al sistema :
, se lo puede plantear :
x 2y = 4
2
Haciendo : A =
1
1 x 4
=
2 y 4
1
x
4
; X = ; B =
2
y
4
(3.1)
a12
a
a
1 a
Si A = 11 12 , entonces A 1 = 22
D a 21 a11
a 21 a 22
simultneas.
100
(3.3)
4
2 1
1 - 2
y B =
, entonces A -1 = -
Dado que A =
3 - 1
4
1 2
Aplicando la relacin (2) :
x
1 - 2
X = = -
3 - 1
y
1 4
1 - 12 4
= -
=
2 4
3 12 4
Mtodo de Gauss
En este mtodo, se puede calcular el determinante de la matriz de coeficientes; se
aplican dos fases:
Eliminacin Gaussiana y
Retrosustitucin.
Fase de Eliminacin Gaussiana:
Consiste en transformar la matriz de coeficientes aplicando transformaciones
elementales, generndose una matriz triangular superior.
Fase de Retrosustitucin
Una vez realizada la eliminacin gaussiana, se despejan desde la ltima ecuacin, la
ltima incgnita, seguidamente de la penltima ecuacin se despeja la penltima
incgnita y se sustituye el valor de la ltima incgnita, y as sucesivamente hasta
llegar a la primera ecuacin.
101
Ejemplo 3.1.
Dado el sistema:
a11 x1 + a12 x2 + a13 x3 = b1
a21 x1 + a22 x2 + a23 x3 = b2
a31 x1 + a32 x2 + a33 x3 = b3
Resolverlo por el mtodo de Gauss y calcular el determinante de la matriz de
coeficientes:
Solucin:
Deben darse los siguientes pasos:
- Planteamiento matricial
a 21 a 22 a 23 x 2 = b 2
a
31 a 32 a 33 x 3 b 3
- Eliminacin Gaussiana:
Por facilidad, se designa a la fila i como: FI. En la primera etapa, se aplican las
siguientes transformaciones elementales:
m1=-a21/a11; F2 F2 + m1 * F1; m2=-a31/a11; F3 F3 + m2 * F1
a21 = a21 + (-a21 / a11)a11 = 0; a22 = a22 + m1*a12 = a22(1) ; b2 = b2 + m1* b1 = b2 (1)
a31 = a31 + (-a31 / a11)a11 = 0; a32 = a32 + m2*a12 = a32(1) ; b3 = b3 + m2* b1 = b2 (1)
Sustituyendo estos transformaciones en el sistema:
a x b
a a
11 12(1) 13(1) 1 1(1)
(k )
0 a 22 a 23 x 2 = b 2 ; a ij = elemento de la fila i, columna j, etapa de clculo k
0 a (1) a (1) x b (1)
32
33 3 3
0 0 a ( 2) x b ( 2)
33
3 3
102
2 - 1 3
- 1 4 2
3 2 - 2
1 1 - 1
1 x1 1
1 x 2 2
=
- 1 x 3 3
2 x 4 - 2
- Eliminacin Gaussiana:
En la primera etapa, se aplican las siguientes transformaciones elementales:
m1 =1/2; F2 F2 + m1 * F1; m2 = -3/2; F3 F3 + m2 * F1; m3 = -1/2; F4 F4 + m3 * F1
a21 = -1 + (1/2)*2 = 0; a22 = 4 + (1/2)*(-1) = 7/2; a23 = 2 + (1/2)*3 = 7/2;
a24 = 1 + (1/2)*1 = 3/2; b2 = 2 + (1/2)*1 = 5/2;
a31 = 3 + (-3/2)*2 = 0; a32 = 2 + (-3/2)*(-1) = 7/2; a33 = -2 + (-3/2)*3 = -13/2;
a34 = -1 + (-3/2)*1 = -5/2; b3 = 3 + (-3/2)*1 = 3/2;
a41 = 1 + (-1/2)*2 = 0; a42 = 1 + (-1/2)*(-1) = 3/2; a43 = -1 + (-1/2)*3 = -5/2;
a44 = 2 + (-1/2)*1 = 3/2; b4 = -2 + (-1/2)*1 = -5/2;
103
0
0
1 x1 1
7/2 7/2 3/2 x 2 5/2
=
7/2 - 13/2 - 5/2 x 3 3/2
3/2 - 5/2 3/2 x 4 - 5/2
-1
3
2 - 1
0 7/2 7/2
0 0 - 10
0 0 - 4
1 x1 1
3/2 x 2 5/2
=
- 4 x 3 - 1
6/7 x 4 - 25/7
3
2 - 1
0 7/2 7/2
0 0 - 10
0 0
0
111/35;
1 x1 1
3/2 x 2 5/2
=
- 4 x 3 - 1
86/35 x 4 - 111/35
104
x1
47
x
2
1 56
X=
=
x 3 86 53
- 111
x 4
1 1
a 21 a 22 a 23 ... a 2k ... a 2j ... a 2N x 2 b 2
a 31 a 32 a 33 ... a 3k ... a 3j ... a 3N x 3 b 3
... ...
...
a k1 a k2 a k3 ... a kk ... a kj ... a kN x k = b k
...
... ...
a i1 a i2 a i3 ... a ik ... a ij ... a iN x i b i
...
... ...
x b
a N1 a N2 a N3 ... a Nk ... a Nj ... a NN N N
105
a 11
106
Fase de retrosustitucin :
b k = (bk -
b )/akk (3.8)
kj j
j = k +1
k = N, N - 1, ..., 3, 2, 1
Algoritmo para el mtodo de Gauss:
a)
b)
c)
107
Planteamiento matricial
a11 a12
a 21 a 22
a
31 a 32
a13 x1 b1
a 23 x 2 = b 2
a 33 x 3 b 3
(1)
0 a 32 a 33(1) x 3 b3(1)
1 0 a13( 2) x1 b1( 2)
0 1 a 23( 2) x 2 = b 2 ( 2)
0 0 a 33( 2) x 3 b3( 2)
108
( 3)
0 1 0 x 2 = b 2
0 0 1 x b (3)
3 3
Ejemplo 3.4.
Resolver el sistema de ecuaciones lineales y simultneas, propuesto en ejemplo 3.3.
Expresar las operaciones y resultados en nmeros racionales.
Solucin:
-Planteamiento matricial:
2 - 1 3
- 1 4 2
3 2 - 2
1 1 - 1
1 x 1 1
1 x 2 2
=
- 1 x 3 3
2 x 4 - 2
109
1
0
0
0
1
0
0
0
0 2 5/7 x1 6/7
1 1 3/7 x 2 = 5/7
0 - 10 - 4 x 3 - 1
0 - 4 6/7 x 4 - 25/7
1
0
0
0
0
1
0
0
0 - 3/35 x1 23/35
0 1/35 x 2 = 43/70
1 2/5 x 3 1/10
0 86/35 x 4 - 111/35
110
x1
47
x
2
1 56
X=
=
x 3 86 53
- 111
x 4
1 1
a 21 a 22 a 23 ... a 2k ... a 2j ... a 2N x 2 b 2
a 31 a 32 a 33 ... a 3k ... a 3j ... a 3N x 3 b 3
... ...
...
a k1 a k2 a k3 ... a kk ... a kj ... a kN x k = b k
...
... ...
a i1 a i2 a i3 ... a ik ... a ij ... a iN x i b i
...
... ...
x b
a N1 a N2 a N3 ... a Nk ... a Nj ... a NN N N
111
1
0
0
m = akk (3.11); akj = akj / m (3.12); aij = aij - aik * akj (3.13); bi = bi - aik * bk (3.14)
k = 1, 2, 3, ..., N ; i = 1, 2, 3, ..., N, i k; j = k, k + 1, k + 2, ..., N;
N
112
(4)
(5)
(6)
x1
0
x1(1)
x1(2)
...
x1(k)
x2
0
x2(1)
x2(2)
...
x2(k)
x3
0
x3(1)
x3(2)
...
x3(k)
MAX(|xi(k+1) - xik |)
valor1
Valor2
...
valork
113
Se analiza el sistema:
x - 2y = -4 (1)
2x - y = 4 (2)
Despejando x de (1) e y de (2): x = 2y - 4; y = 2x - 4 , de esta manera se genera la
tabla:
ITERACION
0
1
2
3
4
5
6
7
8
x
0
-4
-28
-124
-508
-2044
-8188
-32764
-131068
y
0
-12
-60
-252
-1020
-4092
-16380
-65532
-262140
MAX(|xi
(k+1)
- xi |)
24
96
384
1536
6144
24576
98304
393216
x
0
2
3.5
3.875
3.96875
3.9921875
3.998046875
3.999511719
3.99987793
3.999992371
3.999998093
3.999999523
3.999999881
3.99999997
3.999999993
(k+1)
k
y
MAX(|xi
- xi |)
0
3
3
3.75
0.75
3.9375
0.1875
3.984375
0.046875
3.99609375
0.01171875
3.999023438
0.002929688
3.999755859
0.000732422
3.999938965
0.000183105
3.999984741
4.57764E-05
3.999996185
1.14441E-05
3.999999046
2.86102E-06
3.999999762
7.15256E-07
3.99999994
1.78814E-07
3.999999985
4.47035E-08
114
(1)
(2)
(3)
(4)
(k+1)
k
x3
x4
MAX(|xi
- xi |)
0
0
1.91666667
0.083333333 1.91666667
-0.194444444 0.972222222 1.425925926
1.50925926
0.67592593
1.077160494 0.225308642 -0.08436214
1.51028807
1.69341564
0.229080933 1.022805213 1.086019662
1.1703818
1.08293324
x1
0
1
x2
0
0.75
115
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
116
1.35315374
1.31293008
1.23363384
1.35016091
1.24488382
1.31861066
1.27825804
1.29266995
1.29420353
1.28545377
1.29546519
1.28734317
-1.2925444
1.29002879
1.29063895
1.29108849
1.29024585
1.29106616
1.29046056
-1.2908132
1.29066882
1.29067805
1.29073551
1.29066057
-1.2907252
1.29068156
1.29070434
1.29069717
0.71314713
0.3868842
0.17730249
0.11652706
0.13366319
0.12165325
0.08566686
0.04720858
0.02580722
0.00975889
0.01038997
0.01152228
0.00940577
0.00605741
0.00327304
0.00163222
0.00084263
0.00096408
0.00094661
0.00070279
0.00041179
0.00022403
9.5773E-05
7.4942E-05
8.6092E-05
7.475E-05
5.0751E-05
2.7325E-05
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
1.29069512
1.29070105
1.29069475
1.29069966
1.29069663
1.29069801
1.29069776
1.29069741
1.29069796
1.29069745
1.29069781
1.29069761
1.29069769
1.29069769
1.29069765
1.4647E-05
5.9304E-06
6.7628E-06
7.2658E-06
5.698E-06
3.5334E-06
1.9143E-06
9.0147E-07
5.5042E-07
6.3091E-07
5.8686E-07
4.1987E-07
2.3586E-07
1.2874E-07
5.0702E-08
117
a x
kj j
(3.21)
j = 1, j k
118
Mtodo de Crout
Consiste en sustituir la matriz de coeficientes A, por el producto matricial LU y el
vector de trminos independientes B por el producto LC.
Entonces: A = LU (3.31); B = LC (3.32), donde L, es una matriz triangular inferior de
11 0
21 22
31 32
L =
i1 i2
N1 N2
la forma:
0 ... 0 ...
0 ... 0 ...
...
33 ... 0
...
...
i3 ... ii
...
N3 ... Ni ...
0
0
0
;
0
NN
ij
0, para i < j
: =
0, para i j
1 u
0 112
0 0
U =
0 0
0 0
...
cN
0 ... 0 ... 0 ... 1
aij = ikukj
(3.33)
k =1
frmula:
Clculo de la matriz L:
De (3.33) se tiene que:
j -1
ik kj
(3.34)
k = j +1
Pero:
ujj = 1 (3.35);
u =
ik kj
i(j + 1) (j + 1)j
i(j + 2) (j + 2)j
k = j +1
119
+ ... +
u =0
iN Nj
(3.36)
j -1
ij
= aij - ikukj
(3.37)
k =1
= ai1 (3.38)
Clculo de la matriz U:
De (3.33) se tiene que:
i -1
u +
aij =
u +
ii ij
k =1
Pero:
ik kj
(3.39)
ik kj
k = i +1
u =
ik kj
i(i + 1) (i + 1)j
i(i + 2 ) (i + 2 )j
+ ... +
uNj = 0 (3.40)
iN
k = i +1
uij = (aij -
u )/
ik kj
(3.41)
ii
k =1
u1j = a1j /
11
(3.42)
bi = ijcj
(3.43)
j =1
frmula:
bi = ijcj + iici +
j =1
(3.44)
ij j
j = i +1
120
Pero:
N
c=
ij j
i(i + 1) (i + 1)
i(i + 2 ) (i + 2 )
+ ... +
cN = 0 (3.45)
iN
j = i +1
Sustituyendo (3.45) en (3.44) y despejando ci:
i -1
ci = (bi - ijcj)/
(3.46)
ii
j =1
c1 = b1 /
(3.47)
11
(3.48)
(3.49)
(3.50)
(3.51)
ci = uijxj
De (3.53) se tiene que:
i -1
j =1
ci = uijxj + uiixi +
j =1
(3.53)
u x
(3.54)
ij j
j = i +1
121
Pero:
N
u x = u
ij j
i(i + 1) (i + 1)
j = i +1
i -1
u x = 0
ij j
j =1
xi = ci -
u x
ij j
(3.57)
j =i + 1
xN = c N
(3.58)
Ejemplo 3.7.
Dado el sistema:
a11 x1 + a12 x2 + a13 x3 + a14 x4 = b1
a21 x1 + a22 x2 + a23 x3 + a24 x4 = b2
a31 x1 + a32 x2 + a33 x3 + a34 x4 = b3
a41 x1 + a42 x2 + a43 x3 + a44 x4 = b4
Resolverlo por el mtodo de Crout y calcular el determinante de la matriz de
coeficientes:
122
Solucin:
Se plantea el esquema:
a11
l11
a21
l21
a31
l31
a41
l41
a12
u12
a22
l22
a32
l32
a42
l42
a13
u13
a23
u23
a33
l33
a43
l43
a14
u14
a24
u24
a34
u34
a44
u44
b1
c1 x1
b2
c2 x2
b3
c3 x3
b4
c4 x4
43
= a43 -
uk3 = a43 - (
u13 +
4k
41
u 23)
42
k =1
casos:
Es decir, para calcular un elemento de L, a partir de la segunda columna, se resta del
respectivo elemento de A, la suma de los productos de los elementos de L de esa fila
por los de U de esa columna.
Es decir, para calcular un elemento de U, a partir de la segunda fila, se resta del
respectivo elemento de A, la suma de los productos de los elementos de L de esa fila
por los de U de esa columna, este resultado se divide para el elemento de L de la
diagonal principal.
2
u 34 = (a34 -
u )/
3k k4
33
= (a34 - (
31 14 +
u 24))/
32
33
k =1
det =
ii
(3.59)
i =1
123
det =
ii
11
22
33
44
i =1
c3 = (b3 -
c )/
3j j
33
= (b3 - (
c+
31 1
c ))/
32 2
33
j =1
Es decir, para calcular un elemento de C, a partir del segundo, se resta del respectivo
elemento de B, la suma de los productos de los elementos de L de esa fila por los de
C, este resultado se divide para el elemento de L de la diagonal principal.
Se aplica la formula (3.58) para calcular el ltimo elemento de X; y la frmula (3.57)
para calcular los dems elementos de X, partiendo desde el penltimo hasta el
primero. Un clculo tpico sera:
4
x 2 = c2 - u 2jxj = c2 - (u 23 x3 + u 24 x 4)
j =3
Ejemplo 3.8.
Dado el sistema:
-4 1
1 5
-1 -2
2 -4
-1 2
-2 -4
8 3
3 -5
x1
x2
x3
x4
8
-2
4
-6
124
Solucin:
Se plantea el esquema:
-4
-4
1
1
1
-1/4
5
21/4
-1
1/4
-2
-3/7
2
-1/2
-4
-2/3
-1
-1
-2
-9/4
8
51/7
3
7/51
2
2
-4
-7/2
3
1
-5
-110/17
8
-2-178/99
-2
0164/495
4
14/51112/49
5
-6
58/16558/16
5
125
(3.60)
k =1
Clculo de la matriz U:
a1j = a1j/a11
(3.61);
i -1
(3.62)
k =1
b1 = b1/a11
i -1
(3.64)
j =1
bi = bi -
a xb
ij
(3.65)
j =i + 1
126
Para cada lenguaje, se elabora un solo programa, desde el cual se abren archivos
para lectura de datos; se ejecuta cada uno de los mtodos de resolucin de sistemas
de ecuaciones lineales y simultneas y se escribe los resultados en un archivo de
salida, en el caso de Matlab y se despliegan los resultados en una caja lista en el
caso de Visual C#.
Archivos fuente en MATLAB:
%SisEcuac1.m
% APERTURA DE ARCHIVOS:
ent = fopen('D:\AnalNumat\matriz.dat', 'r');
sal = fopen('D:\AnalNumat\rsisecuac1.dat', 'w');
% LECTURA DE LA DIMENSION DEL SISTEMA:
N = fscanf(ent, '%d', 1);
A = LeeMatriz(ent, N);
B = LeeVector(ent, N);
fprintf(sal, 'RESOLUCION DE SISTEMAS DE ECUACIONES LINEALES Y SIMULTANEAS:\r\n');
fprintf(sal, 'MATRIZ DE COEFICIENTES:\r\n');
EscribeMatriz(sal, A, N);
fprintf(sal, 'VECTOR DE TERMINOS INDEPENDIENTES:\r\n');
EscribeVector(sal, B, N);
[det, X] = Gauss(A, B, N);
fprintf(sal, 'METODO DE GAUSS: VECTOR DE SOLUCIONES:\r\n');
EscribeVector(sal, X, N);
fprintf(sal, 'DETERMINANTE = %f\r\n', det);
[det, X] = GaussJordan(A, B, N);
fprintf(sal, 'METODO DE GAUSS - JORDAN: VECTOR DE SOLUCIONES:\r\n');
EscribeVector(sal, X, N);
fprintf(sal, 'DETERMINANTE = %f\r\n', det);
fprintf(sal, 'METODO DE GAUSS - SEIDEL: VECTOR DE SOLUCIONES:\r\n');
[X, estado, msg] = GaussSeidel(A, B, N);
if estado > 0
EscribeVector(sal, X, N);
end;
fprintf(sal, msg);
fprintf(sal, '\r\nUSO DEL METODO MATLAB:\r\n');
Y = A \ B';
fprintf(sal, 'VECTOR DE SOLUCIONES:\r\n');
EscribeVector(sal, Y, N);
fclose(ent);
ent = fopen('D:\AnalNumat\matrizcr.dat', 'r');
N = fscanf(ent, '%d', 1);
NS = fscanf(ent, '%d', 1);
A = LeeMatriz(ent, N);
fprintf(sal, 'RESOLUCION DE SISTEMAS DE ECUACIONES LINEALES Y SIMULTANEAS:\r\n');
fprintf(sal, 'METODO DE CROUT:\r\n');
fprintf(sal, 'MATRIZ DE COEFICIENTES:\r\n');
EscribeMatriz(sal, A, N);
for i = 1 : NS
B = LeeVector(ent, N);
fprintf(sal, 'VECTOR DE TERMINOS INDEPENDIENTES:\r\n');
EscribeVector(sal, B, N);
[det, A, B] = Crout(A, B, N, i);
if i == 1
fprintf(sal, 'DETERMINANTE = %f\r\n', det);
end;
fprintf(sal, 'VECTOR DE SOLUCIONES:\r\n');
EscribeVector(sal, B, N);
end;
127
fclose('all');
% LeeMatriz.m
function A = LeeMatriz(f, m)
for i = 1 : m
for j = 1 : m
A(i, j) = fscanf(f, '%f', 1);
end;
end;
return
% LeeVector.m
function A = LeeVector(f, m)
for i = 1 : m
A(i) = fscanf(f, '%f', 1);
end;
return
% EscribeMatriz.m
function EscribeMatriz(f, A, m)
for i = 1 : m
fprintf(f, 'FILA %d:\r\n',i);
for j = 1 : m
fprintf(f, '%f ',A(i, j));
if (mod(j, 8) == 0) | (j == m)
fprintf(f, '\r\n');
end;
end;
end;
return
% EscribeVector.m
function EscribeVector(f, A, m)
for i = 1 : m
fprintf(f, '%f ',A(i));
if (mod(i, 8) == 0) | (i == m)
fprintf(f, '\r\n');
end;
end;
return
% Gauss.m
function [det, X] = Gauss(A, B, N)
for k = 1 : N - 1
if A(k,k) == 0
% buscar un A(k,k) > 0
for i = k + 1: N
if A(i,k) ~= 0 % intercambiar las filas i con la k:
X = A(i,:); A(i,:) = A(k,:); A(k,:) = X;
det = B(i); B(i) = B(k); B(k) = det;
break;
end;
end;
if A(k,k) == 0
fprintf('Error!, no existe solucin para este sistema');
exit;
end;
end;
% Eliminacin gaussiana:
for i = k + 1 : N
det = -A(i,k)/A(k,k);
for j = k : N
128
129
return;
end;
end;
for iter = 1 : niter
d = 0;
for i = 1 : N
s = 0;
for j = 1 : N
if i ~= j
s = s + A(i, j) * X(j);
end;
end;
t = (B(i) - s) / A(i, i);
max = abs(t - X(i));
if max > d
d = max;
end;
X(i) = t;
end;
if d <= precision
break;
end;
end;
ST = 1;
if d > precision
[msg, er] = sprintf('Atencin!, no se alcanza la precisin requerida: %f\n', d);
fprintf(msg);
ST = -1;
else
[msg, er] = sprintf('Solucin en %d iteraciones\n', iter);
fprintf(msg);
end;
return;
% crout.m
function [det, A, B] = Crout(A, B, N, St)
det = 1;
if St == 1 % calcular las matrices L, U y el determinante:
for i = 1 : N
for j = 2 : N
l = j - 1;
if i < j l = i - 1; end;
for k = 1 : l
A(i, j) = A(i, j) - A(i, k) * A(k,j);
end;
if i < j A(i, j) = A(i, j) / A(i, i); end;
end;
det = det * A(i,i);
end;
end;
% Clculo del vector C:
for i = 1 : N
for j = 1 : i - 1
B(i) = B(i) - A(i, j) * B(j);
end;
B(i) = B(i) / A(i, i);
end;
% Clculo del vector X:
for i = N - 1 : -1 : 1
for j = i + 1 : N
B(i) = B(i) - A(i, j) * B(j);
130
end;
end;
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa, como se haba indicado anteriormente, la entrada salida, se realiza
sobre archivos, cuyos contenidos son:
matriz.dat:
3
3
1
2
1
-4
-1
2
-1
4
1
2
3
matrizcr.dat:
3
2
3
1
2
1
-4
-1
2
-1
4
1
2
3
3
2
1
rsisecuac1.dat
RESOLUCION DE SISTEMAS DE ECUACIONES LINEALES Y SIMULTANEAS:
MATRIZ DE COEFICIENTES:
FILA 1:
3.000000 1.000000 2.000000
FILA 2:
131
132
Lectura Proceso: abre una ventana para registrar el nombre del archivo el cual
permite leer los datos del sistema de ecuaciones lineales y simultneas desde el
medio de almacenamiento permanente; ejecuta los mtodos para calcular las
soluciones y despliega los resultados.
Guardar Como : Abre una ventana para registrar el nombre del archivo en el
que se almacena permanentemente los resultados obtenidos en el cuadro lista.
Limpiar: deja en blanco el cuadro lista; y
Salir: Finaliza la ejecucin del programa.
133
134
}
}
}
MessageBox.Show(ex.Message);
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
135
136
137
138
MATRIZ DE COEFICIENTES:
FILA 1
312
FILA 2
1 -4 -1
FILA 3
2 -1 4
VECTOR DE TERMINOS INDEPENDIENTES:
123
VECTOR SOLUCION:
0.186046541943642 -0.58139532685865 0.511627897313517
139
MATRIZ INVERSA
La matriz inversa de una matriz Anxn se denota por A-1 tambin de nxn en donde
el producto de A x A-1 = A-1 x A = I (Matriz Identidad).
Deduccin de frmulas:
A x A-1 = A-1 x A = I (3.71)
Sea: A-1 = X (3.72)
Entonces, para calcular la matriz inversa, se necesita resolver el sistema:
AX = I
(3.73)
(3.76)
I2
1
0
0
:
0
0
1
0
:
0
140
IN
0
0
0
:
1
Ejemplo 3.9.
Dada la matriz:
A=
-4
2
-1
1
5
-1
-2
8
calcular A-1
Solucin:
Se aplica el mtodo de Crout, por ser el ms apropiado para resolver varios sistemas
de ecuaciones a la vez.
Los resultados se muestran a continuacin:
141
142
143
144
if (CResolSistemaEcuaciones.n > 0)
{
// Definir la matriz de coeficientes:
CResolSistemaEcuaciones.A = new double[CResolSistemaEcuaciones.n][];
for (int i = 0; i < CResolSistemaEcuaciones.n; i++)
CResolSistemaEcuaciones.A[i] = new double[CResolSistemaEcuaciones.n];
// Leer la matriz de coeficientes y generar los resultados en el listbox:
// Clculo de la matriz inversa por el mtodo de Crout:
CLectEscrMatrices.n = CResolSistemaEcuaciones.n;
CLectEscrMatrices.a = CResolSistemaEcuaciones.A;
CLectEscrMatrices.obent = objes;
CLectEscrMatrices.lb = listBox1;
CLectEscrMatrices.LeeMatriz();
listBox1.Items.Add("CALCULO DE LA MATRIZ INVERSA POR EL METODO DE CROUT");
CLectEscrMatrices.EscribeMatriz(" DE COEFICIENTES:");
CResolSistemaEcuaciones.MatrizInversa();
CLectEscrMatrices.EscribeMatriz(" INVERSA:");
// Cerrar flujo de entrada:
objes.CerrarLectura();
listBox1.Items.Add("");
listBox1.Items.Add("");
listBox1.Items.Add("");
}
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
145
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
146
147
aij = aji
i = 1, 2, 3,,N
j = 1, 2, 3,,N
Con el fin de sistematizar el mtodo para resolver este tipo de ecuaciones y optimizar
tiempo de proceso y espacio de memoria, se realiza las siguientes simplificaciones:
uij =
ji
(3.81)
ii
j -1
ij
= aij -
ik jk
(3.82)
kk
k =1
148
xN = cN
(3.87)
= ai1 (3.83)
cj)/
ij
(3.88)
ii
(3.84)
j =1
c1 = b1 /
(3.85)
11
N
xi = ci - (
x )/
ji j
ii
(3.86)
j =i + 1
Algoritmo para el
mtodo de Crout
Simtrico:
1) Aplicando
la
frmula 3.83 se
determinan los
valores de la primera columna de L.
2) Aplicando la frmula 3.82 se calcula la siguiente columna de la matriz L.
3) Repetir el paso 2 hasta calcular todas las columnas de la matriz L.
4) Aplicando la frmula 3.85, se calcula el primer elemento del vector C.
5) Aplicando la frmula 3.84 se calcula el siguiente elemento del vector C.
6) Repetir el paso 5 hasta calcular todos los elementos del vector C.
7) Aplicando la frmula 3.87 se determina el ltimo elemento del vector X.
8) Aplicando la frmula 3.86 se calcula el elemento anterior del vector X.
9) Repetir el paso 8 hasta calcular todos los elementos del vector X.
149
Ejemplo 3.10:
Resolver el sistema de ecuaciones lineales y simultneas del ejemplo 3.8.
Solucin:
Se plantea el esquema:
-4
-4
1
1
-1
-1
2
2
5
21/4
-2
-9/4
-4
-7/2
-2
8
51/7
3
1
0
-5
-110/17
150
14/51
58/465
8
-2
4
-6
-178/99
164/495
112/495
58/165
151
152
matrizsm.dat:
2
4
-4
1
5
-1
-2
8
2
-4
3
-5
8
-2
4
-6
4
3
2
1
rsisecuac2.dat:
RESOLUCION DE SISTEMAS SIMETRICOS DE ECUACIONES LINEALES Y
SIMULTANEAS:
METODO DE CROUT SIMETRICO:
MATRIZ DE COEFICIENTES:
FILA 1:
-4.000000 1.000000 -1.000000 2.000000
FILA 2:
1.000000 5.000000 -2.000000 -4.000000
FILA 3:
-1.000000 -2.000000 8.000000 3.000000
FILA 4:
2.000000 -4.000000 3.000000 -5.000000
VECTOR DE TERMINOS INDEPENDIENTES:
8.000000 -2.000000 4.000000 -6.000000
DETERMINANTE = 990.000000
VECTOR DE SOLUCIONES:
-1.797980 0.331313 0.226263 0.351515
VECTOR DE TERMINOS INDEPENDIENTES:
4.000000 3.000000 2.000000 1.000000
VECTOR DE SOLUCIONES:
-1.424242 0.424242 0.484848 -0.818182
153
154
double det;
int n = CResolSisSimetricosEcuaciones.n;
// Definir el arreglo de trminos independientes:
CResolSisSimetricosEcuaciones.B = new double[n];
// Definir la matriz de coeficientes:
CResolSisSimetricosEcuaciones.A = new double[n * (n + 1) / 2];
/* Leer la matriz de coeficientes, el arreglo de trminos independientes y generar los resultados en
el listbox: */
// Prueba del mtodo de Crout:
CLectEscrMatrices.n = n;
CLectEscrMatrices.c = CResolSisSimetricosEcuaciones.A;
CLectEscrMatrices.b = CResolSisSimetricosEcuaciones.B;
CLectEscrMatrices.obent = objes;
CLectEscrMatrices.lb = listBox1;
CLectEscrMatrices.LeeMatrizSimetrica();
listBox1.Items.Add("RESOLUCION DE SISTEMAS SIMETRICOS DE ECUACIONES LINEALES Y
SIMULTANEAS");
listBox1.Items.Add("POR EL METODO DE CROUT");
CLectEscrMatrices.EscribeMatrizSimetrica(" DE COEFICIENTES:");
for (i = 0; i < ns; i++)
{
CLectEscrMatrices.LeeVector();
CLectEscrMatrices.EscribeVector(" DE TERMINOS INDEPENDIENTES:");
det = CResolSisSimetricosEcuaciones.CroutSimetrico(i == 0);
CLectEscrMatrices.EscribeVector(" SOLUCION:");
if (i == 0) listBox1.Items.Add("DETERMINANTE = " + det.ToString());
}
// Cerrar flujo de entrada:
objes.CerrarLectura();
listBox1.Items.Add("");
listBox1.Items.Add("");
listBox1.Items.Add("");
}
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
Botn Limpiar:
155
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
156
157
(i i + 2 j ) / 2; i m
inb(i, j, m) =
2
(2im 2i + 2 j m + m) / 2; i > m
(3.91)
la funcin:
La cual, transforma los ndices matriciales (i, j) en un ndice lineal inb, convirtiendo de
esta manera la
matriz
de
coeficientes en
un vector. Se
debe tener en
cuenta que i >
0 y j > 0. Para
valores desde
cero pueden
producirse
ndices
repetidos. La
figura 3.7 ilustra esta conversin.
Ejemplo 3.11:
158
Dado el sistema:
6X1 + 3X2 + 2X3 = 4
3X1 + 4X2 + X3 + 2X4 = 3
2X1 + X2 + 5X3 + 3X4 + 4X5 = 1
2X2 + 3X3 + 7X4 + 5X5 + 3X6 = -2
4X3 + 5X4 + 2X5 + 2X6 = 5
3X4 + 2X5 + 8X6 = 6
Resolverlo por el mtodo de Crout y calcular el determinante de la matriz de
coeficientes.
Solucin:
Planteamiento matricial:
3
2
0
0
3
4
1
2
0
0
2
1
5
3
4
0
0
2
3
7
5
3
0
0
4
5
2
2
0 x1 4
0 x 2 3
0 x 3 1
=
3 x 4 2
2
x 5 5
8
x 6 6
159
160
EscribeVector(sal, B, N);
end;
fclose('all');
% Inb.m
% TRANSFORMA LOS INDICES MATRICIALES I,J EN EL INDICE VECTORIAL INB,
% PARA UNA MATRIZ BANDA INFERIOR DE ANCHO M */
function a = Inb(i, j, m)
if (i > m)
a = (m - m ^ 2 + 2 * i * m - 2 * i + 2 * j) / 2;
else
a = (i ^ 2 - i + 2 * j) / 2;
end;
return;
% LeeMatBan.m
function A = LeeMatBan(f, n, m)
for i = 1 : n
l = 1;
if (i > m) l = i - m + 1; end;
for j = l : i
A(Inb(i, j, m)) = fscanf(f, '%f', 1);
end;
end;
return;
% EscribeMatBan.m
function EscribeMatBan(f, A, n, m)
for i = 1 : n
fprintf(f, 'FILA %d:\r\n',i);
for j = 1 : n
fprintf(f, '%f ',RecuBan(A, i, j, m));
if (mod(j, 8) == 0) | (j == n)
fprintf(f, '\r\n');
end;
end;
end;
return
% RecuBan.m
function v = RecuBan(A, i, j, m)
if (j <= i - m) | (j >= i + m)
v = 0;
elseif (i < j)
v = A(Inb(j, i, m));
else
v = A(Inb(i, j, m));
end;
return;
% croutban.m
function [det, A, B] = CroutBan(A, B, N, m, St)
det = 1;
if St == 1 % calcular las matrices L, U y el determinante:
det = A(Inb(1, 1, m));
for i = 2 : N
l = 2;
if (i > m) l = i - m + 2; end;
for j = l : i
for k = l - 1 : j - 1
A(Inb(i, j, m)) = A(Inb(i, j, m)) - A(Inb(i, k, m)) * A(Inb(j, k, m)) / A(Inb(k, k, m));
161
end;
end;
det = det * A(Inb(i,i, m));
end;
end;
% Clculo del vector C:
for i = 1 : N
l = 1;
if (i > m) l = i - m + 1; end;
for j = l : i - 1
B(i) = B(i) - A(Inb(i, j, m)) * B(j);
end;
B(i) = B(i) / A(Inb(i, i, m));
end;
% Clculo del vector X:
for i = N - 1 : -1 : 1
l = m + i - 1;
if (i > N - m) l = N; end;
for j = i + 1 : l
B(i) = B(i) - A(Inb(j, i, m)) * B(j) / A(Inb(i, i, m));
end;
end;
162
3
1
-2
5
6
1
0
0
0
0
0
rsisecuac3.dat:
RESOLUCION DE SISTEMAS SIMETRICOS EN BANDA DE ECUACIONES
LINEALES Y SIMULTANEAS:
MATRIZ DE COEFICIENTES:
FILA 1:
6.000000 3.000000 2.000000 0.000000 0.000000 0.000000
FILA 2:
3.000000 4.000000 1.000000 2.000000 0.000000 0.000000
FILA 3:
2.000000 1.000000 5.000000 3.000000 4.000000 0.000000
FILA 4:
0.000000 2.000000 3.000000 7.000000 5.000000 3.000000
FILA 5:
0.000000 0.000000 4.000000 5.000000 2.000000 2.000000
FILA 6:
0.000000 0.000000 0.000000 3.000000 2.000000 8.000000
VECTOR DE TERMINOS INDEPENDIENTES:
4.000000 3.000000 1.000000 -2.000000 5.000000 6.000000
DETERMINANTE = -3646.000000
VECTOR DE SOLUCIONES:
-0.631377 0.962150 2.450905 -0.702688 -2.211465 1.566374
VECTOR DE TERMINOS INDEPENDIENTES:
1.000000 0.000000 0.000000 0.000000 0.000000 0.000000
VECTOR DE SOLUCIONES:
0.469007 -0.409764 -0.292375 0.262205 0.036753 -0.107515
163
Como se puede
apreciar, existen un
cuadro de lista, el
cual
permitir
presentar
los
resultados; y cuatro
botones, mismos
que
tienen
la
funcionalidad:
- Lectura
Proceso: abre
una
ventana
para registrar el
nombre
del
archivo el cual
permite leer los
datos
del
sistema
simtrico
en
banda
de
ecuaciones lineales y simultneas desde el medio de almacenamiento
permanente; ejecuta los mtodos para calcular la matriz inversa y despliega los
resultados.
- Guardar Como : Abre una ventana para registrar el nombre del archivo en el
que se almacena permanentemente los resultados obtenidos en el cuadro lista.
- Limpiar: deja en blanco el cuadro lista; y
- Salir: Finaliza la ejecucin del programa.
Se describe a continuacin el cdigo fuente:
Botn Lectura Proceso:
private void b_lectura_Click(object sender, EventArgs e)
{
try
{
int i, ns;
// Abrir el archivo de entrada a traves del cuadro de dialogo:
CEntSalArchivos objes = new CEntSalArchivos("c:\\");
objes.AbrirArchivoEnt();
// Leer el nmero de procesamiento de soluciones:
ns = Int16.Parse(objes.LeerLinea());
// Leer el nmero de ecuaciones o de incgnitas:
CResolSisSimBanda.n = Int16.Parse(objes.LeerLinea());
if (CResolSisSimBanda.n > 0)
{
double det;
int n = CResolSisSimBanda.n;
164
Botn Limpiar:
165
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
Las
referencias
a
las
clases
CEntSalArchivos,
CResolSisSimBanda,
CLectEscrMatrices y CListBox, provienen de la definicin de las respectivas clases
incorporadas en la librera del archivo CLibComunClases.cs, el cual se describe en el
Apndice B de este libro.
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
166
167
(3.101)
La cual, transforma los ndices matriciales (i, j) en un ndice lineal insk, convirtiendo
de esta manera la matriz de coeficientes en un vector. Se debe tener en cuenta que i
> 0 y j > 0. Para valores desde cero pueden producirse ndices repetidos.
168
Ejemplo 3.12:
Dado el sistema:
3X1 - 2X2 = -2
-2X1 -3X2 + X3 + 2X4 = 1
X2 - 2X3 + 3X4 = 0
2X2 + 3X3 - X4 = 3
2X5 = 5
-2X6 + 3X7 = 2
3X6 + X7 = -1
Resolverlo por el mtodo de Crout y calcular el determinante de la matriz de
coeficientes.
Solucin:
Planteamiento matricial:
3 - 2 0 0 0 0
- 2 - 3 1 2 0 0
0 1 - 2 3 0 0
0 2 3 - 1 0 0
0 0 0 0 2 0
0 0 0 0 0 2
0 0 0 0 0 3
Como se puede observar,
ancho de banda variable.
0 x1 - 2
0 x 2 1
0 x 3 0
0 x 4 = 3
0 x 5 5
3 x 6 2
1 x 7 - 1
es un sistema simtrico Skyline, de 7 ecuaciones, con un
Se plantea el esquema:
169
170
if i == 1
fprintf(sal, 'DETERMINANTE = %f\r\n', det);
end;
fprintf(sal, 'VECTOR DE SOLUCIONES:\r\n');
EscribeVector(sal, B, N);
end;
fclose('all');
% Insk.m
function a = Insk(i, j, Dr)
a = Dr(i) + j - i;
return;
% LeeMSkyLine.m
function [A, Dir] = LeeMSkyLine(f, n)
PrimerCeroDeFila = 0;
Cont = 0;
for i = 1 : n
for j = 1 : i
Aux = fscanf(f, '%f', 1);
if (i == j)
Cont = Cont + 1;
A(Cont) = Aux;
Dir(i)= Cont;
PrimerCeroDeFila = 1;
elseif (Aux ~= 0)
Cont = Cont + 1;
A(Cont) = Aux;
PrimerCeroDeFila = 0;
elseif (PrimerCeroDeFila == 0)
Cont = Cont + 1;
A(Cont) = Aux;
end;
end;
end;
return;
% EscribeSkyLine.m
function EscribeSkyLine(f, A, n, Dr)
for i = 1 : n
fprintf(f, 'FILA %d:\r\n',i);
for j = 1 : i
fprintf(f, '%f ',RecuSkyLine(A, i, j, Dr));
if (mod(j, 8) == 0) | (j == n)
fprintf(f, '\r\n');
end;
end;
for j = i + 1 : n
fprintf(f, '%f ',RecuSkyLine(A, j, i, Dr));
if (mod(j, 8) == 0) | (j == n)
fprintf(f, '\r\n');
end;
end;
end;
fprintf(f, '\r\nVECTOR DE DIRECCIONES\r\n');
for i = 1 : n
fprintf(f, '%d ', Dr(i));
end;
fprintf(f, '\r\n\r\n');
return
171
% RecuSkyLine.m
function Val = RecuSkyLine(X, I, J, Dir)
if (I + J == 2)| (J > I - Dir(I) + Dir(I - 1))
Val = X(Insk(I, J, Dir));
else
Val = 0;
end;
return;
% croutskyline.m
function [det, A, B] = CroutSkyLine(A, B, N, Dr, St)
det = 1;
if St == 1 % calcular las matrices L, U y el determinante:
det = A(Insk(1, 1, Dr));
for i = 2 : N
lim = i - Dr(i) + Dr(i - 1) + 1;
for j = lim + 1 : i
for k = lim : j - 1
A(Insk(i, j, Dr)) = A(Insk(i, j, Dr)) - A(Insk(i, k, Dr)) * RecuSkyLine(A, j, k, Dr) / A(Insk(k, k, Dr));
end;
end;
det = det * A(Insk(i,i, Dr));
end;
end;
% Clculo del vector C:
for i = 1 : N
d1 = 0;
if i > 1 d1 = Dr(i - 1); end;
for j = i - Dr(i) + d1 + 1 : i - 1
B(i) = B(i) - A(Insk(i, j, Dr)) * B(j);
end;
B(i) = B(i) / A(Insk(i, i, Dr));
end;
% Clculo del vector X:
for i = N - 1 : -1 : 1
for j = i + 1 : N
B(i) = B(i) - RecuSkyLine(A, j, i, Dr) * B(j) / A(Insk(i, i, Dr));
end;
end;
172
1
-2
0
2
3
-1
0
0
0
0
2
0
0
0
0
0
-2
0
0
0
0
0
3
1
-2
1
0
3
5
2
-1
1
0
0
0
0
0
0
rsisecuac4.dat:
RESOLUCION DE SISTEMAS SIMETRICOS SKYLINE DE ECUACIONES
LINEALES Y SIMULTANEAS:
MATRIZ DE COEFICIENTES:
FILA 1:
3.000000 -2.000000 0.000000 0.000000 0.000000 0.000000 0.000000
FILA 2:
173
174
175
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
176
Las
referencias
a
las
clases
CEntSalArchivos,
CResolSisSimSkyline,
CLectEscrMatrices y CListBox, provienen de la definicin de las respectivas clases
incorporadas en la librera del archivo CLibComunClases.cs, el cual se describe en el
Apndice B de este libro.
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
177
178
179
180
181
5.- Ejecutar los programas en Matlab y Visual C#.NET para la resolucin de sistemas
de ecuaciones lineales y simultneas, ingresando los datos del sistema del numeral
1.
182
6.- Dada la matriz de coeficientes del sistema del numeral 1. Calcular la matriz
inversa por el mtodo de Crout. Expresar las operaciones y resultados en nmeros
racionales.
7.- Ejecutar los programas en Matlab y Visual C#.NET para calcular la matriz inversa
por el mtodo de Crout, ingresando los datos del sistema del numeral 6.
183
9.- Ejecutar los programas en Matlab y Visual C#.NET para la resolucin de sistemas
simtricos de ecuaciones lineales y simultneas, ingresando los datos del sistema del
numeral 8.
184
0
0
4 - 3 0 0
0
0
- 3 2 - 1 0
0 - 1 3 4
0
0
0 0 4 5 - 2 0
0 0 0 2 4 2
0 0 0 0
2 3
0
1
0 0 0 0
0 x1 2
0 x 2 - 3
0 x 3 - 1
0 x 4 = 4
0 x 5 1
1 x 6 - 2
6 x 7 3
185
4 - 1
- 1 - 2
0
1
2
0
0
0
0
0
0
0
0
1
-3
0
2
4
0
0
0
4
0
0
-1 3
3 2
0 1
0
0
0
0
1
2
3
0 x1 - 1
0 x 2 2
0 x 3 3
0 x 4 = 1
0 x 5 4
3 x 6 - 2
2 x 7 1
186
187
188
4
INTERPOLACIN
189
190
INTERPOLACION:
Generalidades
La interpolacin permite generar nuevos puntos a partir de un conjuto finito de puntos
conocidos.
Existen varios mtodos para interpolar, los que se tratarn en este captulo son:
Interpolacin Lineal.
Interpolacin de Lagrange.
Interpolacin de Newton.
INTERPOLACION LINEAL:
Consiste en generar el nuevo punto, trazando una secante (segmento AC), sobre la
curva
y = f(x) (ver figura
4.1). Para deducir la frmula
de interpolacin, se establece
la
semejanza
entre
los
tringulos ABC y ADE, de
manera que:
DE BC y y1 y 2 y1
=
=
(4.1)
AD AB x x1 x 2 x1
Despejando de (4.1), el valor
de y:
y=
y 2 y1
( x x1) + y1.; (4.2)
x2 x1
191
Ejemplo 4.1
Dados los puntos A(1; -1) y C(3; 15) tomados de la curva real y = 2x2 3; calcular la
ordenada y; si x = 2.8 adems del error absoluto y relativo.
Solucin:
Aplicando la frmula 4.2:
y 2 y1
15 (1)
( x x1) + y1 =
(2.8 1)) + (1) = 13.4
x 2 x1
3 1
Y = f(2.8) = 2(2.8)2 - 3 = 12.68
EAy = |12.68 - 13.4| = 0.72; y = 0.72 / 12.68 = 0.05678 = 5.68 %
y=
Ejemplo 4.2
Desarrollar una aplicacin Matlab y otra en Visual C#, que tenga una funcin para
calcular la ordenada y por el mtodo de interpolacin lineal; dicha funcin debe
recibir las coordenadas de dos puntos de la curva real y el valor de la abscisa x. la
aplicacin debe calcular tambin los errores absoluto y relativo. Considerar los datos
del ejemplo 4.1.
Solucin:
% PrbInterpolacionLineal.m
% Prueba interpolacin lineal y calcula errores.
x(1) = input('Ingrese la abscisa inicial: ');
while 1
x(2) = input('Ingrese la abscisa final: ');
if x(2) > x(1) break; end;
end;
y(1) = fit(x(1));
y(2) = fit(x(2));
disp(sprintf('y1: %f', y(1)));
disp(sprintf('y2: %f', y(2)));
while 1
xit = input(sprintf('Ingrese la abscisa a interpolar (salir < %f o > %f): ', x(1), x(2)));
% Validacin:
192
193
Como se puede apreciar, existen nueve cuadros de edicin, mismos que permitirn
ingresar los datos; y presentar los resultados; y dos botones, mismos que tienen la
funcionalidad:
- Procesar: Ejecuta el mtodo para interpolar linealmente y desplegar los
resultados. Y
- Salir: Finaliza la ejecucin del programa.
Se describe a continuacin el cdigo fuente:
Funcin Externa:
// Funcion a ser analizada:
static double fit(double x)
{
return 2 * Math.Pow(x, 2) - 3;
}
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
// Prueba interpolacin lineal y calcula errores.
// Declaracion de variables
double x, y, Y, erAbs, erRel;
try
{
// Definicin del vector de abscisas:
194
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Salir:
private void cb_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
195
INTERPOLACION DE LAGRANGE:
Consiste en generar el nuevo punto, a partir de un grupo de puntos (x0; y0), (x1; y1),
(x2; y2), (xn; yn) para lo cual se emplea la frmula de los polinomios de Lagrange:
n
L( x, n) = p j ( x) y j ;4.3
j =0
pj =
i = 0 ,i j
196
x xi
pj es el
;4.4 interpolador;
x j xi
n; n es el grado del
polinomio de interpolacin.
1; j = k
p j ( xk ) =
;4.5
0; j k
La frmula 4.5 garantiza que exista un solo polinomio interpolador y que todos los
puntos dados pasen por la curva real y = f(x).
Ejemplo 4.3.
Obtener el polinomio de intepolacin del problema 4.1 aplicando las frmulas de
interpolacin mediante los polinomios de Lagrange, dados los puntos:
a) A(1, -1) y C(3, 15)
b) A(1, -1), B(2, 5) y C(3, 15)
Solucin:
1
x x0
x x1
a) Dado que L ( x,1) = p j ( x) y j = p 0 ( x) y 0 + p1 ( x) y1 =
y0 +
y1
x
x
x
x
j
=
0
0
1
1
0
se
tienen
dos puntos,
x3
x 1
=
(1) +
(15) = 8 x 9
n=1
2
2
L( x,2) = p j ( x) y j = p 0 ( x) y 0 + p1 ( x) y1 + p 2 ( x) y 2
j =0
x x0 x x1
x x0 x x 2
x x1 x x 2
.
.
.
y2
y1 +
y0 +
x 2 x0 x 2 x1
x1 x0 x1 x 2
x0 x1 x0 x 2
x2 x3
x 1 x 3
x 1 x 2
.
(1) +
.
(5) +
.
(15) = 2 x 2 3
1 2
1
1
2
1
197
Solucin
El proceso se desarrolla en la siguiente tabla:
Punto X
0
6.2831853
1
3.1415926
2
0
3
3.1415927
y = f(X)
x = 1.57079633
Fracciones
Fracciones1 Fracciones2 3
Fracciones4
-1.5
2.5
0.0432139
-0.25
1
0.16666667
23.140693
0.375
-0.5
0.25
0.5
1.5
0.5
0.75
0.75
0.5
1.5
0.5
0.25
-0.5
0.0234375
-0.15625
0.703125
0.46875
-0.0390625
yjpJ(x) = 4.3768E-05
L(/2)
=
9.3740273 y exacta = -4.8104774
-0.0067522
0.703125 10.8471997
-20.917643
0.0018674
Fracciones
0
1.25 0.83333333
0.625
6.2831853 535.49166
pJ(x) =
Error
Error
Abso= 4.5635492 relativo =
0.94866866 = 94.87 %
Ejemplo 4.5
Desarrollar una aplicacin Matlab y otra en Visual C#, que tenga una funcin para
calcular la ordenada y por el mtodo de interpolacin mediante los polinomios de
Lagrange; dicha funcin debe recibir las coordenadas de dos puntos (extremos) de la
curva real; el valor de la abscisa x y el nmero de puntos; la aplicacin debe calcular
tambin los puntos intermedios y los errores absoluto y relativo. Considerar los datos
del ejemplo 4.1.
Solucin:
198
% PrbInterpolacionLagrange.m
% Prueba interpolacin de Lagrange y calcula errores.
x(1) = input('Ingrese la abscisa inicial: ');
while 1
xn = input('Ingrese la abscisa final: ');
if xn > x(1) break; end;
end;
y(1) = fit(x(1));
yn = fit(xn);
disp(sprintf('Ordenada inicial: %f', y(1)));
disp(sprintf('Ordenada final: %f', yn));
while 1
n = input('Ingrese el nmero de puntos (Salir < 2): ');
if n < 2 break; end;
xit = input(sprintf('Ingrese la abscisa a interpolar (salir < %f o > %f): ', x(1), xn));
% Validacin:
if (xit < x(1) | xit > xn) break; end;
% Clculo del resto de puntos:
if n > 2
h = (xn - x(1)) / (n - 1);
for i = 2 : n - 1
x(i) = x(i-1) + h;
y(i) = fit(x(i));
end
end;
x(n) = xn;
y(n) = yn;
% Clculo de la ordenada interpolada:
yit = InterpolacionLagrange(x, y, xit, n);
% Clculo de Y:
Y = fit(xit);
% Calcular los errores:
erAbs = abs(Y - yit);
erRel = erAbs / abs(Y);
% Desplegar resultados:
disp(sprintf('Ordenada (exacta): %f', Y));
disp(sprintf('Ordenada (interpolada con polinomios de Lagrange): %f', yit));
disp(sprintf('Error absoluto: %f', erAbs));
disp(sprintf('Error relativo: %f\n', erRel));
end;
% InterpolacionLagrange.m
function y = InterpolacionLagrange(vx, vy, x, n)
y = 0;
199
for j = 1 : n
L = 1;
for i = 1 : n
if i ~= j
L = L * (x - vx(i)) / (vx(j) - vx(i));
end;
end;
y = y + vy(j) * L;
end;
return
% fit.m
function y = fit(x)
y = 2 * x ^2 - 3;
return
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
200
Como se puede apreciar, existen diez cuadros de edicin, mismos que permitirn
ingresar los datos; y presentar los resultados; y dos botones, mismos que tienen la
funcionalidad:
- Procesar: Ejecuta el mtodo para interpolar mediante Polinomios de Lagrange y
desplegar los resultados. Y
- Salir: Finaliza la ejecucin del programa.
Se describe a continuacin el cdigo fuente:
Funcin Externa:
// Funcion a ser analizada:
static double fit(double x)
{
return 2 * Math.Pow(x, 2) - 3;
}
Evento que carga el formulario:
private void Form1_Load(object sender, EventArgs e)
{
CInterpolacion.AsignaMemoria(false, 50);
}
Botn Procesar:
201
202
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Salir:
private void cb_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
203
INTERPOLACION DE NEWTON:
Consiste en generar el nuevo punto, a partir de un grupo de puntos (x0; y0), (x1; y1),
(x2; y2), (xn; yn) para lo cual se emplea las frmulas de las diferencias divididas
definidas por:
Conocida una funcin real f(x), definida sobre las abscisas x0, x1, x2, ..., xn no
necesariamente equidistantes; entonces:
f ( x j , x j +1 ) =
f ( x j +1 ) f ( x j )
x j +1 x j
f ( x j , x j +1 , x j + 2 ) =
;4.6
f ( x j +1 , x j + 2 ) f ( x j , x j +1 )
x j+2 x j
;4.7
...
f ( x0 , x1 , x 2 ,...xi ,...x n ) =
j = 0,1,2,3,..., n 1
Las frmulas 4.6, 4.7 y 4.8 corresponden a las diferencias finitas divididas de primer,
segundo y n-esimo grado, respectivamente.
Para obtener el polinomio de interpolacin de Newton, se emplea la frmula:
P( x) = f ( x0 ) + f ( x0 , x1 )( x x0 ) + f ( x0 , x1 , x 2 )( x x0 )( x x1 ) + ...
+ f ( x0 , x1 , x 2 ,...xi ,...x n )( x x0 )( x x1 )...( x x n 1 );4.9
Ejemplo 4.6.
Obtener el polinomio de intepolacin del problema 4.1 aplicando las frmulas de
interpolacin de Newton mediante diferencias divididas, dados los puntos:
a) A(1, -1) y C(3, 15)
b) A(1, -1), B(2, 5) y C(3, 15)
Solucin:
a) El clculo de la primera diferencia dividida se desarrolla en la siguiente tabla:
2
Punto
y = f(x) = 2x 3
IDD
0
1
1
3
-1
15
204
De la tabla anterior: f(x0) = -1; f(x0, x1) = 8; con estos valores, se aplica 4.9:
P(x) = -1 + 8(x - 1) = 8x - 9
b) El clculo de las diferencias divididas se desarrolla en la siguiente tabla:
2
Punto
y = f(x) = 2x 3
IDD
IIDD
0
1
1
2
-1
5
6
10
15
De la tabla anterior: f(x0) = -1; f(x0, x1) = 6 y f(x0, x1, x2) = 2; con estos valores, se
aplica 4.9:
P(x) = -1 + 6(x - 1) + 2(x - 1) (x - 2) = -1 + 6x - 6 + 2x2 - 6x + 4 = 2x2 - 3.
Como se puede observar, el mtodo de Newton, ha producido los mismos resultados
que el de Lagrange.
Ejemplo 4.7
Calcular el polinomio de interpolacin de cuarto grado P(x, 4) y P(/2, 4) para y = f(x)
= excos(2x), en el dominio -2 x 2. Tomar cinco abscisas equidistantes.
Mediante interpolacin por diferencias divididas de Newton.
Solucin
El clculo de las diferencias divididas, se desarrolla en la siguiente tabla:
x =1.570796327
Punto
0
1
2
3
4
X
y = f(x)
IDD
IIDD
IIIDD
IVDD
-6.28318531 0.00186744 0.013160992
-3.14159265 0.04321392 0.304554469 0.04637671 0.10894819
0
1
7.047601352 1.07318924 2.52113665 0.191955859
3.14159265 23.1406926 163.0863767 24.8343424
6.28318531 535.491656
Por tanto f(x0) = 0.00186744; f(x0, x1) = 0.013160992; f(x0, x1, x2) = 0.04637671;
f(x0, x1, x2, x3) = 0.108948193; y f(x0, x1, x2, x3, x4) = 0.191955859.
Entonces el polinomio de interpolacin es:
P(x, 4) = 0.00187+0.01316(x + 2)+0.04638(x + 2)(x + )+0.10895(x + 2)(x + )x
+ 0.19196(x + 2)(x + )x(x - )
Sustituyendo valores y operando P(/2, 4) = -9.374026529.
Como se puede observar los resultados son bastante parecidos en los mtodos de
Lagrange y Newton. IDD, IIDD, corresponden a la primera y segunda diferencia
dividida, etc..
205
Ejemplo 4.8
Desarrollar una aplicacin Matlab y otra en Visual C#, que tenga una funcin para
calcular la ordenada y por el mtodo de interpolacin de Newton, mediante
diferencias divididas; dicha funcin debe recibir las coordenadas de dos puntos
(extremos) de la curva real; el valor de la abscisa x y el nmero de puntos; la
aplicacin debe calcular tambin los puntos intermedios y los errores absoluto y
relativo. Considerar los datos del ejemplo 4.1.
Solucin:
% PrbInterpolacionNewton.m
% Prueba interpolacin de Newton mediante diferencias divididas y calcula errores.
x(1) = input('Ingrese la abscisa inicial: ');
while 1
xn = input('Ingrese la abscisa final: ');
if xn > x(1) break; end;
end;
y(1) = fit(x(1));
yn = fit(xn);
disp(sprintf('Ordenada inicial: %f', y(1)));
disp(sprintf('Ordenada final: %f', yn));
while 1
n = input('Ingrese el nmero de puntos (Salir < 2): ');
if n < 2 break; end;
xit = input(sprintf('Ingrese la abscisa a interpolar (salir < %f o > %f): ', x(1), xn));
% Validacin:
if (xit < x(1) | xit > xn) break; end;
% Clculo del resto de puntos:
if n > 2
h = (xn - x(1)) / (n - 1);
for i = 2 : n - 1
x(i) = x(i-1) + h;
y(i) = fit(x(i));
end
end;
206
x(n) = xn;
y(n) = yn;
% Clculo de la ordenada interpolada:
yit = InterpolacionNewton(x, y, xit, n);
% Clculo de Y:
Y = fit(xit);
% Calcular los errores:
erAbs = abs(Y - yit);
erRel = erAbs / abs(Y);
% Desplegar resultados:
disp(sprintf('Ordenada (exacta): %f', Y));
disp(sprintf('Ordenada (interpolada mediante diferencias divididas de Newton): %f',
yit));
disp(sprintf('Error absoluto: %f', erAbs));
disp(sprintf('Error relativo: %f\n', erRel));
end;
% InterpolacionNewton.m
% Interpolacin de Newton mediante diferencias divididas.
function [y, vdd] = InterpolacionNewton(vx, vy, x, n)
vdd = vy;
y = vy(1);
p = 1;
for j = 1 : n - 1
k = 1;
p = p * (x - vx(j));
for i = 1 : n - j
vdd(i) = (vdd(i + 1) - vdd(i)) / (vx(k + j) - vx(k));
k = k + 1;
end;
y = y + vdd(1) * p;
end;
return
% fit.m
function y = fit(x)
y = 2 * x ^2 - 3;
return
207
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
208
Como se puede apreciar, existen diez cuadros de edicin, mismos que permitirn
ingresar los datos; y presentar los resultados; y dos botones, mismos que tienen la
funcionalidad:
- Procesar: Ejecuta el mtodo para interpolar mediante Diferencias Divididas de
Newton y desplegar los resultados. Y
- Salir: Finaliza la ejecucin del programa.
Se describe a continuacin el cdigo fuente:
Funcin Externa:
// Funcion a ser analizada:
static double fit(double x)
{
return 2 * Math.Pow(x, 2) - 3;
}
Evento que carga el formulario:
private void Form1_Load(object sender, EventArgs e)
{
CInterpolacion.AsignaMemoria(true, 50);
}
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
209
210
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Salir:
private void cb_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
211
212
2.- Ejecutar los programas en Matlab y Visual C# para interpolacin mediante los
polinomios de Lagrange, incorporando la funcin del numeral 1.
3.-Dadas las condiciones de la funcin del problema del numeral 1, calcular el
polinomio de interpolacin y P(/4, 4), as como los errores absoluto y relativo; a
travs de interpolacin por diferencias divididas de Newton. Analizar los resultados.
213
5.- Conocidos los puntos (-2; -35), (1; 1), (2; 29), (3; 125), interpolar para x = 1.5,
usando el mtodo de Aitken.
214
5
AJUSTE DE CURVAS
215
216
CAPITULO V
AJUSTE DE CURVAS
En este captulo, se analizarn los siguientes temas:
Generalidades.
Deduccin del mecanismo para ajustar una muestra de puntos a una funcin
polinmica a travs del mtodo de mnimos cuadrados.
Implementacin del mtodo anteriormente descrito en Matlab y Visual C#.NET.
Ajuste de curvas por interpolacin mediante Polinomios de Lagrange.
Implementacin del mtodo anteriormente descrito en Matlab y Visual C#.NET.
Ajuste de curvas por interpolacin mediante Diferencias Divididas de Newton.
Implementacin del mtodo anteriormente descrito en Matlab y Visual C#.NET.
Generalidades.
Introduccin:
En el presente capitulo, se presenta un modelo estadstico, mismo que permitir
ajustar una muestra obtenida de la medicin de un experimento del mundo real (en
medicina, biologa, informtica, fsica, qumica, administracin, mercadeo, economa,
entre otras disciplinas); a una funcin polinmica de la forma:
y = f(x, n) = A0 + A1x + A2x2 + + An-1xn-1 + Anxn
(5.1)
Donde: y es la variable dependiente; x es la variable independiente; n es el grado de
la funcin polinmica; y A0 , A1, A2,, An-1, An son los coeficientes de la funcin
polinmica; dicho modelo estadstico se denomina de los mnimos cuadrados y
previo a su descripcin, se realizan las siguientes definiciones:
Muestra estadstica:
(o muestra aleatoria o simplemente muestra) es un subconjunto de casos o
individuos de una poblacin estadstica. Se obtienen con la intencin de inferir
propiedades de la totalidad de la poblacin, para lo cual deben ser representativas de
la misma. Para cumplir esta caracterstica la inclusin de sujetos en la muestra debe
seguir una tcnica de muestreo. Por otra parte, en ocasiones, el muestreo puede ser
ms exacto que el estudio de toda la poblacin porque el manejo de un menor
217
1
85
2
93
3
61
4
73
5
110
6
93
7
89
8
54
9
91
10
60
11
79
12
87.5
218
S = (i y )
(5.5)
i =1
S = (i A0 - A1x - A2x
- - An - 1 x
n -1
- Anx )
(5.6)
i =1
(5.7)
pA0 + A1 xi + A2 xi + ... + An - 1 xi
i =1
i =1
i =1
219
n 1
+ An xi = i
i =1
i =1
i =1
i =1
n +1
A0 xi + A1 xi + A2 xi + ... + An - 1 xi + An xi
i =1
i =1
i =1
= xii
i =1
A0 xi + A1 xi + A2 xi + ... + An - 1 xi
i =1
i =1
i =1
i =1
220
n +1
+ An xi
i =1
n+2
= xi 2 i
i =1
n +1
A0 xi + A1 xi
i =1
n+ 2
+ A2 xi
i =1
2n 1
+ ... + An - 1 xi
i =1
i =1
2n
+ An xi
i =1
= x i n i
i =1
2
n
p
p
p
x
x
...
x
i
i
i
i =1
i =1
i =1
i = 1
n +1 A
3
2
p
p
p
p
p
0
xi xi
xi
... xi A
x
i i
i = 1
1
i =1
i =1
i =1
i =1
+
2
3
4
n
2
p
A 2 = p
p
p
p
(5.8)
2
xi
x
i i
x
x
...
x
i
i
i
i = 1
... i = 1
i =1
i =1
i =1
A
...
...
p
2n
p n p n +1 p n + 2
p
xi n i
xi
x
x
...
x
i
i
i
i = 1
i = 1
i =1
i =1
i =1
Como se puede observar, se genera un sistema de n + 1 ecuaciones por n+1
incgnitas, plantendolo matricialmente
221
12
12
xi
i = 1
12 2
xi
i = 1
12 3
xi
i = 1 4
12
xi
i = 1 5
12
xi
i = 1
12
x x
i
i =1
12
i =1
12
xi
i =1
12
i
5
i
6
xi
i =1
i =1
12
i =1
12
i
7
xi
i =1
12
i =1
12
7
i
i =1
12
xi
i =1
i =1
12
xi
x
i =1
12
i =1
12
i =1
i =1
12
xi
i =1
12
i =1
i = 1
6
12
12
xi A xi i
0
i =1
i = 1
7 A
12
1 12 2
xi
xi i
i =1
A 2 = i = 1
8
12
A 3 12 3
xi xi i
i =1
A 4 i =p 1
9
12
4
A 5
xi i
x
i
i12= 1
i =1
10
5
12
xi i
x
i
i = 1
i =1
12
i =1
12
xi
i =1
12
12
i =1
12
i =1
12
xi
x
i =1
12
i =1
12
12
i =1
12
i =1
12
12
975.5
x^2
1
4
9
16
25
36
49
64
81
100
121
144
x^3
1
8
27
64
125
216
343
512
729
1000
1331
1728
x^4
1
16
81
256
625
1296
2401
4096
6561
10000
14641
20736
x^5
1
32
243
1024
3125
7776
16807
32768
59049
100000
161051
248832
x^6
1
64
729
4096
15625
46656
117649
262144
531441
1000000
1771561
2985984
650
222
x^7
1
128
2187
16384
78125
279936
823543
2097152
4782969
10000000
19487171
35831808
x^8
1
256
6561
65536
390625
1679616
5764801
16777216
43046721
100000000
214358881
429981696
x^9
1
512
19683
262144
1953125
10077696
40353607
134217728
387420489
1000000000
2357947691
5159780352
x^10
1
1024
59049
1048576
9765625
60466176
282475249
1073741824
3486784401
1E+10
2.5937E+10
6.1917E+10
y*x
85
186
183
292
550
558
623
432
819
600
869
1050
y*x
85
372
549
1168
2750
3348
4361
3456
7371
6000
9559
12600
y*x
85
744
1647
4672
13750
20088
30527
27648
66339
60000
105149
151200
y*x
85
1488
4941
18688
68750
120528
213689
221184
597051
600000
1156639
1814400
y*x
85
2976
14823
74752
343750
723168
1495823
1769472
5373459
6000000
12723029
21772800
6247
51619
481849
4817443
50294137
650
143
6084
1859
60710
21164
630708
236093
6735950
2636348
60710
1.3347e+003
630708
6735950
26026
1.1583e+004
6735950
73399404
812071910
3.7828e+005 3.0116e+005 9.4135e+004
73399404
812071910 9092033028 1.02769E+11
4963530
5.3346e+006 3.0594e+006 7.0720e+005
975.5
81.2917
147.8863636
6247
-0.6556
-87.47589896
51619
-0.0014
35.36357851
481849
0.1020
-5.840924088
4817443
0.0724
0.41567899
50294137
-0.010562783
0.010562783
223
5
2
224
225
226
227
ajuste.dat:
AJUSTE DE POLINOMIOS POR EL METODO DE LOS MINIMOS CUADRADOS
VECTOR DE ABSCISAS:
1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000
9.000000 10.000000 11.000000 12.000000
VECTOR DE ORDENADAS:
85.000000 93.000000 61.000000 73.000000 110.000000 93.000000 89.000000
54.000000
91.000000 60.000000 79.000000 87.500000
COEFICIENTES DEL POLINOMIO DE GRADO N = 11 DE LA FORMA:
F(X)=A[0]+A[1]*X+...+A[N-1]*X**(N-1)+A[N]*X**N
DADOS EN EL ORDEN A[0],A[1],...,A[N-1],A[N]
2349.158228
-7404.022421
9968.630885
-7339.915015
3323.094358
984.274611 196.774301 -26.766037
2.441228 -0.142745 0.004830 -0.000072
COEFICIENTES DEL POLINOMIO DE GRADO N = 5 DE LA FORMA:
F(X)=A[0]+A[1]*X+...+A[N-1]*X**(N-1)+A[N]*X**N
DADOS EN EL ORDEN A[0],A[1],...,A[N-1],A[N]
228
229
Para que funcione el proceso de graficacin, debe activarse el evento paint, en este
caso debe accederse a las propiedades del pictureBox del formulario, realizar los
pasos que se detalla grficamente a continuacin:
:
Luego de hacer doble clic sobre el evento Paint, se genera el siguiente cdigo fuente:
private void pictureBox1_Paint(object sender, PaintEventArgs e)
{
}
Se describe a continuacin el cdigo fuente:
Variables y funciones miembro de la clase Form1:
// Datos Miembro:
Random r = new Random(DateTime.Now.Millisecond);
private Bitmap offScreenBitMap = null;
double EH, EV, MAXX, MAXY, MINX, MINY, X0, X1, Y0, Y1;
double MAXPX, MAXPY, CX, CY;
int n;
// Funciones miembro:
public Form1()
{
InitializeComponent();
}
230
private void MaxMin (double [] v, int n, ref double max, ref double min)
{
max = v[0];
min = v[0];
for (int i = 1; i < n; i++)
{
if (max < v[i]) max = v[i];
if (min > v[i]) min = v[i];
}
}
private void calculaGraficaPolinomio()
{
int i;
// Desplegar resultados del polinomio:
listBox1.Items.Add("COEFICIENTES DEL POLINOMIO DE GRADO N = " + n.ToString() + " DE LA
FORMA");
listBox1.Items.Add("F(X)=A[0]+A[1]*X+...+A[N-1]*X**(N-1)+A[N]*X**N");
listBox1.Items.Add("DADOS EN EL ORDEN A[0],A[1],...,A[N-1],A[N]");
// Calcular el polinomio de grado n:
CGenMatrizYVector.n = ++n;
CGenMatrizYVector.GeneraMatrizYVector();
CResolSisSimetricosEcuaciones.n = n;
CResolSisSimetricosEcuaciones.A = CGenMatrizYVector.T;
CResolSisSimetricosEcuaciones.B = CGenMatrizYVector.A;
CResolSisSimetricosEcuaciones.CroutSimetrico(true);
CLectEscrMatrices.n = n;
CLectEscrMatrices.b = CGenMatrizYVector.A;
CLectEscrMatrices.EscribeVector(" POLINOMIO:");
listBox1.Items.Add("");
/* DIBUJAR EL POLINOMIO */
// Invertir el orden del vector A, para calcular puntos:
double aux;
for (i = 0; i < n / 2; i++)
{
aux = CGenMatrizYVector.A[i];
CGenMatrizYVector.A[i] = CGenMatrizYVector.A[n - i - 1];
CGenMatrizYVector.A[n - i - 1] = aux;
}
CHorner.a = CGenMatrizYVector.A;
CHorner.n = n - 1;
X1 = X0;
int red = r.Next(255);
int green = r.Next(255);
int blue = r.Next(255);
Graphics g = Graphics.FromImage(offScreenBitMap);
Pen p = new Pen(Color.FromArgb(red, green, blue), 2);
for (i = 0; i < MAXPX; i++)
{
g.DrawEllipse(p, (float)(Math.Round((X1 - X0) * EH) + CX), (float)(2 * MAXPY - Math.Round(EV *
(CHorner.horner(X1) - Y0)) + CY), 2, 2);
X1 = X1 + 1 / EH;
}
pictureBox1.Invalidate();
p.Dispose();
231
g.Dispose();
Botn Graficar:
private void b_Graficar_Click(object sender, EventArgs e)
{
try
{
n = Int32.Parse(txtGra.Text);
if (n < 1 || n >= CGenMatrizYVector.p) {
MessageBox.Show("Error: el grado " + n.ToString() + " debe estar entre 1 y " +
(CGenMatrizYVector.p - 1).ToString());
return;
}
calculaGraficaPolinomio();
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
232
objes.AbrirArchivoEnt();
// Leer el tamao de la muestra:
CGenMatrizYVector.p = Int16.Parse(objes.LeerLinea());
if (CGenMatrizYVector.p > 0)
{
// Definir el arreglo de abscisas de la muestra:
CGenMatrizYVector.X = new double[CGenMatrizYVector.p];
CLectEscrMatrices.b = CGenMatrizYVector.X;
CLectEscrMatrices.n = CGenMatrizYVector.p;
CLectEscrMatrices.obent = objes;
CLectEscrMatrices.LeeVector();
// Definir el arreglo de ordenadas de la muestra:
CGenMatrizYVector.Y = new double[CGenMatrizYVector.p];
CLectEscrMatrices.b = CGenMatrizYVector.Y;
CLectEscrMatrices.LeeVector();
n = CGenMatrizYVector.p;
CGenMatrizYVector.A = new double[n];
CGenMatrizYVector.T = new double[n * (n + 1) / 2];
n--;
CLectEscrMatrices.lb = listBox1;
listBox1.Items.Add("AJUSTE DE POLINOMIOS POR EL METODO DE LOS MINIMOS
CUADRADOS");
CLectEscrMatrices.b = CGenMatrizYVector.X;
CLectEscrMatrices.EscribeVector(" DE ABSCISAS:");
listBox1.Items.Add("");
CLectEscrMatrices.b = CGenMatrizYVector.Y;
CLectEscrMatrices.EscribeVector(" DE ORDENADAS:");
listBox1.Items.Add("");
// Calcular las constantes de graficacin:
int i;
if (offScreenBitMap == null)
{
offScreenBitMap = new Bitmap(pictureBox1.Width, pictureBox1.Height);
}
Graphics g = Graphics.FromImage(offScreenBitMap);
Pen p = new Pen(Color.FromArgb(0, 0, 0), 2);
MAXPX = Math.Round(0.75 * pictureBox1.Width, 0);
MAXPY = Math.Round(0.40 * pictureBox1.Height, 0);
CX = Math.Round(0.11 * pictureBox1.Width, 0);
CY = Math.Round(0.08 * pictureBox1.Height, 0);
/* CALCULA EL MAXIMO Y MINIMO DE LAS COORDENADAS DE LOS PUNTOS */
MaxMin(CGenMatrizYVector.X, CGenMatrizYVector.p, ref MAXX, ref MINX);
MaxMin(CGenMatrizYVector.Y, CGenMatrizYVector.p, ref MAXY, ref MINY);
/* CALCULA ESCALAS Y PUNTOS INICIALES DE REFERENCIA */
EH = MAXPX / (MAXX - MINX);
EV = 2 * MAXPY / (MAXY - MINY);
X0 = 0;
Y0 = 0;
X1 = MAXX;
Y1 = MAXY;
if (MAXX < 0)
{
EH = MAXPX / Math.Abs(MINX);
X0 = MINX;
233
X1 = 0;
}
else if (MINX >= 0)
EH = MAXPX / MAXX;
else
{
X0 = MINX;
X1 = MAXX;
}
if (MAXY < 0)
{
EV = 2 * MAXPY / Math.Abs(MINY);
Y0 = MINY;
Y1 = 0;
}
else if (MINY >= 0)
EV = 2 * MAXPY / MAXY;
else
{
Y0 = MINY;
Y1 = MAXY;
}
/* DIBUJAR LOS EJES DEL SISTEMA DE REFERENCIA */
g.DrawLine(p, (float)CX, (float)(Math.Round(2 * MAXPY + Y0 * EV) + CY),
(float)(Math.Round((X1 - X0) * EH) + Math.Round(1.25 * CX)),
(float)(Math.Round(2
* MAXPY + Y0 * EV) + CY));
Font df = new Font("Arial", 14);
SolidBrush db = new SolidBrush(Color.Black);
g.DrawString("X", df, db, (float)(Math.Round((X1 - X0) * EH) + Math.Round(1.4375 * CX)),
(float)(Math.Round(2 * MAXPY + Y0 * EV) + CY - 3));
g.DrawLine(p, (float)(Math.Round(-X0 * EH) + CX), (float)(Math.Round(2 * MAXPY) + CY),
(float)(Math.Round(-X0 * EH) + CX),
(float)(Math.Round(2 * MAXPY - (Y1 - Y0) * EV) + Math.Round(0.5 * CY)));
g.DrawString("Y", df, db, (float)(Math.Round(-X0 * EH) + Math.Round(0.9375 * CX)),
(float)(Math.Round(2 * MAXPY - (Y1 - Y0) * EV)));
/* DIBUJAR LOS PUNTOS */
p.Dispose();
int red = r.Next(255);
int green = r.Next(255);
int blue = r.Next(255);
p = new Pen(Color.FromArgb(red, green, blue), 2);
for (i = 0; i < CGenMatrizYVector.p; i++)
g.DrawEllipse(p, (float)(Math.Round(EH * (CGenMatrizYVector.X[i] - X0)) + CX),
(float)(Math.Round(2 * MAXPY - EV * (CGenMatrizYVector.Y[i] - Y0)) + CY), 6, 6);
p.Dispose();
g.Dispose();
calculaGraficaPolinomio();
b_Graficar.Enabled = true;
}
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
234
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
Las
referencias
a
las
clases
CGenMatrizYVector,
CEntSalArchivos,
CResolSisSimetricosEcuaciones, CLectEscrMatrices, CHorner y CListBox, provienen
de la definicin de las respectivas clases incorporadas en la librera del archivo
CLibComunClases.cs, el cual se describe en el Apndice B de este libro.
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
235
236
237
Ajuste.txt:
AJUSTE DE POLINOMIOS POR EL METODO DE LOS MINIMOS CUADRADOS
VECTOR DE ABSCISAS:
1 2 3 4 5 6 7 8 9 10 11 12
VECTOR DE ORDENADAS:
85 93 61 73 110 93 89 54 91 60 79 87.5
COEFICIENTES DEL POLINOMIO DE GRADO N = 11 DE LA FORMA
F(X)=A[0]+A[1]*X+...+A[N-1]*X**(N-1)+A[N]*X**N
DADOS EN EL ORDEN A[0],A[1],...,A[N-1],A[N]
VECTOR POLINOMIO:
772.193092935719 -2798.38646674622 4447.78536029696 -3685.18347312394
1818.62937746471 -575.887139274014 121.605801964184 -17.3267057107536
1.64470659239433
-0.0995500372432104
0.00347066808153631
5.29594954157183E-05
COEFICIENTES DEL POLINOMIO DE GRADO N = 5 DE LA FORMA
F(X)=A[0]+A[1]*X+...+A[N-1]*X**(N-1)+A[N]*X**N
DADOS EN EL ORDEN A[0],A[1],...,A[N-1],A[N]
VECTOR POLINOMIO:
147.886363638063 -87.4758989808835 35.3635785181271 -5.8409240883619
0.415678990143525 -0.0105627828059102
238
239
240
end;
% Grafica el polinomio de interpolacin:
plot(xc, yc, 'm:');
pos = pos + 1;
if fdis ~= 1
% Generacin de los puntos de la curva real:
xr(1) = vx(1);
yr(1) = fit(xr(1));
delta = (vx(n) - vx(1)) / ni;
for i = 2 : ni + 1
xr(i) = xr(i - 1) + delta;
yr(i) = fit(xr(i));
end;
% Grafica la funcin real:
plot(xr, yr, colr(mod(pos, 6) + 1));
pos = pos + 1;
end;
xlabel('Ajuste de la curva mediante interpolacin por Polinomios de Lagrange');
hold off;
end;
% fit.m
function y = fit(x)
y = exp(x) * cos(2 * x);
return
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa.,
241
242
243
double EH, EV, MAXX, MAXY, MINX, MINY, X0, X1, Y0, Y1;
double MAXPX, MAXPY, CX, CY;
int n;
// Funciones miembro:
public Form1()
{
InitializeComponent();
}
private void MaxMin (double [] v, int n, ref double max, ref double min)
{
max = v[0];
min = v[0];
for (int i = 1; i < n; i++)
{
if (max < v[i]) max = v[i];
if (min > v[i]) min = v[i];
}
}
// Funcion a ser analizada:
static double fit(double x)
{
return Math.Exp(x) * Math.Cos(2 * x);
}
private void calculaGraficaPolinomio(bool real)
{
int i;
X1 = X0;
int red = r.Next(255);
int green = r.Next(255);
int blue = r.Next(255);
Graphics g = Graphics.FromImage(offScreenBitMap);
Pen p = new Pen(Color.FromArgb(red, green, blue), 2);
if (real)
for (i = 0; i < MAXPX; i++)
{
g.DrawEllipse(p, (float)(Math.Round((X1 - X0) * EH) + CX), (float)(2 * MAXPY - Math.Round(EV *
(fit(X1) - Y0)) + CY), 3, 3);
X1 = X1 + 1 / EH;
}
else
for (i = 0; i < MAXPX; i++)
{
g.DrawEllipse(p, (float)(Math.Round((X1 - X0) * EH) + CX), (float)(2 * MAXPY - Math.Round(EV
* (CInterpolacion.InterpolacionLagrange(X1, n) - Y0)) + CY), 1, 1);
X1 = X1 + 1 / EH;
}
pictureBox1.Invalidate();
p.Dispose();
g.Dispose();
}
244
Opcin Continua:
private void rbContinua_CheckedChanged(object sender, EventArgs e)
{
label1.Show();
txtNPuntos.Show();
txtLIni.Show();
txtLFin.Show();
b_Graficar.Enabled = false;
try
{
245
int np = Int32.Parse(txtNPuntos.Text);
if (np > 1)
b_Graficar.Enabled = true;
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Opcin Discreta:
private void rbDiscreta_CheckedChanged(object sender, EventArgs e)
{
label1.Hide();
txtNPuntos.Hide();
txtLIni.Hide();
txtLFin.Hide();
b_Graficar.Enabled = true;
}
Botn Graficar:
private void b_Graficar_Click(object sender, EventArgs e)
{
offScreenBitMap = null;
try
{
if (rbDiscreta.Checked) {
n = vdx.Length;
for (int i = 0; i < n; i++)
{
CInterpolacion.vx[i] = vdx[i];
CInterpolacion.vy[i] = vdy[i];
}
}
else
{
double h, li, lf;
n = Int32.Parse(txtNPuntos.Text);
li = double.Parse(txtLIni.Text);
lf = double.Parse(txtLFin.Text);
if (li >= lf) return;
h = (lf - li) / (n - 1);
CInterpolacion.vx[0] = li;
CInterpolacion.vy[0] = fit(CInterpolacion.vx[0]);
for (int i = 1; i < n; i++)
{
CInterpolacion.vx[i] = CInterpolacion.vx[i - 1] + h;
CInterpolacion.vy[i] = fit(CInterpolacion.vx[i]);
}
}
// Calcular las constantes de graficacin:
if (offScreenBitMap == null)
246
247
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
248
249
250
251
% fit.m
function y = fit(x)
y = exp(x) * cos(2 * x);
return
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa.,
252
253
254
255
b_Graficar.Enabled = false;
try
{
int np = Int32.Parse(txtNPuntos.Text);
if (np > 1)
b_Graficar.Enabled = true;
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Opcin Continua:
private void rbContinua_CheckedChanged(object sender, EventArgs e)
{
label1.Show();
txtNPuntos.Show();
txtLIni.Show();
txtLFin.Show();
b_Graficar.Enabled = false;
try
{
int np = Int32.Parse(txtNPuntos.Text);
if (np > 1)
b_Graficar.Enabled = true;
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
Opcin Discreta:
private void rbDiscreta_CheckedChanged(object sender, EventArgs e)
{
label1.Hide();
txtNPuntos.Hide();
txtLIni.Hide();
txtLFin.Hide();
b_Graficar.Enabled = true;
}
Botn Graficar:
private void b_Graficar_Click(object sender, EventArgs e)
{
offScreenBitMap = null;
try
{
if (rbDiscreta.Checked) {
n = vdx.Length;
256
}
else
{
double h, li, lf;
n = Int32.Parse(txtNPuntos.Text);
li = double.Parse(txtLIni.Text);
lf = double.Parse(txtLFin.Text);
if (li >= lf) return;
h = (lf - li) / (n - 1);
CInterpolacion.vx[0] = li;
CInterpolacion.vy[0] = fit(CInterpolacion.vx[0]);
for (int i = 1; i < n; i++)
{
CInterpolacion.vx[i] = CInterpolacion.vx[i - 1] + h;
CInterpolacion.vy[i] = fit(CInterpolacion.vx[i]);
}
}
// Calcular las constantes de graficacin:
if (offScreenBitMap == null)
{
offScreenBitMap = new Bitmap(pictureBox1.Width, pictureBox1.Height);
}
Graphics g = Graphics.FromImage(offScreenBitMap);
Pen p = new Pen(Color.FromArgb(0, 0, 0), 2);
MAXPX = Math.Round(0.75 * pictureBox1.Width, 0);
MAXPY = Math.Round(0.40 * pictureBox1.Height, 0);
CX = Math.Round(0.11 * pictureBox1.Width, 0);
CY = Math.Round(0.08 * pictureBox1.Height, 0);
/* CALCULA EL MAXIMO Y MINIMO DE LAS COORDENADAS DE LOS PUNTOS */
MaxMin(CInterpolacion.vx, n, ref MAXX, ref MINX);
MaxMin(CInterpolacion.vy, n, ref MAXY, ref MINY);
/* CALCULA ESCALAS Y PUNTOS INICIALES DE REFERENCIA */
EH = MAXPX / (MAXX - MINX);
EV = 2 * MAXPY / (MAXY - MINY);
X0 = 0;
Y0 = 0;
X1 = MAXX;
Y1 = MAXY;
if (MAXX < 0)
{
EH = MAXPX / Math.Abs(MINX);
X0 = MINX;
X1 = 0;
}
else if (MINX >= 0)
EH = MAXPX / MAXX;
else
{
X0 = MINX;
X1 = MAXX;
257
}
if (MAXY < 0)
{
EV = 2 * MAXPY / Math.Abs(MINY);
Y0 = MINY;
Y1 = 0;
}
else if (MINY >= 0)
EV = 2 * MAXPY / MAXY;
else
{
Y0 = MINY;
Y1 = MAXY;
}
/* DIBUJAR LOS EJES DEL SISTEMA DE REFERENCIA */
g.DrawLine(p,(float)CX,(float)(Math.Round(2*MAXPY+Y0*EV) + CY),
(float)(Math.Round((X1 - X0) * EH) + Math.Round(1.25 * CX)),
(float)(Math.Round(2 * MAXPY + Y0 * EV) + CY));
Font df = new Font("Arial", 14);
SolidBrush db = new SolidBrush(Color.Black);
g.DrawString("X", df, db, (float)(Math.Round((X1 - X0) * EH) +
Math.Round(1.4375 * CX)),
(float)(Math.Round(2 * MAXPY + Y0 * EV) + CY - 3));
g.DrawLine(p, (float)(Math.Round(-X0 * EH) + CX),
(float)(Math.Round(2 * MAXPY) + CY),
(float)(Math.Round(-X0 * EH) + CX),
(float)(Math.Round(2*MAXPY-(Y1-Y0) * EV) + Math.Round(0.5 * CY)));
g.DrawString("Y", df, db, (float)(Math.Round(-X0 * EH) +
Math.Round(0.9375 * CX)),
(float)(Math.Round(2 * MAXPY - (Y1 - Y0) * EV)));
/* DIBUJAR LOS PUNTOS */
p.Dispose();
int red = r.Next(255);
int green = r.Next(255);
int blue = r.Next(255);
p = new Pen(Color.FromArgb(red, green, blue), 2);
for (int i = 0; i < n; i++)
g.DrawEllipse(p, (float)(Math.Round(EH * (CInterpolacion.vx[i] X0)) + CX),
(float)(Math.Round(2 * MAXPY - EV * (CInterpolacion.vy[i] - Y0))
+ CY), 8, 8);
p.Dispose();
g.Dispose();
calculaGraficaPolinomio(false);
if (rbContinua.Checked)
calculaGraficaPolinomio(true);
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Salir:
258
259
260
261
262
6
DIFERENCIACIN (DERIVACIN) NMERICA
263
264
CAPITULO VI
Generalidades
Este captulo se centra en el clculo de una aproximacin de la derivada de una
funcin contnua o discreta en un punto dado.
Se revisarn los siguientes temas:
f ' ( x) = lim
h 0
f ( x + h) f ( x ) f ( x + h) f ( x )
; (6.1)
h
h
265
Solucin:
El proceso se describe en la tabla siguiente:
266
cual debe calcularse la derivada aplicando la frmula 6.1; calcular adems el error
absoluto y el error relativo Asumir un valor de h = 1e-8.
Solucin:
Programa en Matlab:
% DerivadaIOrden.m
clc; % Limpia pantalla
clear; % Libera espacio de variables
syms x % Define como simbolo a x
f = 4/(2 - exp(x^2)) % Define la funcin f(x)
derivada = diff(f) % Genera simblicamente f'(x)
h = 1e-8;
disp('Clculo de la derivada de primer orden en un punto dado');
while 1
x = input('Ingresar x (salir < -5 o > 5): ');
if (abs(x) > 5) break; end;
fpe = eval(derivada);
fx = eval(f);
x1 = x;
x = x + h;
fxh = eval(f);
fp = (fxh - fx) / h;
era = abs(fpe - fp);
err = era / abs(fpe) * 100;
disp(sprintf('fp(%f) = %f, fpexacta(%f) = %f, Error absoluto = %f, Error relativo =
%f%%\n', x1, fp, x1, fpe, era, err));
end;
Los resultados, se presentan a continuacin:
267
268
269
Botn Salir:
270
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
double x, fx, fxh, fp;
try
{
x = double.Parse(txtX.Text);
txtFp.Text = "";
txtFx.Text = "";
txtFxh.Text = "";
if (x < CInterpolacion.vx[0] || x > CInterpolacion.vx[5]) return;
fx = CInterpolacion.InterpolacionNewton(x, n);
fxh = CInterpolacion.InterpolacionNewton(x + h, n);
fp = (fxh - fx) / h;
txtFp.Text = Math.Round(fp, 12).ToString();
txtFx.Text = Math.Round(fx, 12).ToString();
txtFxh.Text = Math.Round(fxh, 12).ToString();
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
271
h2
h3
h 4 ( 4)
f ' ' ( x) +
f ' ' ' ( x) +
f ( x) + ...(6.2)
2!
3!
4!
Sustituyendo h por -h en 6.2, se obtiene:
h2
h3
h 4 ( 4)
f ( x h) = f ( x) hf ' ( x) +
f ' ' ( x)
f ' ' ' ( x) +
f ( x) ...(6.3)
2!
3!
4!
f ( x + h) = f ( x) + hf ' ( x) +
272
4h 2
8h 3
16h 4 ( 4)
f ' ' ( x) +
f ' ' ' ( x) +
f ( x) + ...(6.4)
2!
3!
4!
4h 2
8h 3
16h 4 ( 4)
f ' ' ( x)
f ' ' ' ( x) +
f ( x) ...(6.5)
2!
3!
4!
h2
h 4 ( 4)
f ' ' ( x) +
f ( x) + ...); (6.6)
2!
4!
Realizando 6.2 - 6.3, se obtiene:
f ( x + h) + f ( x h) = 2( f ( x) +
h3
h 5 ( 5)
f ' ' ' ( x) +
f ( x) + ...); (6.7)
3!
5!
Realizando 6.4 + 6.5, se obtiene:
4h 2
16h 4 ( 4)
f ( x + 2h) + f ( x 2h) = 2( f ( x) +
f ' ' ( x) +
f ( x) + ...); (6.8)
2!
4!
f ( x + h) f ( x h) = 2(hf ' ( x) +
8h 3
32h 5 (5)
f ' ' ' ( x) +
f ( x) + ...), (6.9)
3!
5!
f ( x)
f ( x) ...; (6.12)
12
360
h2
Para un valor pequeo de h, se puede descartar desde el tercer trmino de 6.12 y
sustituyendo el segundo por O(h2):
f ( x + h ) 2 f ( x ) + f ( x h)
f ' ' ( x) =
+ O(h 2 ); (6.13)
h2
Las derivadas superiores a
273
segundo orden, se las puede obtener, efectundo operaciones entre las frmulas de
6.2 a 6.9 resultando ser:
f ( x + 2 h) 2 f ( x + h) + 2 f ( x h) f ( x 2 h)
+ O(h 2 ); (6.14)
3
2h
f ( 4 ) ( x) =
f ( x + 2h ) 4 f ( x + h ) + 6 f ( x ) + 4 f ( x h ) f ( x 2 h )
+ O(h 2 ); (6.15)
h4
Ejemplo 6.7
Implementar un programa en Matlab para resolver el problema del ejemplo 6.6,
generarando las derivadas simblicamente, ingresar el valor en el que se desean
evaluar las derivadas de orden superior, las cuales deben estar en el dominio -5 x
5, luego de lo cual deben calcularse las derivadas aplicando las frmulas 6.13, 6.14 y
6.15; calcular adems los errores absoluto y relativo.
Solucin:
Programa en Matlab:
% DerivadaOrdenSup.m
% Calculo de derivadas de orden superior usando
% aproximacin de diferencias centradas.
clc; % Limpia pantalla
clear; % Libera espacio de variables
syms x % Define como simbolo a x
f = exp(x)*cos(2*x) % Define la funcin f(x)
274
275
f ' ( x) =
f ( x + 2h) + 4 f ( x + h) 3 f ( x) h 2
h 3 ( 4)
+
f ' ' ' ' ( x) +
f ( x) + ...; (6.16)
2h
3
4
f ' ( x) =
f ( x + 2h) + 4 f ( x + h) 3 f ( x )
+ O(h 2 ); (6.17)
2h
276
Las derivadas superiores a f(x), se las obtiene al efectuar operaciones entre f(x) , f(x
+ h) ,
f(x + 2h) , f(x + 3h) , f(x + 4h). etc. y se las describe a continuacin:
f ( x + 3h) + 4 f ( x + 2h) 5 f ( x + h) + 2 f ( x)
+ O(h 2 ); (6.18)
2
h
3 f ( x + 4h) + 14 f ( x + 3h) 24 f ( x + 2h) + 18 f ( x + h) 5 f ( x)
f ' ' ' ( x) =
+ O(h 2 ); (6.19)
3
2h
2 f ( x + 5h) + 11 f ( x + 4h) 24 f ( x + 3h) + 26 f ( x + 2h) 14 f ( x + h) + 3 f ( x)
f ( 4 ) ( x) =
+ O(h 2 ); (6.20)
h4
f ' ' ( x) =
Ejemplo 6.8
Resolver el problema del ejemplo 6.6 usando aproximacin de diferencias
adelantadas.
Solucin:
El proceso se encuentra desarrollado en la siguiente tabla:
Ejemplo 6.9
Implementar un programa en Matlab para resolver el problema del ejemplo 6.8,
generarando las derivadas simblicamente, ingresar el valor en el que se desean
evaluar las derivadas de orden superior, las cuales deben estar en el dominio -5 x
5, luego de lo cual deben calcularse las derivadas aplicando las frmulas 6.18, 6.19 y
6.20; calcular adems los errores absoluto y relativo.
Solucin:
Programa en Matlab:
% DerivadaOrdenSupDifAdelan.m
% Calculo de derivadas de orden superior usando
% aproximacin de diferencias adelantadas.
clc; % Limpia pantalla
clear; % Libera espacio de variables
syms x % Define como simbolo a x
f = exp(x)*cos(2*x) % Define la funcin f(x)
derivada2 = diff(f,2) % Genera simblicamente f''(x)
derivada3 = diff(f,3) % Genera simblicamente f'''(x)
277
278
f ( x 2h) 4 f ( x h) + 3 f ( x) h 2
h 3 ( 4)
f ' ( x) =
f ( x 2h ) 4 f ( x h) + 3 f ( x )
+ O(h 2 ); (6.22)
2h
Las derivadas superiores a f(x), se las obtiene al efectuar operaciones entre f(x) , f(x
- h) ,
f(x - 2h) , f(x - 3h) , f(x - 4h). etc. y se las describe a continuacin:
f ' ( x) =
f ( x 3h) + 4 f ( x 2h) 5 f ( x h) + 2 f ( x)
+ O(h 2 ); (6.23)
2
h
3 f ( x 4h) 14 f ( x 3h) + 24 f ( x 2h) 18 f ( x h) + 5 f ( x)
f ' ' ' ( x) =
+ O(h 2 ); (6.24)
3
2h
2 f ( x 5h) + 11 f ( x 4h) 24 f ( x 3h) + 26 f ( x 2h) 14 f ( x h) + 3 f ( x)
f ( 4 ) ( x) =
+ O(h 2 ); (6.25)
h4
f ' ' ( x) =
Las diferencias adelantadas y retrasadas permiten calcular f'(x), f''(x), f'''(x), etc.
teniendo como entradas tres, cuatro, cinco, etc. puntos, respectivamente.
Ejemplo 6.10
Resolver el problema del ejemplo 6.6 usando aproximacin de diferencias
retrasadas.
Solucin:
El proceso se encuentra desarrollado en la siguiente tabla:
279
Ejemplo 6.11
Implementar un programa en Matlab para resolver el problema del ejemplo 6.10,
generarando las derivadas simblicamente, ingresar el valor en el que se desean
evaluar las derivadas de orden superior, las cuales deben estar en el dominio -5 x
5, luego de lo cual deben calcularse las derivadas aplicando las frmulas 6.23, 6.24 y
6.25; calcular adems los errores absoluto y relativo.
Solucin:
Programa en Matlab:
% DerivadaOrdenSupDifRetras.m
% Calculo de derivadas de orden superior usando
% aproximacin de diferencias retrasadas.
clc; % Limpia pantalla
clear; % Libera espacio de variables
syms x % Define como simbolo a x
f = exp(x)*cos(2*x) % Define la funcin f(x)
derivada2 = diff(f,2) % Genera simblicamente f''(x)
derivada3 = diff(f,3) % Genera simblicamente f'''(x)
derivada4 = diff(f,4) % Genera simblicamente f(4)(x)
h = 0.1;
disp('Clculo de la derivada de orden superior en un punto dado');
disp(' usando aproximacin de diferencias retrasadas.');
while 1
x = input('Ingresar x (salir < -5 o > 5): ');
if (abs(x) > 5) break; end;
fpe2 = eval(derivada2);
fpe3 = eval(derivada3);
fpe4 = eval(derivada4);
fx = eval(f);
x1 = x;
x = x - h;
fxmh = eval(f);
x = x - h;
fxm2h = eval(f);
x = x - h;
fxm3h = eval(f);
280
x = x - h;
fxm4h = eval(f);
x = x - h;
fxm5h = eval(f);
fp2 = (-fxm3h + 4*fxm2h - 5*fxmh + 2*fx) / (h^2);
fp3 = (3*fxm4h - 14*fxm3h + 24*fxm2h - 18*fxmh + 5*fx) / (2*h^3);
fp4 = (-2*fxm5h + 11*fxm4h - 24*fxm3h + 26*fxm2h - 14*fxmh + 3*fx) / (h^4);
era2 = abs(fpe2 - fp2);
err2 = era2 / abs(fpe2) * 100;
era3 = abs(fpe3 - fp3);
err3 = era3 / abs(fpe3) * 100;
era4 = abs(fpe4 - fp4);
err4 = era4 / abs(fpe4) * 100;
disp(sprintf('fp2(%f) = %f, fp2exacta(%f) = %f, Error absoluto = %f, Error relativo =
%f%%', x1, fp2, x1, fpe2, era2, err2));
disp(sprintf('fp3(%f) = %f, fp3exacta(%f) = %f, Error absoluto = %f, Error relativo =
%f%%', x1, fp3, x1, fpe3, era3, err3));
disp(sprintf('fp4(%f) = %f, fp4exacta(%f) = %f, Error absoluto = %f, Error relativo =
%f%%', x1, fp4, x1, fpe4, era4, err4));
end;
Los resultados, se presentan a continuacin:
281
Ejemplo 6.12
Dado el polinomio f(x) = 5x4 3x2 + 2x -20; calcular f(3) y f(4)(3).
Solucin:
Aplicando sucesivamente dos veces la regla de derivacin y sustituyendo x = 3:
f(x) = 20x3 6x + 2; f(x) = 60x2 6; para x = 3:
f(3) = 60(3)2 6 = 534.
Aplicando sobre f(x) sucesivamente dos veces la regla de derivacin y sustituyendo
x = 3:
f(x) = 120x; f(4)(x) = 120; f(4)(3) = 120.
Ejemplo 6.13
Implementar en Matlab y Visual C# un programa para ingresar el grado y los
coeficientes del polinomio, el orden de la derivada y el valor a evaluar y generar
simblicamente el polinomio ingresado y la derivada resultante, calcular adems el
valor de la derivada.
Solucin:
Programa en Matlab:
% PDerivNumPol.m
% Diferenciacin numrica de polinomios.
% FUNCION PRINCIPAL:
clc; % Limpia pantalla
clear; % Libera espacio de variables
global A N
disp('Derivacin de polinomios: ');
[N, A] = LeePolinomio;
disp('Polinomio ingresado: ');
EscribePolinomio(A, N + 1);
while 1
od = input('Ingrese el orden de la derivada (salir < 1): ');
if (od < 1)
break
end;
x = input('Ingrese el valor a evaluar: ');
fprintf(sprintf('Derivada %d(%f) = %f\n', od, x, DerivNumPol(od, x)));
end;
% EscribePolinomio.m
function EscribePolinomio(A, m)
gr = m - 1;
282
ep = [];
for i = 1 : m
if (A(i) ~= 0)
if (i == 1)
if (gr > 0)
ep = [ep, sprintf('%5.2f*', A(i))];
else
ep = [ep, sprintf('%5.2f', A(i))];
end;
else
if (A(i) < 0)
ep = [ep, ' -'];
else
ep = [ep, ' +'];
end;
if (abs(A(i)) ~= 1)
if (gr > 0)
ep = [ep, sprintf('%5.2f*', abs(A(i)))];
else
ep = [ep, sprintf('%5.2f', abs(A(i)))];
end;
end;
end;
if (gr > 0)
if (gr == 1)
ep = [ep, 'x'];
else
ep = [ep, sprintf('x^%d', gr)];
end;
end;
end;
gr = gr - 1;
end;
disp(ep);
return
% DerivNumPol.m
% Diferenciacin numrica de polinomios.
function d = DerivNumPol(nd, x)
global A N;
if (nd > N)
d = 0;
return;
end;
% Obtener una copia de A:
for i = 1 : N + 1
283
B(i) = A(i);
end;
% Aplicar la regla de derivacin:
gd = N;
for i = 1 : nd
gd = gd - 1;
for j = 1 : gd + 1
B(j) = (gd - j + 2) * B(j);
end;
end;
disp('Polinomio resultante: ');
EscribePolinomio(B, gd + 1);
% Evaluar para x al polinomio resultante:
d = EvaluaPolinomio(B, gd, x);
return
% EvaluaPolinomio.m
function b = EvaluaPolinomio(A, N, x)
b = A(1);
for i = 1 : N
b = A(i + 1) + b * x;
end
return
La definicin de la
funcin
LeePolinomio,
se defini en un captulo
anterior. Los resultados
de la ejecucin de este
programa, se presenta
en la siguiente pantalla:
284
285
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
286
}
res += "(" + x.ToString() + ") = " + rd.ToString();
listBox1.Items.Add(res);
listBox1.Items.Add("");
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
287
288
Conforme a los resultados obtenidos, se puede concluir que para la f(x) indicada, se
comete menos errores en f(2) y f(2) para la aproximacin de diferencias centrada y
en f(4)(2), los errores son menores para la aproximacin de diferencias adelantada.
Para valores de h ms pequeos (por ejemplo, h = 0.01), se obtienen resultados de
la siguiente tabla:
Como se puede observar, son resultados, bastante interesantes, en todo caso, para
calcular la derivada mediante aproximacin de diferencias, se debe realizar un
proceso iterativo, haciendo variar h, hasta lograr la convergencia.
289
2.- Dada la siguiente muestra de ventas (en miles de USD) de un producto en los
doce meses del ao:
Mes
Venta
1
85
2
93
3
61
4
73
5
110
6
93
7
89
8
54
9
91
10
60
11
79
12
87.5
290
3.- Obtener f''(3), f'''(3) y f(4)(3) y sus correspondientes errores absoluto y relativo
mediante aproximacin de diferencias centradas, para h = 0.01 a partit de la funcin
del problema 1.
4.- Obtener f''(3), f'''(3) y f(4)(3) y sus correspondientes errores absoluto y relativo
mediante aproximacin de diferencias adelantadas, para h = 0.01 a partit de la
funcin del problema 1.
5.- Obtener f''(3), f'''(3) y f(4)(3) y sus correspondientes errores absoluto y relativo
mediante aproximacin de diferencias retrasadas, para h = 0.01 a partit de la funcin
del problema 1.
6.- Implementar los programas en Matlab y Visual C# para calcular f''(x), f'''(x) y f(4)(x)
mediante aproximacin de diferencias centradas, mediante un proceso iterativo,
haciendo variar a h = 10-K para todo 1 k 10, k es un nmero entero; tomar como
entradas el argumento x, y la funcin del problema 1. Especificar para que valor de h,
se tiene la mejor aproximacin.
291
292
7
INTEGRACIN NUMRICA
293
294
CAPITULO VII
INTEGRACION NUMERICA.
Generalidades
La integracin numrica consiste en determinar el valor aproximado de la integral
b
definida:
F ( x) = f ( x)dx; (7.1) .donde f(x) es una funcin integrable en el dominio a x b,
a
F(x) es la primitiva o antiderivada; aplicando un mtodo
numrico. Tiene la ventaja de calcular integrales que implican un proceso
analticamente complejo o imposible de resolver.
Definicin: Dada la funcin continua
f(x) definida en el dominio a x b,
entonces la frmula 7.1 se define
como el area bajo (sobre) la curva
f(x),
comprendida
entre
las
ordenadas que pasan por las
abscisas a, b y el eje X como lo
muestra la figura 7.1.
Para poder calcular el area de esa
regin se la divide en un grupo de
rectngulos de base hi y altura yi.
Para conseguir una mayor precisin,
la base debe ser lo suficientemente pequea, de manera que:
hi se le denomina tambin Longitud de
b
n
n
f ( x)dx = lim h y = lim h y ; (7.2) paso; n es el nmero de areas pequeas.
hi 0
i =1
i =1
295
de
Solucin:
1
0 i =1
n
n
( x 2 ) 2i 1
(1) i +1
(1) i +1
4i 2
=
dx =
x
dx
(2i 1)!
i =1 ( 2i 1)! 0
i =1 ( 2i 1)!( 4i 1)
1 1
1
1
1
1
= 0.3102683
3 42 1320 75600 6894720 918086400
Ejemplo 7.2
Implementar un programa en Matlab, que tenga como entradas la funcin f(x) a
integrarse, el intervalo y el nmero de trminos de la serie de Taylor y que calcule el
valor de la integral, adems el error absoluto y relativo entre el valor que produce
Matlab y el calculado. Desarrollar para f(x) = sen(x2), repetir el proceso ingresando
valores en el intervalo [-5, 5].
Solucin:
Conocido el intervalo [a, b], la expresin a implementarse es:
1
i =1
F ( x) = sen( x 2 )dx =
(1) i +1
(b 4i 1 a 4i 1 )
(2i 1)!(4i 1)
Programa en Matlab:
% IntegralSerieTaylor.m
% Calculo de integrales usando series de Taylor, para f(x) = sin(x^2).
clc; % Limpia pantalla
clear; % Libera espacio de variables
disp('Calculo de integrales usando series de Taylor');
while 1
% Ingreso de valores:
a = input('Ingresar el valor inicial (salir < -5 o > 5): ');
if (abs(a) > 5) break; end;
b = input(sprintf('Ingresar el valor final (salir < -5 o > 5 o < %f): ', a));
if (abs(b) > 5 | b < a) break; end;
n = input('Ingresar el nmero de trminos (salir < 1 o > 10): ');
296
297
n 1 xi +1
i = 0 xi
xi +1
n 1
n 1
f ( x)dx = y dx = y dx = y ( x
i =0 xi
i =0
xi
i =0
i +1
n 1
n 1
i =0
i =0
xi ) = y i h = h y i
= h( y0 + y1 + y 2 + ... + y n1 ); (7.4)
Cabe destacar que yi sale de la integral, debido a que es una constante para un
rectngulo.
Regla trapezoidal:
Se utilizan trapecios en lugar de
rectngulos. Para utilizar esta regla se
debe deducir la ecuacin de la recta
que se indica en la figura 7.2, en este
caso,
y = f(x) = [(yi+1 yi) / h] x + yi (7.5)
Sustituyendo 7.5 en 7.3:
b
n 1 xi +1
i =0
f ( x)dx =
= h(
n 1
xi +1
x2
h
{[( yi +1 yi ) / h] + yi x}| = ( yi +1 + yi )
x {[( yi+1 yi ) / h]x + yi }dx =
xi
2
i =0
i =0 2
i
n 1
y0 + yn
+ y1 + y 2 + ... + y n1 ); (7.6)
2
298
Regla de Simpson:
Usando
dos
trapecios
adyacentes, se tiene lo que se
indica en la figura 7.3, para
tres puntos, debe plantearse la
siguiente ecuacin:
y = Ax2 + Bx + C (7.7).
Reemplazando los valores de
los puntos de los puntos en
(7.7), se obtiene el siguiente
sistema de tres ecuaciones por
tres incgnitas:
yi = C; (7.8)
2
yi +1 = Ah + Bh + C; (7.9)
2
yi + 2 = 4 Ah + 2 Bh + C; (7.10)
Resolviendo el sistema:
yi + 2 2 yi +1 + yi
; (7.11)
A =
2h 2
yi + 2 + 4 yi +1 3 yi
; (7.12)
B =
2
2
h
C = y i ; (7.13)
n / 2 xi + 2
n/2
xi + 2
x3
x2
h3
+
B
+
Cx
)
=
(
8
A
+ 2 Bh 2 + 2Ch); (7.14)
|
a
x
i
3
2
3
i =0
i = 0 xi
i =0
Sustituyendo 7.11, 7.12 y 7.13 en 7.14 y simplificando:
b
n/2
f ( x)dx = ( Ax 2 + Bx + C )dx = ( A
n/2
h
1
4
1
f ( x)dx = h( yi + 2 + yi +1 + yi ) = ( y0 + 4 y1 + 2 y2 + 4 y3 + ... + 2 yn 2 + 4 yn 1 + yn ); (7.15)
3
3
3
3
i =0
Ejemplo 7.3
1
299
300
b = input(sprintf('Ingresar el valor final (salir < -5 o > 5 o < %f): ', a));
if (abs(b) > 5 | b < a) break; end;
n = input('Ingresar el nmero de elementos (salir < 3 o > 1000): ');
if (n < 3 | n > 1000) break; end;
% Clculo de la longitud de paso:
h = (b - a) / n;
% Clculo de la integral usando la regla del rectngulo:
vinterec = 0;
x = a;
for i = 1 : n
vinterec = vinterec + eval(f);
x = x + h;
end;
vinterec = h * vinterec;
% Clculo de la integral usando la regla del trapecio:
x = a;
vintetra = eval(f) / 2;
x = x + h;
for i = 1 : n - 1
vintetra = vintetra + eval(f);
x = x + h;
end;
vintetra = h * (vintetra + eval(f) / 2);
% Clculo de la integral usando la regla de Simpson:
x = a;
vinteSim = eval(f);
x = x + h;
for i = 1 : n - 1
fac = 2 * (mod(i, 2) + 1);
vinteSim = vinteSim + fac * eval(f);
x = x + h;
end;
vinteSim = h * (vinteSim + eval(f)) / 3;
% Clculo de la integral usando Matlab:
integrando = @fsin; % Se hacemos manipulable la funcin f(x) = sin(x^2).
intMat = quad(integrando, a, b); % Calcula la integral
% Calculo de los errores absoluto y relativo y despliegue de resultados:
era = abs(intMat - vinterec);
err = era / abs(intMat) * 100;
disp(sprintf('Integral rectangulos = %f, Integral Matlab = %f, Error absoluto = %f,
Error relativo = %f%%', vinterec, intMat, era, err));
era = abs(intMat - vintetra);
err = era / abs(intMat) * 100;
disp(sprintf('Integral trapecios = %f, Integral Matlab = %f, Error absoluto = %f, Error
relativo = %f%%', vintetra, intMat, era, err));
era = abs(intMat - vinteSim);
301
j =0
Ejemplo 7.5
Deducir la regla de Simpson, usando la frmula 7.16.
Solucin:
De acuerdo a la figura 7.3, para la deduccin, se requieren tres puntos, entonces:
p0 ( x) =
1
x x1 x x2
x h x 2h
.
=
.
= 2 ( x 2 3hx + 2h 2 )
x0 x1 x0 x2 (h) (2h) 2h
2h
1 2h 2
1 x 3 3hx 2
2
(
3
+
2
)
=
(
+ 2h 2 x )|
x
hx
h
dx
2
2
0
0
a
2h 0
2h 3
2
3
3
1 8h
1 2h
h
= 2(
6h 3 + 4h 3 ) = 2 (
) = ; (7.17 )
2h
3
2h
3
3
b
2h
p0 ( x)dx = p0 ( x)dx =
302
p1 ( x) =
1
x x0 x x2
x x 2h
.
= .
= 2 ( x 2 2hx )
x1 x0 x1 x2 h (h)
h
2h
1 2h 2
1 x3
(
2
)
=
( hx 2 )|
x
hx
dx
2
2
0
0
a
h 0
h 3
3
3
1 8h
1
4h
4h
= 2(
4h 3 ) = 2 (
)=
; (7.18)
3
3
3
h
h
1
x x0 x x1
x xh
.
=
.
= 2 ( x 2 hx )
p2 ( x) =
2h
x2 x0 x2 x1 2h (h)
b
2h
2h
p1 ( x)dx = p1 ( x)dx =
p2 ( x)dx = p2 ( x)dx =
3
1
2h 2
2h
( x 2 hx )dx =
1 x3 hx 2 2 h
(
)|
2h 2 3
2 0
1 8h
1 2h
h
(
2h 3 ) = 2 (
) = ; (7.19)
2
2h
3
2h
3
3
n/2
F ( x) = y j p j ( x) dx = (
j =0
j =0
h
4h
h
yj +
y j +1 + y j + 2 )
3
3
3
h
( y0 + 4 y1 + 2 y2 + 4 y3 + 2 y4 + ... + 2 yn 2 + 4 yn 1 + yn ); (7.20)
3
Ejemplo 7.6
Deducir la siguiente regla de Newton Cotes, con respecto a la de Simpson, usando
la frmula 7.16.
Solucin:
Para la deduccin, se requieren de cuatro puntos, entonces:
p0 ( x) =
x x1 x x2 x x3
x h x 2h x 3h
.
.
.
.
=
x0 x1 x0 x2 x0 x3 (h) (2h) (3h)
1
( x 3 6hx 2 + 11h 2 x 6h 3 )
6h 3
b
3h
1 3h
p0 ( x)dx = p0 ( x)dx = 3 ( x 3 6hx 2 + 11h 2 x 6h3 )dx
a
0
6h 0
4
2 2
3h
1 x
11h x
1 81h 4
99 h 4
= 3 ( 2hx 3 +
6h 3 x )| = 3 (
54 h 4 +
18h 4 )
0
6h 4
2
6h
4
2
4
1
9h
3h
= 3 (
)=
; (7.21)
6h
4
8
=
303
p1 ( x) =
1
x x0 x x2 x x3
x x 2h x 3h
= 3 ( x 3 5hx 2 + 6h 2 x)
= .
.
.
.
x1 x0 x1 x2 x1 x3 h ( h) ( 2h) 2h
1 x 4 5hx 3
2 2 3h
p1 ( x) dx =
+
=
(
+
3
(
5
6
)
x
hx
h
x
dx
h
x )|
3
0
0
0
a
2h 4
3
1 81h 4
1 9h 4
9h
= 3(
45h 4 + 27 h 4 ) = 3 (
)=
; (7.22)
2h
4
2h
4
8
b
1
p1 ( x) dx = 3
2h
3h
p2 ( x) =
1
x x0 x x1 x x3
x x h x 3h
=
= 3 ( x3 4hx 2 + 3h 2 x)
.
.
.
.
2h
x2 x0 x2 x1 x2 x3 2h ( h) ( h)
3h
p2 ( x) dx = p2 ( x) dx =
=
3h
1
2h3
3h
( x3 4hx 2 + 3h 2 x) dx =
1 x 4 4hx 3 3h 2 x 2 3h
+
(
)|
0
2h3 4
3
2
1 81h 4
27 h 4
1
9h 4
9h
4
(
36
+
)
=
)=
; (7.23)
h
3
3
2h
4
2
2h
4
8
p3 ( x) =
1
x x0 x x1 x x2
x x h x 2h
.
.
= .
.
= 3 ( x 3 3hx 2 + 2h 2 x)
6h
x3 x0 x3 x1 x3 x2 3h (2h) ( h)
3h
1 3h 3
1 x4
2
2
(
3
+
2
)
=
( hx 3 + h 2 x 2 )|
x
hx
h
x
dx
3
3
0
0
a
6h 0
6h 4
4
4
1 81h
1 9h
3h
= 3(
27 h 4 + 9h 4 ) = 3 (
)=
; (7.24)
6h
4
6h
4
8
b
3h
p3 ( x) dx = p3 ( x) dx =
n/3
F ( x) = y j p j ( x)dx = (
j =0
j =0
3h
9h
9h
3h
yj +
y j +1 +
y j+2 +
y j+2 )
8
8
8
8
3h
( y0 + 3 y1 + 3 y2 + 2 y3 + 3 y4 + ... + 3 yn 2 + 3 yn 1 + yn ); (7.25)
8
F ( x) =
2h
(7 y0 + 32 y1 + 12 y2 + 32 y3 + 14 y4 + 32 y5 + 12 y6 + ... + 12 yn 2 + 32 yn 1 + 7 yn ); (7.26)
45
Ejemplo 7.7
1
304
305
% Ingreso de valores:
a = input('Ingresar el valor inicial (salir < -5 o > 5): ');
if (abs(a) > 5) break; end;
b = input(sprintf('Ingresar el valor final (salir < -5 o > 5 o < %f): ', a));
if (abs(b) > 5 | b < a) break; end;
n = input('Ingresar el nmero de elementos (salir < 3 o > 1000): ');
if (n < 3 | n > 1000) break; end;
% Clculo de la longitud de paso:
h = (b - a) / n;
% Clculo de la integral usando la regla de tres octavos:
x = a;
vinteTroc = eval(f);
x = x + h;
for i = 1 : n - 1
vinteTroc = vinteTroc + to(mod(i - 1, 3) + 1) * eval(f);
x = x + h;
end;
vinteTroc = 3 * h * (vinteTroc + eval(f)) / 8;
% Clculo de la integral usando la regla de dos cuarentaycincoavos:
x = a;
vinteDocu = 7 * eval(f);
x = x + h;
for i = 1 : n - 1
vinteDocu = vinteDocu + dc(mod(i - 1, 4) + 1) * eval(f);
x = x + h;
end;
vinteDocu = 2 * h * (vinteDocu + 7 * eval(f)) / 45;
% Clculo de la integral usando Matlab:
integrando = @fsin; % Se hacemos manipulable la funcin f(x) = sin(x^2).
intMat = quad(integrando, a, b); % Calcula la integral
% Calculo de los errores absoluto y relativo y despliegue de resultados:
era = abs(intMat - vinteTroc);
err = era / abs(intMat) * 100;
disp(sprintf('Integral Tres octavos = %f, Integral Matlab = %f, Error absoluto = %f,
Error relativo = %f%%', vinteTroc, intMat, era, err));
era = abs(intMat - vinteDocu);
err = era / abs(intMat) * 100;
disp(sprintf('Integral Dos cuarentaycincoavos = %f, Integral Matlab = %f, Error
absoluto = %f, Error relativo = %f%%\n', vinteDocu, intMat, era, err));
end;
Los resultados se presentan en la siguiente pantalla:
306
Ejemplo 7.9
Implementar un programa en Visual C#, que tenga como entradas la funcin f(x) a
integrarse, el intervalo y el nmero de elementos que forman F(x) y que calcule el
valor de la integral aplicando todas las reglas vistas de Newton - Cotes, como
alternativa calcular la integral mediante series de Taylor. Desarrollar para f(x) =
sen(x2), validar ingresando valores para el nmero de elementos en el intervalo [3,
1000].
Solucin:
Archivos fuente en Visual C#:
307
Opcin Series:
private void rbSeries_CheckedChanged(object sender, EventArgs e)
{
lvar.Text = "Nmero de trminos:";
}
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
// Calcula integrales por interpolacin de Newton de una funcin continua.
// Declaracion de variables
double vinteg, limi, limf;
int n; // Nmero de elementos.
try
{
// Transferencia de captura de datos a variables:
n = int.Parse(txtNun.Text);
// Validacin:
if (n < 3 || n > 1000) return;
308
limi = double.Parse(txtLIni.Text);
limf = double.Parse(txtLFin.Text);
// Validacin:
if (limf < limi) return;
CIntegracion.n = n;
CIntegracion.li = limi;
CIntegracion.lf = limf;
if (rbNC.Checked)
{
CIntegracion.punf = f;
vinteg = CIntegracion.RRectangulos();
listBox1.Items.Add("Integral Rectangulos = " + Math.Round(vinteg, 9).ToString()
+ ", para " + CIntegracion.n.ToString() + " elementos");
vinteg = CIntegracion.RTrapecios();
listBox1.Items.Add("Integral Trapecios = " + Math.Round(vinteg, 9).ToString()
+ ", para " + CIntegracion.n.ToString() + " elementos");
vinteg = CIntegracion.RSimpson();
listBox1.Items.Add("Integral Simpson = " + Math.Round(vinteg, 9).ToString()
+ ", para " + CIntegracion.nc.ToString() + " elementos");
vinteg = CIntegracion.RTresOctavos();
listBox1.Items.Add("Integral Tres octavos = " + Math.Round(vinteg, 9).ToString()
+ ", para " + CIntegracion.nc.ToString() + " elementos");
vinteg = CIntegracion.RDosCuarentaycincoavos();
listBox1.Items.Add("Integral Dos cuarentaycincoavos = " + Math.Round(vinteg, 9).ToString()
+ ", para " + CIntegracion.nc.ToString() + " elementos");
listBox1.Items.Add(" ");
}
else
{
vinteg = CIntegracion.IntSenCuad();
listBox1.Items.Add("Integral Sen(x^2) por series Taylor = " + Math.Round(vinteg, 9).ToString()
+ ", para " + CIntegracion.n.ToString() + " elementos");
listBox1.Items.Add(" ");
}
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
309
CListBox.Limpiar();
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
310
F ( x) =
3.3
f ( x)dx
1.9
Solucin:
El area a integrarse es la que se muestra en la
figura 7.5
F(X) = A1 + A2 + A3, donde:
A1 es el area en el tramo 1.9 x 2.5
A2 es el area en el tramo 2.5 x 3
A3 es el area en el tramo 3 x 3.3
Para A1 se tomarn los siete puntos, h = 0.1 y se
aplicar dos veces la regla de tres octavos.
Para A2 se tomarn los once puntos para calcular
siete ordenadas por interpolacin,
h = (3 2.5) / 6 = 1 / 12 y se aplicar dos veces
la regla de tres octavos.
Para A3 se tomarn los cuatro puntos, h = 0.1 y se aplicar una vez la regla de tres
octavos.
Desarrollando el proceso:
A1 = 3*0.1*(2.89 + 3*6.35 + 3*10.23 + 2*14.40 + 3*18.70 + 3*22.92 + 26.82) / 8 = 8.74
Al efectuar la interpolacin de los puntos de A2, se obtiene la siguiente tabla:
entonces,
A2 = 3*(1/12)*(26.82 + 3*29.52 + 3*31.90 + 2*33.22 + 3*33.59 + 3*32.68 + 30.51) / 8
= 15.84
311
Ejemplo 7.11
Implementar un programa en Matlab y otro en Visual C#, que tenga como entradas la
muestra a integrarse, y el nmero de elementos que forman F(x) y que calcule el
valor de la integral aplicando las reglas de tres octavos y de dos cuarentaicincoavos
e interpolacin mediante diferencias divididas de Newton. Desarrollar para la muestra
del ejemplo 7.8, repetir el proceso ingresando el nmero de elementos en el intervalo
[11, 1000]. El programa debe calcular apropiadamente el nmero de puntos y los
errores absoluto y relativo.
312
Solucin:
Programa en Matlab:
IntegracionXInterpolacionNewtonDis.m
% Calcula integral por interpolacin de Newton mediante diferencias divididas
% para una muestra.
clear;
clc;
% Definicin de la muestra
syms x;
F = exp(x)*cos(2*x);
vx = [1.9 2 2.1 2.2 2.3 2.4 2.5 3 3.1 3.2 3.3];
vy = [2.89 6.35 10.23 14.40 18.70 22.92 26.82 30.51 25.81 18.65 8.87];
to = [3, 3, 2]; % coeficientes regla tres octavos
dc = [32, 12, 32, 14]; % coeficientes regla dos cuarentaycincoavos
tf = length(vx);
% Proceso:
while 1
n = input(sprintf('Ingrese el nmero de elementos (salir < %d o > %d): ', tf, 1000));
% Validacin:
if (n < tf) | (n > 1000) break; end;
% regla tres octavos:
rg = 3;
% Ajuste de n:
if (mod(n, rg) ~= 0) n = n + rg - mod(n, rg); end;
h = (vx(tf) - vx(1)) / n;
x = vx(1);
vinteTroc = InterpolacionNewton(vx, vy, x, tf);
x = x + h;
for i = 1 : n - 1
vinteTroc = vinteTroc + to(mod(i - 1, rg) + 1) * InterpolacionNewton(vx, vy, x, tf);
x = x + h;
end;
vinteTroc = 3 * h * (vinteTroc + InterpolacionNewton(vx, vy, x, tf)) / 8;
% Calculo de la integral exacta:
x = vx(tf);
vintexa = eval(F);
x = vx(1);
vintexa = vintexa - eval(F);
% Calculo de los errores absoluto y relativo y despliegue de resultados:
era = abs(vintexa - vinteTroc);
err = era / abs(vintexa) * 100;
disp(sprintf('Integral Tres octavos = %f, Integral Exacta = %f, Error absoluto = %f,
Error relativo = %f%%, para %d elementos', vinteTroc, vintexa, era, err, n));
% regla dos cuarentaycincoavos:
313
rg = 4;
% Ajuste de n:
if (mod(n, rg) ~= 0) n = n + rg - mod(n, rg); end;
h = (vx(tf) - vx(1)) / n;
x = vx(1);
vinteDocu = 7*InterpolacionNewton(vx, vy, x, tf);
x = x + h;
for i = 1 : n - 1
vinteDocu = vinteDocu + dc(mod(i - 1, rg) + 1) * InterpolacionNewton(vx, vy, x, tf);
x = x + h;
end;
vinteDocu = 2 * h * (vinteDocu + 7*InterpolacionNewton(vx, vy, x, tf)) / 45;
% Calculo de los errores absoluto y relativo y despliegue de resultados:
era = abs(vintexa - vinteDocu);
err = era / abs(vintexa) * 100;
disp(sprintf('Integral Dos cuarentaycincoavos = %f, Integral Exacta = %f, Error
absoluto = %f, Error relativo = %f%%, para %d elementos\n', vinteDocu, vintexa, era,
err, n));
end;
314
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
// Calcula integrales por interpolacin de Newton de una muestra y calcula errores.
// Declaracion de variables
double vinteTroc, vinteDocu, vintexa, era, err;
int n; // Nmero de elementos.
try
{
// Transferencia de captura de datos a variables:
n = int.Parse(txtNroEl.Text);
// Validacin:
if (n < tf || n > 1000) return;
CIntegracion.n = n;
CIntegracion.li = vdx[0];
CIntegracion.lf = vdx[tf - 1];
CIntegracion.punf = InterpolacionNewton;
vinteTroc = CIntegracion.RTresOctavos();
// Calculo de la integral exacta:
vintexa = F(vdx[tf - 1]) - F(vdx[0]);
// Calculo de los errores absoluto y relativo y despliegue de resultados:
315
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Limpiar:
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
316
317
F2 ( x) =
4
Solucin:
F1 ( x) =
3 0
0.316238 0
0
0.316238
f ( x)dxdx y
f ( x)dxdxdxdx
0
0
0.316238
= (60 x 2 6)|
0
0.316238
F2 ( x) =
4
3 0
0
0
0.316238
0.316238
+ (60 x 2 )|
0
0.316238
)|
0
0.316238
= 6
0.316238 0
3
3
3
((20 x 6 x) + (20 x 6 x)| )dx = (20 x 6 x + 508)dx =
3
Ejemplo 7.13
Implementar un programa en Matlab y otro en Visual C#, que tenga como entradas
f(x) a integrarse, y el intervalo de integracin [li, lf] y que calcule el valor de la integral
aplicando integracin sucesiva. Desarrollar para los casos del ejemplo 7.12, repetir el
proceso ingresando valores de li y lf en el intervalo [-2999, 2999].
Solucin:
Programa en Matlab:
% PIntegralNumPol.m
% Prueba Integracin numrica de polinomios.
% FUNCION PRINCIPAL:
clc; % Limpia pantalla
clear; % Libera espacio de variables
global A N
disp('Integracin de polinomios: ');
[N, A] = LeePolinomio;
disp('Polinomio ingresado: ');
EscribePolinomio(A, N + 1);
while 1
li = input('Ingrese el limite inicial (salir < -2999 o > 2999): ');
if (abs(li) > 2999)
break
end;
318
319
320
// Datos miembro:
static bool pv = true;
// Funciones miembro:
321
COprPolinomios.lb = listBox1;
COprPolinomios.EscribePolinomio("Polinomio ingresado: ");
it = COprPolinomios.IntegraPolinomio(limi, limf);
listBox1.Items.Add("Valor de la integral = " + it.ToString());
Botn Limpiar:
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
322
323
324
Para el calculo del error relativo de la integral de la muestra del ejemplo 7.11, se ha
considerado que la muestra tiende a la funcin f(x) = ex(cos(2x) 2sen(2x)), cuya
antiderivada es F(x) = excos(2x).
Analizando los resultados puede apreciarse de que el error cometido al aplicar
integracin por interpolacin de Newton y las reglas de tres octavos y dos
cuarentaicincoavos para un nmero de elementos pequeo no supera el 0.053% y
para un n grande de 0.046%, en todo caso la regla de dos cuarentaicincoavos es
ms precisa.
325
2.- Modificar los programas de los ejemplos 7.4 y 7.8 para implementar el ejercicio 1,
desplegar resultados ingresando algunos intervalos.
326
4.- Implementar un programa en Matlab, que tenga como entradas la funcin f(x) a
integrarse, el intervalo y el nmero de trminos de la serie de Taylor y que calcule el
valor de la integral, adems el error absoluto y relativo entre el valor que produce
Matlab y el calculado. Desarrollar para el f(x) del ejercicio 3, repetir el proceso
ingresando valores del lmite inicial y final en el intervalo [-7, 7] y del nmero de
trminos en el intervalo [3, 15].
327
5.- Dada la siguiente muestra de ventas (en miles de USD) de un producto en ocho
meses del ao:
Mes
Venta
1
85
2
93
3
61
4
73
5
110
9
91
10
60
12
87.5
6.- Implementar un programa en Matlab y otro en Visual C#, que tenga como
entradas la muestra a integrarse, y el nmero de elementos que forman F(x) y que
calcule el valor de la integral aplicando las reglas de tres octavos y de dos
cuarentaicincoavos e interpolacin mediante los polinomios de Lagrange. Desarrollar
para la muestra del ejercicio 5, repetir el proceso ingresando el nmero de elementos
en el intervalo [11, 1000]. El programa debe calcular apropiadamente el nmero de
puntos.
328
F2 ( x) =
0
4 3 2
4 1
f ( x)dxdxdxdx
8.- Implementar un programa en Matlab y otro en Visual C#, que tenga como
entradas la funcin a integrarse f(x), y el nmero de elementos que forman F(x) y
que calcule el valor de la integral aplicando el mtodo de Romberg. Desarrollar para
la funcin del ejercicio 3, repetir el proceso ingresando el nmero de elementos en el
intervalo [3, 1000]. El programa debe calcular apropiadamente el nmero de puntos.
329
330
8
RESOLUCIN DE ECUACIONES
DIFERENCIALES DE PRIMER ORDEN
331
332
CAPITULO VIII
RESOLUCION DE ECUACIONES DIFERENCIALES DE PRIMER ORDEN.
Generalidades
Las ecuaciones diferenciales de primer orden tienen su aplicacin en la fsica,
quimica, mecnica, economa entre otros.
En este captulo se resolvern ecuaciones diferenciales de primer orden de la forma:
f'(x) = g(x; f(x)), (x0; f(x0)), y = f(x) (8.1)
Ejemplos:
1) y' = (x - y) / (x + y), (x0; f(x0))
2) y' = -3y / (2x), (1; 1)
3) y' = x / y, (0; 2)
4) y' = y, (2; e2)
La frmula (8.1), tambin puede adoptar la siguiente forma:
e(x, y)dx + h(x, y)dy = 0, (x0; f(x0)), y = f(x) (8.2)
Ejemplos:
1) (5x2y - xy)dx + (2x3y2 + x3y4)dy = 0, (x0; f(x0))
2) ydx - xdy = 0, (1; 2)
3) excos(2x)ydx + y3dy = 0, (x0; f(x0))
Los temas que se revisarn en este captulo son:
yk + 1 = yk + hg(xk ; yk ) (8.6)
333
4 + x 2 ; f(2) =
334
4 + 22 =
8 = 2.82842712
Mtodo de Taylor:
Se basa en la utilizacin de la frmula (6.2), en este caso se puede sustituir en la
frmula anterior f(x+h) por yk+1, f(x) por yk, de manera que:
yk +1 = yk + hf ' ( x) +
h2
h3
h4 ( 4)
f ' ' ( x) +
f ' ' ' ( x) +
f ( x) + ...(8.7)
2!
3!
4!
yk +1 = yk + hg ( x, y) +
h2
h2 g ( x, y)
g ( x, y)
)(8.8)
f ' ' ( x) = yk + hg ( x, y) + (
+ g ( x, y)
2!
2
x
y
335
Solucin:
g(x, y) = x / y; gx(x, y) = 1 / y; gy(x, y) = -x / y2; x0 = 0; y0 = 2; x1 = 2; n = 8; h = (2 0) / 8 = 0.25; el resto del proceso se detalla en la siguiente tabla:
pk + 1 = yk + hg(xk ; yk ) (8.10)
ck +1 = yk + 1 = yk +
h
(g(x k ; y k ) + g(x k +1 ; p k +1 )) (8.11)
2
336
m1 = hg(x k ; y k )
h
m
; yk + 1 )
2b
2b
y k + 1 = y k + (1 - b)m1 + bm 2 ; (8.12)
m 2 = hg(x k +
manera:
337
b) xk = x0
c) yk = y0
d) m1 = h*g(xk, yk)
e) m2 = h*g(xk + h/2/b, yk + m1/2/b)
f) yk = yk + (1 - b)*m1 + b*m2
g) xk = xk + h
h) Repetir desde el paso (d) hasta que xk >= xn
i) Desplegar yn
Ejemplo 8.4
Calcular f(2), para y' = x / y, (0; 2), usando el mtodo R-K-2 y ocho elementos.
Solucin:
g(x, y) = x / y; x0 = 0; y0 = 2; x1 = 2; n = 8; h = (2 - 0) / 8 = 0.25; b = 0.5; el resto del
proceso se detalla en la siguiente tabla:
338
m1 = hg(x k ; y k )
h
m
m 2 = hg(x k + ; y k + 1 )
2h
2
m
m 3 = hg(x k + ; y k + 2 )
2
2
m 4 = hg(x k + h; y k + m3 )
yk +1 = yk +
manera:
m1 + 2m 2 + 2m 3 + m 4
; (8.13)
6
339
% PrbEcDiferencial.m
clc; % Limpia pantalla
clear; % Libera espacio de variables
syms x y % Define como simbolo a x y
g = input('Ingresar la funcin g(x, y): ');
340
341
y = y0;
for k = 1: n
f = eval(g);
a = eval(derivadax);
b = eval(derivaday);
y = y + h * f + h^2 * (a + b * f) / 2;
x = x + h;
end;
% EcDifEulerModificado.m
function y = EcDifEulerModificado(g, x0, y0, xn, n)
h = (xn - x0) / n;
syms x y; % Define como simbolo a x y
x = x0;
y = y0;
for k = 1: n
f = eval(g);
pk1 = y + h * f;
x = x + h;
y1 = y;
y = pk1;
m = eval(g);
y = y1 + h * (f + m) / 2;
end;
% EcDifRungeKuttaDos.m
function y = EcDifRungeKuttaDos(g, x0, y0, xn, n, b)
h = (xn - x0) / n;
syms x y; % Define como simbolo a x y
x = x0;
y = y0;
for k = 1: n
m1 = h * eval(g);
x1 = x;
x = x + h/2/b;
y1 = y;
y = y + m1/2/b;
m2 = h * eval(g);
x = x1 + h;
y = y1 + (1 - b) * m1 + b * m2;
end;
% EcDifRungeKuttaCuatro.m
function y = EcDifRungeKuttaCuatro(g, x0, y0, xn, n)
h = (xn - x0) / n;
syms x y; % Define como simbolo a x y
342
x = x0;
y = y0;
for k = 1: n
m1 = h * eval(g);
x1 = x;
x = x + h/2;
y1 = y;
y = y + m1/2;
m2 = h * eval(g);
y = y1 + m2/2;
m3 = h * eval(g);
x = x1 + h;
y = y1 + m3;
m4 = h * eval(g);
y = y1 + (m1 + 2*m2 + 2*m3 + m4)/6;
end;
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
343
344
Funciones externas:
// Funcion a ser analizada:
static double g(double x, double y)
{
return x / y;
}
// Solucin de la EDO:
static double slg(double x)
{
return Math.Sqrt(Math.Pow(x, 2) + 4);
}
// Derivada parcial de g respecto de x:
static double gx(double x, double y)
{
return 1 / y;
}
// Derivada parcial de g respecto de y:
static double gy(double x, double y)
{
return -1 / Math.Pow(y, 2);
}
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
// Resuelve ecuaciones diferenciales ordinarias de primer orden.
// Declaracion de variables
double x0, y0, xn, yn, b, era, err, solexa;
int n; // Nmero de elementos.
try
{
// Transferencia de captura de datos a variables:
n = int.Parse(txtNun.Text);
345
// Validacin:
if (n < 3 || n > 1000) return;
x0 = double.Parse(txtX0.Text);
y0 = double.Parse(txtY0.Text);
xn = double.Parse(txtXn.Text);
CEcDifOrdIOrden.n = n;
CEcDifOrdIOrden.x0 = x0;
CEcDifOrdIOrden.y0 = y0;
CEcDifOrdIOrden.xn = xn;
CEcDifOrdIOrden.pung = g;
solexa = slg(xn);
if (rbTd.Checked)
{
CEcDifOrdIOrden.pungx = gx;
CEcDifOrdIOrden.pungy = gy;
yn = CEcDifOrdIOrden.Euler();
era = Math.Abs(yn - solexa);
err = era / Math.Abs(solexa) * 100;
listBox1.Items.Add("f(" + xn.ToString() + ") (Mtodo de Euler) = "
+ Math.Round(yn, 7).ToString()
+ ", n = " + n.ToString()
+ ", Error absoluto = " + Math.Round(era, 7).ToString()
+ ", Error relativo = " + Math.Round(err, 7).ToString() + "%");
yn = CEcDifOrdIOrden.Taylor();
era = Math.Abs(yn - solexa);
err = era / Math.Abs(solexa) * 100;
listBox1.Items.Add("f(" + xn.ToString() + ") (Mtodo de Taylor) = "
+ Math.Round(yn, 7).ToString()
+ ", n = " + n.ToString()
+ ", Error absoluto = " + Math.Round(era, 7).ToString()
+ ", Error relativo = " + Math.Round(err, 7).ToString() + "%");
yn = CEcDifOrdIOrden.EulerModificado();
era = Math.Abs(yn - solexa);
err = era / Math.Abs(solexa) * 100;
listBox1.Items.Add("f(" + xn.ToString() + ") (Mtodo de Euler
Modificado) = " + Math.Round(yn, 7).ToString()
+ ", n = " + n.ToString()
+ ", Error absoluto = " + Math.Round(era, 7).ToString()
+ ", Error relativo = " + Math.Round(err, 7).ToString() + "%");
}
else if (rbRk2.Checked)
{
b = double.Parse(txtB.Text);
if (b <= 0 || b >= 1) return;
CEcDifOrdIOrden.b = b;
yn = CEcDifOrdIOrden.RungeKuttaOrdenDos();
era = Math.Abs(yn - solexa);
err = era / Math.Abs(solexa) * 100;
listBox1.Items.Add("f(" + xn.ToString() + ") (Mtodo de Runge
Kutta Orden dos) = " + Math.Round(yn, 7).ToString()
+ ", n = " + n.ToString()
+ ", b = " + Math.Round(b, 7).ToString()
+ ", Error absoluto = " + Math.Round(era, 7).ToString()
+ ", Error relativo = " + Math.Round(err, 7).ToString() + "%");
}
346
else
{
yn = CEcDifOrdIOrden.RungeKuttaOrdenCuatro();
era = Math.Abs(yn - solexa);
err = era / Math.Abs(solexa) * 100;
listBox1.Items.Add("f(" + xn.ToString() + ") (Mtodo de Runge
Kutta Orden cuatro) = " + Math.Round(yn, 7).ToString()
+ ", n = " + n.ToString()
+ ", Error absoluto = " + Math.Round(era, 7).ToString()
+ ", Error relativo = " + Math.Round(err, 7).ToString() + "%");
}
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void cb_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
347
348
Para calcular el error relativo, se ha tomado como referencia la solucin exacta del
ejemplo 8.1.
De acuerdo a los resultados obtenidos, se puede concluir que para valores pequeos
de n, el mtodo ms preciso el de Runge-Kutta de orden cuatro. Para valores
grandes, se aaden a esta lista los mtodos de Runge-Kutta de orden dos y de Euler
Modificado. Para Runge-Kutta de orden dos, cuando n es pequeo, el mtodo es
ms preciso cuando b < 0.5; en cambio cuando n es grande, es ms preciso para b
= 0.5. Independientemente del mtodo con excepcin del de los rectngulos se
obtiene mayor precisin cuando se toman ms elementos en el clculo.
349
350
351
APENDICE A
EJERCICIOS EN MATLAB
Ejercicio a.1:
Uso de la sentencia if:
Implementar un programa en Matlab para ingresar los lados de un tringulo en orden
ascendente y determinar si lo es; en el caso de que lo sea, determinar y desplegar su
tipo y calcular desplegando el permetro y la superficie.
Solucin:
Se representa un tringulo en la figura a.1:
Se tiene un tringulo si
a+bc
Y para determinar su tipo:
si a = c entonces
es un tringulo equiltero
caso contrario si (a = b) o
(b = c) entonces
es un tringulo isceles
caso contrario
es un tringulo escaleno.
permetro = a + b + c; s = permetro
/ 2;
superficie = s (s - a) (s - b) (s - c)
352
elseif (a == b) || (b == c)
disp('Es un tringulo isceles');
else
disp('Es un tringulo escaleno');
end;
per = a + b + c;
s = per / 2;
sup = sqrt(s * (s - a) * (s - b) * (s - c));
disp(sprintf('Permetro = %8.4f; Superficie = %8.4f', per, sup));
else
disp('No es un tringulo');
end
end
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
353
Ejercicio a.2:
Uso de ciclos:
Implementar un programa en Matlab para ingresar el ngulo en grados
sexagesimales;
Si el ngulo se encuentra entre -180 y 180 grados, se convierte a radianes y se
calcula el coseno trigonomtrico de ese ngulo y se ingresa otro ngulo. El proceso
se repite hasta cuando se ingresa un ngulo fuera del rango indicado.
Para el clculo coseno trigonomtrico, se usa la siguiente serie de Taylor:
2
Cos(x) = 1 x / 2! + x / 4! x / 6! +
Donde x es el ngulo en radianes.
Para finalizar el clculo tomar hasta el trmino >= 1e-5.
Solucin:
Archivo fuente en MATLAB:
% coseno.m
precision = 1e-5;
while 1
a = input('Ingresar un ngulo en grados sexagesimales (salir < -180 0 > 180):');
if abs(a) > 180
break;
end;
x = pi * a / 180;
suma = 0;
signo = 1;
for par = 0 : 2: 1000
%Calcular el factorial:
fac = 1;
for i = 2 : par
fac = fac * i;
end;
termino = (x ^ par) / fac;
suma = suma + signo * termino;
signo = -signo;
if termino < precision
break;
end
end;
disp(sprintf('Coseno(%8.4f) = %8.4f', a, suma));
end
354
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
Ejercicio a.3:
Grfico en tres dimensiones:
Implementar un programa en Matlab para ingresar la base y la altura de un cilindro y
graficarlo
Solucin:
Se usa el mtodo de sobre posicin de discos elementales.
Archivo fuente en MATLAB:
% cilindro.m
% PROGRAMA PARA INGRESAR DESDE EL TECLADO LAS DIMENSIONES DE UN
% CILINDRO DE REVOLUCION (DIAMETRO DE LA BASE Y ALTURA DEL CILINDRO
% Y DIBUJAR LA PARED LATERAL DEL CILINDRO DE COLOR VERDE,
% LA BASE DE COLOR ROJO Y LA TAPA DE COLOR AZUL. SE USA EL METODO DE
% LA SOBREPOSICION DE SEMIDISCOS
rango = 100;
d = input('INGRESE EL DIAMETRO DE LA BASE: ');
h = input('INGRESE LA ALTURA DEL CILINDRO: ');
r = d/2.;
x = -r : d/rango : r;
n = length(x);
% DIBUJO DE LA CARA LATERAL FRONTAL:
y = sqrt(r*r - x.*x);
% ARTIFICIO PARA QUE SE DEFINA LA DIMENSION DE z:
355
i = 0;
for j = 1 : n
z(j) = i;
i = i + h/rango;
end;
for i = 0 : h/rango : h
for j = 1 : n
z(j) = i;
end;
plot3(x, y, z, 'g');
hold on;
end;
% DIBUJO DE LA CARA LATERAL POSTERIOR:
y = -sqrt(r*r - x.*x);
for i = 0 : h/rango : h
for j = 1 : n
z(j) = i;
end;
plot3(x, y, z, 'g');
hold on;
end;
% DIBUJO DE LA BASE:
for i = -r : d/rango : r
plot3([i, i], [-sqrt(r*r - i*i), sqrt(r*r - i*i)], [0, 0], 'r');
hold on;
end;
% DIBUJO DE LA TAPA:
for i = -r : d/rango : r
plot3([i, i], [-sqrt(r*r - i*i), sqrt(r*r - i*i)], [h, h], 'b');
hold on;
end;
xlabel('EJE X');ylabel('EJE Y');zlabel('EJE Z')
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
356
Ejercicio a.4:
Grfico en tres dimensiones:
Implementar un programa en Matlab para ingresar el dimetro de una esfera y
graficarla
Solucin:
Se usa el mtodo de sobre posicin de discos elementales.
Archivo fuente en MATLAB:
% esfera.m
% PROGRAMA PARA INGRESAR DESDE EL TECLADO EL DIAMETRO DE UNA
% ESFERA Y DIBUJAR CUATRO PARTES DE LA MISMA CON LOS COLORES,
% DE LA BANDERA NACIONAL. SE USA EL METODO DE LA SOBREPOSICION DE
% DISCOS. ANTES DE PLOTEAR, ALMACENA EN ARREGLOS
resolhor = 100;
resolvert = 50;
d = input('INGRESE EL DIAMETRO DE LA ESFERA: ');
r = d/2.;
% GRAFICACION DE LA PARTE INFERIOR
k = 0;
for i = -r : d/resolvert : -r/2.
rp = sqrt(r*r - i*i);
l = 0;
k = k + 1;
for j = -rp : 2*rp/resolhor : rp
l = l + 1;
x(k, l, 1) = j;
357
x(k, l, 2) = j;
y(k, l, 1) = -sqrt(rp*rp - j*j);
y(k, l, 2) = sqrt(rp*rp - j*j);
z(k, l, 1) = i;
z(k, l, 2) = i;
end;
end;
plot3(x, y, z, 'r');
hold on;
% GRAFICACION DE LA PARTE CENTRAL
k = 0;
for i = -r/2. : d/resolvert : 0.
rp = sqrt(r*r - i*i);
l = 0;
k = k + 1;
for j = -rp : 2*rp/resolhor : rp
l = l + 1;
x1(k, l, 1) = j;
x1(k, l, 2) = j;
y1(k, l, 1) = -sqrt(rp*rp - j*j);
y1(k, l, 2) = sqrt(rp*rp - j*j);
z1(k, l, 1) = i;
z1(k, l, 2) = i;
end;
end;
plot3(x1, y1, z1, 'b');
hold on;
% GRAFICACION DE LA PARTE SUPERIOR
k = 0;
for i = 0 : d/resolvert : r
rp = sqrt(r*r - i*i);
l = 0;
k = k + 1;
for j = -rp : 2*rp/resolhor : rp
l = l + 1;
x2(k, l, 1) = j;
x2(k, l, 2) = j;
y2(k, l, 1) = -sqrt(rp*rp - j*j);
y2(k, l, 2) = sqrt(rp*rp - j*j);
z2(k, l, 1) = i;
z2(k, l, 2) = i;
end;
end;
plot3(x2, y2, z2, 'y');
hold on;
xlabel('EJE X');ylabel('EJE Y');zlabel('EJE Z')
Resultados obtenidos:
Se presenta a continuacin los resultados que produce la ejecucin de este
programa:
358
359
APENDICE B
EJERCICIOS EN VISUAL C#.NET
Ejercicio b.1:
Uso de la sentencia if:
Implementar el programa del ejercicio a.1
Solucin:
Archivos fuente en Visual C#:
Antes
de
describir
el
cdigo fuente,
en Visual C#,
se
ha
implementado
un
proyecto
tipo Windows
Application, la
cual contiene
el diseo de la
ventana de la
figura b.1
Como
se
puede
apreciar,
existen
seis
cuadros
de
texto, tres que
permiten
ingresar
los
lados del tringulo en orden ascendente y los otros tres para desplegar los
resultados; y dos botones, mismos que tienen la funcionalidad:
- Procesar: valida los lados del tringulo, determina y despliega el tipo de tringulo;
calcula y despliega el permetro y superficie.
- Salir: Finaliza la ejecucin del programa.
360
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
// Declaracion de variables
double a, b, c; // Lados del triangulo
double perimetro, superfifie, s;
try
{
a = double.Parse(txtLado1.Text);
b = double.Parse(txtLado2.Text);
c = double.Parse(txtLado3.Text);
if (a + b >= c)
{ // Es un triangulo
// Determinar el tipo:
if (a == c)
txtTipo.Text = "EQUILATERO";
else if (a == b || b == c)
txtTipo.Text = "ISOCELES";
else
txtTipo.Text = "ESCALENO";
// Calcular el perimetro y superficie:
perimetro = a + b + c;
s = perimetro / 2;
superfifie = Math.Sqrt(s * (s - a) * (s - b) * (s - c));
txtPerimetro.Text = perimetro.ToString();
txtSuperficie.Text = superfifie.ToString();
}
else
txtTipo.Text = "NO ES UN TRIANGULO";
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
361
Ejercicio b.2:
Uso de ciclos:
Implementar el programa del ejercicio a.2
Solucin:
362
363
Botn Procesar:
private void b_procesar_Click(object sender, EventArgs e)
{
// Declaracion de variables
double a; // Angulo en grados
double x; // Angulo en radianes
double resfuncion; // Resultado de la funcion
try
{
// Transferir del textbox a memoria:
a = double.Parse(txtAGrados.Text);
// Convertir el angulo de grados sexagecimales a radianes:
x = Math.PI * a / 180;
// Calcular el coseno:
resfuncion = coseno(x);
// Desplegar resultados:
txtARadianes.Text = x.ToString();
txtFuncion.Text = resfuncion.ToString();
txtNTer.Text = numterminos.ToString();
}
catch (Exception ex)
{
MessageBox.Show("Error " + ex.Message);
}
}
Botn Salir:
364
365
Ejercicio b.3:
Uso de clases:
Implementar una versin del programa del ejercicio b.2, que utilice orientacin a
objetos.
Solucin:
Archivos
fuente
Visual C#:
en
Antes
de
describir
el
cdigo fuente,
en Visual C#,
se
ha
implementado
un
proyecto
tipo Windows
Application, la
cual contiene
366
367
368
return;
}
listBox1.Items.Add(sfun + "(" + ang.ToString() + ") = " + vfun.ToString() + ", Nmero de trminos
= " + nt.ToString());
}
// Cerrar flujo de entrada:
objes.CerrarLectura();
}
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}
Botn Limpiar:
private void b_limpiar_Click(object sender, EventArgs e)
{
CListBox.lb = listBox1;
CListBox.Limpiar();
}
Botn Salir:
private void b_salir_Click(object sender, EventArgs e)
{
this.Close();
this.Dispose();
}
369
370
371
class CEntSalArchivos
{
// Datos Miembro:
private String NPath, NomArchEnt, NomArchSal;
private StreamReader stLectura;
private StreamWriter stEscritura;
// Metodos:
// Constructor
public CEntSalArchivos(String p)
{
NPath = p;
}
372
373
}
// Metodo para leer toda la linea desde stLectura
public String LeerLinea()
{
return stLectura.ReadLine();
}
// Metodo para escribir toda la linea sobre stEscritura
public void EscribirLinea(String s)
{
stEscritura.WriteLine(s);
}
// Metodo para escribir sin saltar la linea sobre stEscritura
public void Escribir(String s)
{
stEscritura.Write(s);
}
// Metodo para cerrar el flujo de entrada:
public void CerrarLectura()
{
stLectura.Close();
}
class CListBox
{
// Datos miembro:
static public ListBox lb;
// Funciones miembro:
static public void Limpiar()
{
// Recuperar desde el ultimo elemento del listbox e ir eliminando:
int n = lb.Items.Count;
for (int i = n - 1; i >= 0; i--)
{
lb.Items.Remove(lb.Items[i].ToString());
}
}
static public void GuardarComo(String spa)
374
try
{
}
catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
375
xd += dx;
YD = punf(xd);
}
while ((P = YI * YD) > 0 && xd < b);
if (P == 0)
if (YI == 0)
lb.Items.Add("RAIZ EXACTA (LOCALIZACION CAMBIO DE SIGNO) = " +
xi.ToString());
else
{
lb.Items.Add("RAIZ EXACTA (LOCALIZACION CAMBIO DE SIGNO) = " +
xd.ToString());
xd += dx;
}
return (P);
}
// Funcion que permite calcular una raiz por el metodo de la biseccion:
static public double Biseccion()
{
int CITER = 0;
double Xm, Ym, YI;
String m = "";
YI = punf(xi);
lb.Items.Add("METODO DE LA BISECCION, RAIZ ENTRE " + xi.ToString() + " Y " +
xd.ToString());
do
{
Xm = (xi + xd) / 2;
Ym = punf(Xm);
if (Ym * YI > 0)
{
xi = Xm;
YI = Ym;
}
else xd = Xm;
}
while ((Ym = Math.Abs(Ym)) > precision && ++CITER < maxIter);
if (Ym > precision)
m = "NO SE ALCANZA LA PRESICION REQUERIDA: " + Ym.ToString();
m += "RAIZ = " + Xm.ToString() + ", EN " + CITER.ToString() + " ITERACIONES";
lb.Items.Add(m);
return (Xm);
}
// Funcion que permite calcular una raiz por el metodo de aproximaciones sucesivas:
static public double AproxSucesivas()
{
int CITER = 0;
double Xm, Ym = (xi + xd) / 2, FP = 0, Prec;
376
String m = "";
lb.Items.Add("METODO DE APROXIMACIONES SUCESIVAS, RAIZ ENTRE " +
xi.ToString() + " Y " + xd.ToString());
do
{
Xm = Ym;
Ym = punfm(Xm);
}
while ((Prec = Math.Abs(Ym - Xm)) > precision && (FP = Math.Abs(punfpm(Xm))) <= 1
&& ++CITER < maxIter);
if (FP > 1)
m = "NO CUMPLE CON LA CONDICION DE CONVERGENCIA";
else
{
if (Prec > precision)
m = "NO SE ALCANZA LA PRESICION REQUERIDA: " + Prec.ToString();
m += " RAIZ = " + Xm.ToString() + " EN " + CITER.ToString() + " ITERACIONES";
}
lb.Items.Add(m);
return (Xm);
}
// Funcion que permite calcular una raiz por el metodo de aproximaciones sucesivas
modificado:
static public double AproxSucesivasMod()
{
int CITER = 0;
double YD = punfm(xi), YI, TGT, BETA;
String m = "";
lb.Items.Add("METODO DE APROXIMACIONES SUCESIVAS MODIFICADO, RAIZ
ENTRE " + xi.ToString() + " Y " + xd.ToString());
do
{
YI = YD;
YD = punfm(xd);
TGT = (YD - YI) / (xd - xi);
BETA = 1 / (1 - TGT);
xi = xd;
xd += BETA * (YD - xd);
}
while ((TGT = Math.Abs(xd - xi)) > precision && ++CITER < maxIter);
if (TGT > precision)
m = "NO SE ALCANZA LA PRESICION REQUERIDA: " + TGT.ToString();
m += " RAIZ = " + xd.ToString() + " EN " + CITER.ToString() + " ITERACIONES";
lb.Items.Add(m);
return (xd);
}
// Funcion que permite calcular una raiz por el metodo de Newton - Raphson:
static public double NewtonRaphson()
377
int CITER = 0;
double Y;
String m = "";
lb.Items.Add("METODO DE NEWTON RAPHSON, RAIZ ENTRE " + xi.ToString() + " Y "
+ xd.ToString());
xi = (xi + xd) / 2;
do
xi -= (Y = punf(xi)) / punfp(xi);
while ((Y = Math.Abs(Y)) > precision && ++CITER < maxIter);
if (Y > precision)
m = "NO SE ALCANZA LA PRESICION REQUERIDA: " + Y.ToString();
m += " RAIZ = " + xi.ToString() + " EN " + CITER.ToString() + " ITERACIONES";
lb.Items.Add(m);
return (xi);
}
// Funcion que permite calcular una raiz por el metodo de Regula - Falsi:
static public double RegulaFalsi()
{
int CITER = 0;
double YI = punf(xi),
YD = punf(xd),
Xm, Ym;
String m = "";
lb.Items.Add("METODO DE REGULA FALSI, RAIZ ENTRE " + xi.ToString() + " Y " +
xd.ToString());
do
{
Xm = xi - YI * (xd - xi) / (YD - YI);
Ym = punf(Xm);
if (YI * Ym > 0)
{
xi = Xm;
YI = Ym;
}
else
{
xd = Xm;
YD = Ym;
}
}
while ((Ym = Math.Abs(Ym)) > precision && ++CITER < maxIter);
if (Ym > precision)
m = "NO SE ALCANZA LA PRESICION REQUERIDA: " + Ym.ToString();
m += " RAIZ = " + Xm.ToString() + " EN " + CITER.ToString() + " ITERACIONES";
lb.Items.Add(m);
return (Xm);
}
378
379
x0 = x1;
x1 = x2;
x2 = x3;
y0 = y1;
y1 = y2;
y2 = y3;
}
y3 = Math.Abs(y3);
if (y3 > precision)
m = "No se alcanza la precisin requerida: " + y3.ToString();
m += "RAIZ = " + x3.ToString() + ", EN " + CITER.ToString() + " ITERACIONES";
lb.Items.Add(m);
return (x3);
}
class CLectEscrMatrices
{
// Datos Miembro:
static public int n; // Dimension del vector, matriz o sistema
static public int m; // Ancho de banda de un sistema simetrico
static public double[][] a; // Matriz de coeficientes
static public double[] c; // vector matriz de coeficientes
static public double[] b; // Vector de coeficientes o terminos independientes
static public int[] Dir; // vector de direcciones de la matriz de coeficientes (Skyline)
static public CEntSalArchivos obent; // Objeto de flujo de entrada
static public ListBox lb; // listbox para despliegue de resultados
// Metodos:
static public void LeeVector()
{
for (int i = 0; i < n; i++)
try
{
b[i] = double.Parse(obent.LeerLinea());
}
catch { }
}
static public void LeeMatriz()
{
int j;
for (int i = 0; i < n; i++)
for (j = 0; j < n; j++)
try
{
a[i][j] = double.Parse(obent.LeerLinea());
}
catch { }
380
}
static public void EscribeVector(String ob)
{
String cad = "";
lb.Items.Add("VECTOR " + ob);
for (int i = 0; i < n; i++)
cad += Math.Round(b[i], 8).ToString() + " ";
lb.Items.Add(cad);
}
static public void EscribeMatriz(String ob)
{
String cad;
int i, j;
lb.Items.Add("MATRIZ " + ob);
for (i = 0; i < n; i++)
{
cad = "";
lb.Items.Add("FILA " + (i + 1).ToString());
for (j = 0; j < n; j++)
cad += a[i][j].ToString() + " ";
lb.Items.Add(cad);
}
}
public static void LeeMatrizSimetrica()
{
int i, j;
for (i = 0; i < n; i++)
for (j = 0; j <= i; j++)
try
{
c[idc.ind(i, j)] = double.Parse(obent.LeerLinea());
}
catch { }
}
public static void EscribeMatrizSimetrica(String ob)
{
String cad;
int i, j;
lb.Items.Add("MATRIZ " + ob);
for (i = 0; i < n; i++)
{
lb.Items.Add("FILA " + (i + 1).ToString());
cad = "";
for (j = 0; j < n; j++)
if (i >= j)
cad += c[idc.ind(i, j)].ToString() + " ";
381
else
cad += c[idc.ind(j, i)].ToString() + " ";
lb.Items.Add(cad);
382
}
else if (Aux != 0) {
c[++ Cont] = Aux;
PrimerCeroDeFila = 0;
}
else if (PrimerCeroDeFila == 0)
c[++ Cont] = Aux;
}
catch { }
383
384
return b;
385
}
na = gd;
EscribePolinomio("Polinomio resultante:");
// Evaluar para x al polinomio resultante:
return EvaluaPolinomio(gd, x);
class CIntegracion
{
// Datos miembro:
public static int n, nc; // Numero de elementos o de trminos
public static double li, lf; // Lmite inicial y final de la integral
private static int[] to = { 3, 3, 2 }; // coeficientes regla tres octavos
private static int[] dc = { 32, 12, 32, 14 }; // coeficientes regla dos cuarentaycincoavos
static public pf punf;
// Metodos:
// Integral de sen(x^2) usando series de Taylor:
public static double IntSenCuad()
{
double val = 0, fac;
int i, j, signo = 1;
for (i = 1; i <= n; i ++) {
fac = 1;
for (j = 2; j <= 2*i - 1; j ++)
fac = fac * j;
val += signo * (Math.Pow(lf, 4*i - 1) - Math.Pow(li, 4*i - 1)) / fac / (4*i - 1);
signo = -signo;
}
return val;
}
// Regla de los rectngulos:
public static double RRectangulos()
{
double h = (lf - li) / n, x, val = 0;
x = li;
for (int i = 0; i < n; i++)
{
val += punf(x);
x += h;
}
return val * h;
}
// Regla de los trapecios:
public static double RTrapecios()
{
386
// Regla de Simpson:
public static double RSimpson()
{
double h, x, val;
int neb = 2;
nc = n;
if (n % neb != 0) nc = n + neb - n % neb;
h = (lf - li) / nc;
x = li;
val = punf(x);
x += h;
for (int i = 1; i < nc; i++)
{
val += neb * (i % neb + 1) * punf(x);
x += h;
}
return (val + punf(x)) * h / 3;
}
// Regla Tres octavos:
public static double RTresOctavos()
{
double h, x, val;
int neb = 3;
nc = n;
if (n % neb != 0) nc = n + neb - n % neb;
h = (lf - li) / nc;
x = li;
val = punf(x);
x += h;
for (int i = 0; i < nc - 1; i++)
{
val += to[i % neb] * punf(x);
x += h;
}
return 3 * (val + punf(x)) * h / 8;
}
387
388
a = pungx(xk, yk);
b = pungy(xk, yk);
yk += h * f + Math.Pow(h, 2) * (a + b * f) / 2;
xk += h;
}
return yk;
389
}
}
class CRaicesRealesComplejas
{
// Datos miembro:
public static int n; // Grado del polinomio
public static double[] A; // Coeficientes del polinomio
public static double[] B; // Coeficientes del polinomio
static public int maxIter; // Maximo de iteraciones
static public double precision; // Precisin requerida
static public ListBox lb; // Lisbox para despliegue de resultados
// Metodos:
// Metodo de Newton - Bairstow:
public static void NewtonBairstow()
{
int CITER = 0, I;
double AMax = Math.Abs(A[0]), CTS, C1 = 0, C2, C3, DT, DP,
DQ, DPMin = 0, P1, Q1;
/* CALCULAR LAS COTAS DEL MODULO DE LA RAIZ */
for (I = 1; I <= n; I++)
if (AMax < (P1 = Math.Abs(A[I]))) AMax = P1;
CTS = (Math.Abs(A[0]) + AMax) / Math.Abs(A[0]);
/* CALCULAR LOS COEFICIENTES INICIALES P1 Y Q1 */
P1 = (CTS > 100) ? A[1] / A[0] : A[n - 1] / A[n - 2];
Q1 = (CTS > 100) ? A[2] / A[0] : A[n] / A[n - 2];
do {
/* DEFLACIONAR EL POLINOMIO AL GRADO 2*M-2: */
B[0] = A[0]; B[1] = A[1] - P1 * B[0];
for (I = 2; I <= n; I++)
B[I] = A[I] - Q1 * B[I - 2] - P1 * B[I - 1];
/* CALCULAR C1, C2 Y C3: */
C3 = B[0]; C2 = B[1] - P1 * B[0];
for (I = 2; I < n; I++) {
C1 = B[I] - Q1 * C3 - P1 * C2;
if (I < n - 1) {
C3 = C2; C2 = C1;
}
}
/* RESOLVER EL SISTEMA:
C1*DP+C2*DQ=B[n]
C2*DP+C3*DQ=B[n-1] */
DT = C1 * C3 - C2 * C2; DP = (B[n] * C3 - B[n - 1] * C2) / DT;
DQ = (C1 * B[n - 1] - C2 * B[n]) / DT;
/* CALCULAR LOS NUEVOS COEFICIENTES P Y Q: */
P1 = P1 + DP; Q1 = Q1 + DQ;
DPMin = (DPMin < Math.Abs(DQ)) ? DPMin = Math.Abs(DQ) : Math.Abs(DP);
390
}
while (DPMin > precision && ++CITER < maxIter);
if (DPMin > precision)
lb.Items.Add("NO SE ALCANZA LA PRECISION REQUERIDA: " + DPMin.ToString());
/* VERIFICAR SI LAS RAICES SON REALES O COMPLEJAS: */
Q1 = P1 * P1 - 4 * Q1;
if (Q1 < 0)
lb.Items.Add("RAICES COMPLEJAS = " + (-P1 / 2).ToString() + " +/- " +
(Math.Sqrt(-Q1) / 2).ToString() + " * I");
else
lb.Items.Add("RAICES REALES MULTIPLES O CERCANAS: " +
((-P1 - Math.Sqrt(Q1)) / 2).ToString() + " Y " +
((-P1 + Math.Sqrt(Q1)) / 2).ToString());
lb.Items.Add("EN " + CITER.ToString() + " ITERACIONES");
/* TRANSFERIR LOS COEFICIENTES B[I] A LOS A[I]: */
for (I = 0 ; I < n - 1; I++) A[I] = B[I];
n -= 2;
}
class CResolSisSimetricosEcuaciones
{
// Datos miembro:
public static int n; // Nmero de ecuaciones = Nmero de incgnitas
public static double[] A; // matriz de coeficientes
public static double[] B; // Vector de trminos independientes
// metodos:
/* FUNCION QUE RESUELVE SISTEMAS DE ECUACIONES LINEALES
SIMETRICOS, POR EL METODO DE CROUT */
public static double CroutSimetrico(Boolean ps)
{
int i, j, k;
double det = 1;
if (ps)
{ // Primera solucion, transformar la matriz y calcular el determinante
det = A[idc.ind(0, 0)];
for (i = 1; i < n; i++)
{
for (j = 1; j <= i; j++)
for (k = 0; k <= j - 1; k++)
A[idc.ind(i, j)] -= A[idc.ind(i, k)] * A[idc.ind(j, k)] / A[idc.ind(k, k)];
det *= A[idc.ind(i, i)];
}
}
// Calculo del vector C:
for (i = 0; i < n; i++)
{
for (k = 0; k <= i - 1; k++) B[i] -= A[idc.ind(i, k)] * B[k];
B[i] /= A[idc.ind(i, i)];
391
}
// Clculo del vector X:
for (i = n - 2; i >= 0; i--)
for (k = i + 1; k < n; k++) B[i] -= A[idc.ind(k, i)] * B[k] / A[idc.ind(i, i)];
return det;
class CResolSisSimBanda
{
// Datos miembro:
public static int n; // Nmero de ecuaciones = Nmero de incgnitas
public static int m; // Ancho de banda
public static double[] A; // matriz de coeficientes
public static double[] B; // Vector de trminos independientes
// metodos:
/* FUNCION QUE RESUELVE SISTEMAS DE ECUACIONES LINEALES
SIMETRICOS EN BANDA, POR EL METODO DE CROUT */
public static double CroutSimBanda(Boolean ps)
{
int i, j, k, l;
double det = 1;
if (ps)
{ // Primera solucion, transformar la matriz y calcular el determinante
det = A[idc.inb(1, 1, m)];
for (i = 2; i <= n; i++)
{
l = 2;
if (i > m) l = i - m + 2;
for (j = 1; j <= i; j++)
for (k = l - 1; k <= j - 1; k++)
A[idc.inb(i, j, m)] -= A[idc.inb(i, k, m)] * A[idc.inb(j, k, m)] / A[idc.inb(k, k, m)];
det *= A[idc.inb(i, i, m)];
}
}
// Calculo del vector C:
for (i = 1; i <= n; i++)
{
l = 1;
if (i > m) l = i - m + 1;
for (k = l; k <= i - 1; k++) B[i - 1] -= A[idc.inb(i, k, m)] * B[k - 1];
B[i - 1] /= A[idc.inb(i, i, m)];
}
// Clculo del vector X:
for (i = n - 1; i >= 1; i--) {
l = m + i - 1;
if (i > n - m) l = n;
for (k = i + 1; k <= l; k++) B[i - 1] -= A[idc.inb(k, i, m)] * B[k - 1] / A[idc.inb(i, i, m)];
}
392
return det;
class CResolSisSimSkyline
{
// Datos miembro:
public static int n; // Nmero de ecuaciones = Nmero de incgnitas
public static int[] Dir; // Vector de direcciones
public static double[] A; // matriz de coeficientes
public static double[] B; // Vector de trminos independientes
// metodos:
/* FUNCION QUE RESUELVE SISTEMAS DE ECUACIONES LINEALES
SIMETRICOS SKYLINE, POR EL METODO DE CROUT */
public static double CroutSimSkyline(Boolean ps)
{
int i, j, k, lim;
double det = 1;
if (ps)
{ // Primera solucion, transformar la matriz y calcular el determinante
det = A[idc.insk(1, 1, Dir)];
for (i = 2; i <= n; i++)
{
lim = i - Dir[i] + Dir[i - 1] + 1;
for (j = lim + 1; j <= i; j ++)
for (k = lim; k <= j - 1; k ++)
A[idc.insk(i, j, Dir)] -= A[idc.insk(i, k, Dir)] *
CLectEscrMatrices.RecuSkyLine(j, k) / A[idc.insk(k, k, Dir)];
det = det * A[idc.insk(i, i, Dir)];
}
}
// Calculo del vector C:
for (i = 1; i <= n; i++)
{
for (j = i - Dir[i] + Dir[i-1] + 1; j <= i - 1; j ++)
B[i-1] -= A[idc.insk(i, j, Dir)] * B[j-1];
B[i-1] /= A[idc.insk(i, i, Dir)];
}
// Clculo del vector X:
for (i = n - 1; i >= 1; i--)
{
for (j = i + 1; j <= n; j ++)
B[i-1] -= CLectEscrMatrices.RecuSkyLine(j, i) * B[j-1] / A[idc.insk(i, i, Dir)];
}
return det;
}
393
// Datos miembro:
public static double[] vx; // Vector de abscisas de los puntos conocidos.
public static double[] vy; // Vector de ordenadas de los puntos conocidos.
public static double[] vdd; // Vector para calcular las diferencias divididas.
// metodos:
/* ASIGNA MEMORIA A LOS ARREGLOS */
public static void AsignaMemoria(bool t, int n)
{
vx = new double[n];
vy = new double[n];
if (t) vdd = new double[n];
}
394
double y, p = 1;
int i, j, k;
for (i = 0; i < n; i ++) vdd[i] = vy[i];
y = vy[0];
for (j = 0; j < n - 1; j ++) {
k = 0;
p *= x - vx[j];
for (i = 0; i < n - j - 1; i ++) {
vdd[i] = (vdd[i + 1] - vdd[i]) / (vx[k + j + 1] - vx[k]);
k ++;
}
y += vdd[0] * p;
}
return y;
395
DIAGRAMAS UML:
A fin de motivar al lector, en el diseo orientado a objetos, se presenta a continuacin
un diagrama de clases y otro de secuencia, mismos que fueron diseados para el
clculo de races reales, vistos en el captulo 1 de este libro. Estos diagramas ha sido
tomados del artculo Modelado Orientado a Objetos para evaluar Mtodos
Numricos utilizando Interfaces Visuales F. Meneses, W. Fuertes de la Revista
Tendencias en Computacin DECC Report.
396
CEntSalArchivos
-
CRaicesReales
+
+
+
+
+
+
+
+
+
+
+
+
a
b
xi
xd
dx
maxIter
precision
lb
punf
punfp
punfm
punfpm
:
:
:
:
:
:
:
:
:
:
:
:
double
double
double
double
double
int
double
ListBox
pf
pf
pf
pf
+
+
+
+
+
+
+
+
+
LocCambioSigno ()
Biseccion ()
AproxSucesivas ()
AproxSucesivasMod ()
NewtonRaphson ()
RegulaFalsi ()
Secante ()
FijaValores ()
MtdCalRaiz ()
+
+
+
+
+
+
+
+
Form : 1
NPath
NomArchEnt
NomArchSal
stLectura
stEscritura
:
:
:
:
:
<<Constructor>>
objes
:
:
:
:
:
:
:
:
:
double
double
double
double
double
double
double
void
int
String
String
String
StreamReader
StreamWriter
CEntSalArchivos (String p)
AbrirArchivoEnt ()
AbrirArchivoSal ()
LeerLinea ()
EscribirLinea (String s)
Escribir (String s)
CerrarLectura ()
CerrarEscritura ()
obent
objcalrai
CFuncionesT rascendentes
+
+
+
+
objes
objcalrai
objlb
objfunciones
:
:
:
:
ListBox
CEntSalArchivos
CRaicesReales
CListBox
double
b_lectura_click ()
b_salir_click ()
b_guardar_como_click ()
b_limpiar_clic ()
:
:
:
:
void
void
void
void
lb
CListBox
objlb
lb
: ListBox
+
+
+
+
Limpiar ()
GuardarComo (String spa)
Actualiza ()
Despliega ()
objfunciones
<<Delegate>>
DelegateDeclarations
objhor
pf (double x)
: double
+
+
+
+
objes
objhor
objcalrai
objlb
objvec
:
:
:
:
:
:
:
:
:
void
void
void
void
objlb
CFuncionesPolinomicas
<<Delegate>>
void
void
String
void
void
void
void
objes
objes
objcalrai
:
:
:
:
:
:
:
CEntSalArchivos
CHorner
CRaicesReales
CListBox
CLectEscrMatrices
b_lectura_click ()
b_salir_click ()
b_guardar_como_click ()
b_limpiar_clic ()
:
:
:
:
void
void
void
void
CLectEscrMatrices
objhor
CHorner
+
+
n
a
+
+
horner (double x)
fphorner (double x)
Form : 2
: int
: double[]
objvec
: double
: double
+
+
+
+
+
n
a
b
obent
lb
:
:
:
:
:
int
double[]
double[]
CEntSalArchivos
ListBox
+
+
LeeVector ()
EscribeVector (String ob)
397
: void
: void
user
CFuncionesTrascendentes
o
CFuncionesPolinomicas
Executa
Requerim.
Programa
CRaicesReales
CListBox
Captura datos
CRaicesReales.FijaValores()
CRaicesReales. MtdCalRaiz()
CListBox.Actualiza()
Siguiente
Mtodo
Despliega
Resultados
10 Requiere
Grabar
9
Siguiente Intervalo
8 CListBox.Despliega()
11 CListBox.GuardarComo()
12
CListBox.GuardarComo()
398
14
399
2.- Implementar un programa en Matlab y en Visual C# para ingresar un valor entre 1 y 1. De ser asi, el programa debe calcular el arco seno del valor ingresado,
aplicando la siguiente serie de Taylor:
sen -1(x) = x + (x3/3)(1/2) + (x5/5)(1*3)/(2*4) + (x7/7)(1*3*5)/(2*4*6) +
desplegando el resultado por pantalla e ingresando otro valor. El proceso se repite
hasta cuando se ingresa un valor fuera del rango indicado.
Para finalizar el clculo tomar hasta el trmino >= 1e-5 o se hayan alcanzado los 100
primeros trminos de la serie indicada.
400
BIBLIOGRAFIA:
1.- B.P. Demidovich, I.A. Maron. Clculo Numrico Fundamental. Paraninfo, ISBN:
84.283.0887 X, 1977.
2.- Smith, W Allen. Anlisis Numrico. Prentice-Hall, ISBN: 519.4 B896, 1998.
3.- Fausto Meneses, Walter Fuertes, Modelado Orientado a Objetos para evaluar
Mtodos Numricos utilizando Interfaces Visuales F. Meneses, W. Fuertes de la
Revista Tendencias en Computacin DECC Report, ISSN 1390 - 5236.
4.- Dr. Flix Calderon Solorio
http://lc.fie.umich.mx/~calderon/programacion/Mnumericos/Muller.html.
5.- www.umng.edu.co/www/resources/INTERPOLACION.doc.
6.- Cndido Pieiro Gmez "Apuntes de Matlab "
www.uhu.es/08003/aula_virtual/modulo_didactico/matlab.pdf.
7.- J. M. Gonzlez de Durana "Introduccin a MATLAB"
http://www.vc.ehu.es/depsi/jg/imatlab.pdf
8.- Pal Medina, Ph.D., "Mtodos Numricos, Derivacin e Integracin".
9.- JAAN KIUSALAAS. "Numerical Methods in Engineering with MATLAB".
10.- John H. Mathews, Kurtis D. Fink . "Numerical Methods using MATLAB".
11.- Frank Ayres Jr. Ph. D. "Ecuaciones Diferenciales TEORIA y 560 problemas
resueltos"
401
Publicaciones Cientficas