Instituto Politécnico Nacional Escuela Superior de Ingeniería Mecánica y Eléctrica Unidad Zacatenco Ingeniería en Comunicaciones y Electrónica
ASIGNATURA:
ANALISIS NUMERICOS
NOMBRE DEL TRABAJO:
TRABAJO DE INVESTIGACION SEGUNDO DEPARTAMENTAL
NOMBRE DEL PROFESOR:
ING. JOSE ANTONIO BERNAL MENDOZA
INTEGRANTES:
AVILA MACIAS JESUS EDUARDO GONZALEZ LOPEZ FERNANDO GUARNEROS DE LA TORRE HECTOR GUTIERREZ VICTOR RAUL MICHEL JULIAN ALEJANDRO
GRUPO: 4CV7
10/ NOVIEMBRE/ 2001
INDICE
INTRODUCCION ................................................................................................................................... 3 OBJETIVOS………………………………………………………………………………………………………………………………………4
MARCO TEORICO……………………………………………………………………………………………………………………………4
METODO METODO DE LU ............................................ ..................................................... ............ 5 PROBLEMAS ................................................................................................................................ 7 CODIGO FUENETE.................. ................. ................. .................. ................. ................. .............. 15
METODO METODO DE CHOLESKY CHOLESKY .......................................................................................... 16 CODIGO FUENTE........................................................................................................................ 22
METODO METODO DE CHEBYSHEV CHEBYSHEV.................................................... .................................... 26 PROBLEMAS…………………………………………………………………………………………………………………………..27
CODIGO FUENTE……………………………………………………………………………………………………………………..30 FUENTE……………………………………………………………………………………………………………………..30
METODO METODO DE HERMITE ................................................ ............................................. 32 PROBLEMAS .............................................................................................................................. 36
INTERPOLACION DE NEWTON…………………………………………………...38 PROBLEMAS……………………………………………………………………………………………………………………39.
CODIGO FUENTE………………………………………………………………………………………………………………41
CONCLUSIONES CONCLUSIONES ................................................................................................. .......... 60 GLOSARIO GLOSARIO .......................................................................................... ........................... 60 BOBLIOGRAFIA BOBLIOGRAFIA ...................................................... .................................................... . 62
Página 2
INDICE
INTRODUCCION ................................................................................................................................... 3 OBJETIVOS………………………………………………………………………………………………………………………………………4
MARCO TEORICO……………………………………………………………………………………………………………………………4
METODO METODO DE LU ............................................ ..................................................... ............ 5 PROBLEMAS ................................................................................................................................ 7 CODIGO FUENETE.................. ................. ................. .................. ................. ................. .............. 15
METODO METODO DE CHOLESKY CHOLESKY .......................................................................................... 16 CODIGO FUENTE........................................................................................................................ 22
METODO METODO DE CHEBYSHEV CHEBYSHEV.................................................... .................................... 26 PROBLEMAS…………………………………………………………………………………………………………………………..27
CODIGO FUENTE……………………………………………………………………………………………………………………..30 FUENTE……………………………………………………………………………………………………………………..30
METODO METODO DE HERMITE ................................................ ............................................. 32 PROBLEMAS .............................................................................................................................. 36
INTERPOLACION DE NEWTON…………………………………………………...38 PROBLEMAS……………………………………………………………………………………………………………………39.
CODIGO FUENTE………………………………………………………………………………………………………………41
CONCLUSIONES CONCLUSIONES ................................................................................................. .......... 60 GLOSARIO GLOSARIO .......................................................................................... ........................... 60 BOBLIOGRAFIA BOBLIOGRAFIA ...................................................... .................................................... . 62
Página 2
INTRODUCCIÓN
Muchas veces en el campo de la ingeniería, ciencia y vida diaria, se encuentran diversos problemas, los cuales deben ser analizados desde el enfoque que los métodos numéricos pueden dar, más aun, si dicho método numérico, es enfrentado a polinomios, de difícil trato. Sabemos que para resolver una ecuación, existen diferentes métodos, algunos de ellos ya los hemos utilizado en repetidas ocasiones, estos métodos son Método de sustitución
1 Se d espeja una incógnita en una u na de las ecuaciones. 2 Se sustituye la expresión de esta incógnita en la otra ecuación, obteniendo un ecuación con una sola incógnita. 3 Se resuelve la ecuación. 4 El valor obtenido se sustituye en la ecuación en la que aparecía la incógnita despejada. 5 Los dos valores obtenidos constituyen la solución del sistema. Método de igualación
1 Se despeja la misma incógnita en ambas ecuaciones. ecu aciones. 2 Se igualan las expresiones, con lo que obtenemos una ecuación con una incógnita. 3 Se resuelve la ecuación. 4 El valor obtenido se sustituye en cualquiera de las dos expresiones en las que aparecía despejada la otra incógnita. 5 Los dos valores obtenidos constituyen la solución del sistema. Método de reducción
1 Se preparan las dos ecuaciones, multiplicándolas por los números que convenga. 2 La restamos, y desaparece una de las incógnitas.
Página 3
3 Se resuelve la ecuación resultante. 4 El valor obtenido se sustituye en una de las ecuaciones iniciales y se resuelve. 5 Los dos valores obtenidos constituyen la solución del sistema. En este trabajo se llevo a cabo la investigación de nuevos métodos matemáticos, para poder resolver ecuaciones de cualquier tipo. Era lógico esperar un grado mayor de razonamiento, para poder aplicar estos métodos, ya que conforme avanzamos, encontraremos cada vez un nivel mayor de dificultad para resolver este tipo de problemas, sin embargo con los conocimientos previos que tenemos, será más que suficiente para poder aplicar estos métodos y no tener problema alguno. Recordamos que los métodos que se aplicaron en esta investigación son:
El método de LU Cholesky Chebyshev Hermite Y la interpolacion de Newton para 2 ecuaciones
Los cuales fueron detalladamente explicados de tal forma que fueran claros para nuestros lectores. OBJETIVOS
-GENERAL Este trabajo busca enseñar la implementación de distintos métodos para la resolución de sistemas de ecuaciones tales como Lu, Cholesky, Chebyshev y Hermite; además del estudio de métodos de interpolación, tales como el método de Newton. -PARTICULAR El objetivo particular de este trabajo tiene como principal meta, el lograr explicar en base a un marco teórico, desarrollo de método, problemas de ejemplo y una implementación en el compilador Turbo o C o DevC++, para cada uno de los métodos que más adelante se explicaran. MARCO TEORICO
Muchas veces, de una función sólo conocemos un conjunto de valores. Esto puede suceder, por ejemplo, porque son los resultados de un experimento gobernado por una ley que desconocemos. Si queremos calcular el valor de la función para una abscisa diferente de las conocidas, debemos utilizar otra función que la aproxime y, naturalmente, el valor que obtengamos será una aproximación del valor real. También puede suceder que sepamos la expresión analítica de la Página 4
función, pero sea lo suficientemente complicada como para calcular aproximaciones a los valores de la función a partir de otros ya conocidos. Existen varias formas de hacer esto, pero la más sencilla y una de las más utilizadas es la interpolación 1, que consiste en construir una función que pase por los valores conocidos (llamados polos) y utilizar ésta como aproximación de la función primitiva. Si se utilizan polinomios como funciones de aproximación, hablamos de interpolación polinómica. DESARROLLO METODO DE LU
Su nombre se deriva de las palabras inglesas ―Lower" y ―Upper‖, que en español se traducen como ―Inferior‖ y ―Superior‖. Estudiando el proceso que se
sigue en la descomposición LU es posible comprender el por qué de este nombre, analizando cómo una matriz original se descompone en dos matrices triangulares, una superior y otra inferior. La descomposición LU involucra solo operaciones sobre los coeficientes de la matriz [A], proporcionando un medio eficiente para calcular la matriz inversa o resolver sistemas de álgebra lineal. Primeramente se debe obtener la matriz [L] y la matriz [U]. [L] es una matriz diagonal inferior con números 1 sobre la diagonal. [U] es una matriz diagonal superior en la que sobre la diagonal no necesariamente tiene que haber números 1. El primer paso es descomponer o transformar [A] en [L] y [U], es decir obtener la matriz triangular inferior [L] y la matriz triangular superior [U]. Página 5
PASOS PARA ENCONTRAR LA MATRIZ TRIANGULAR SUPERIOR (MATRIZ [U] 1. Hacer cero todos los valores abajo del pivote sin convertir este en 1. 2. Para lograr lo anterior se requiere obtener un factor el cual es necesario para convertir a cero los valores abajo del pivote. 3. Dicho factor es igual al número que se desea convertir en cero entre el número pivote. 4. Este factor multiplicado por -1 se multiplica luego por el pivote y a ese resultado se le suma el valor que se encuentra en la posición a cambiar (el valor en la posición que se convertirá en cero). Esto es: - factor * pivote + posición a cambiar PASOS PARA ENCONTRAR LA MATRIZ TRIANGULAR NFERIOR (MATRIZ [L]) Para encontrar la matriz triangular inferior se busca hacer ceros los valores de arriba de cada pivote, así como también convertir en 1 cada pivote. Se utiliza el mismo concepto de ―factor‖ exp licado anteriormente y se ubican todos los ―factores‖ debajo de la diagonal según corresponda en cada uno.
Esquemáticamente se busca lo siguiente:
Originalmente se tenía:
Debido a que [A] = [L][U], al encontrar [L] y [U] a partir de [A] no se altera en nada la ecuación y se tiene lo siguiente
Página 6
Por lo tanto, si Ax = b, entonces LUx = b, de manera que Ax = LUx = b. PASOS PARA RESOLVER UN SISTEMA DE ECUACIONES POR EL MÉTODO DE DESCOMPOSICIÓN L 1. Obtener la matriz triangular inferior L y la matriz triangular superior U. 2. Resolver Ly = b (para encontrar y). 3. El resultado del paso anterior se guarda en una matriz nueva de nombre ―y‖.
4. Realizar Ux = y (para encontrar x). 5. El resultado del paso anterior se almacena en una matriz nueva llamada ―x‖, la cual brinda los valores correspondientes a las incógnitas de la ecuación.
PROBLEMAS RESUELTOS
1.- Encontrar los valores de x1, x2 y x3 para el siguiente sistema de ecuaciones:
NOTA: Recuérdese que si la matriz es 2x2 se hará 1 iteración; si es 3x3, 2 iteraciones; si es 4x4, 3 iteraciones; y así sucesivamente.
SOLUCIÓN: 4
-2
-1
Página 7
9
[A] = 5
1
-1
1
2
-4
[B] =
12
ITERACIÓN 1
factor 1 = (a21 / a11) = 5 / 4 = 1.25 factor 2 = (a31 / a11) = 1 / 4 = 0.25 Encontrando [U] fila 2 = - (factor 1) * (fila 1) + (fila 2) fila 3 = - (factor 2) * (fila 1) + (fila 3) a11 = a11 a12 = a12 a13 = a13 a21 = - (1.25) * (4) + (5) = 0 a22 = - (1.25) * (- 2) + (1) = 3.5 a23 = - (1.25) + (- 1) + (- 1) = 0.25 a31 = - (0.25) * (4) + (1) = 0 a32 = - (0.25) * (- 2) + (2) = 2.5 a33 = - (0.25) * (- 1) + (- 1) = - 0.75
[U] =
7
4
-2
-1
0
3.5
0.25
0
2.5
- 0.75
Página 8
Encontrando [L] 1
0
0
[L] = 1.25 0
0
0.25 0
0
ITERACIÓN 2
factor 3 = (u32 / u22) = 2.5 / 3.5 = 0.7142857143 Encontrando [U] fila 3 = - (factor 3) * (fila 2) + (fila 3) a31 = - (2.5 / 3.5) * (0) + (0) = 0 a32 = - (2.5 / 3.5) * (3.5) + (2.5) = 0 a33 = - (2.5 / 3.5) * (0.25) + (- 0.75) = - 0.9285714286
[U] =
4
-2
-1
0
3.5
0.25
0
0
- 0.9285714286
Encontrando [L]
1
0
0
[L] = 1.25 1
0
0.25 0.7142857143 1
Página 9
Ahora ya se tiene la matriz [U] y la matriz [L]. El siguiente paso es resolver Ly = b para encontrar la matriz y. En pocas palabras es como que se pidiera resolver el siguiente sistema de ecuaciones, encontrando los valores de y1, y2 y y3:
Al resolver el sistema anterior, se obtienen los siguientes valores para y 1, y2 y y3:
El último paso es resolver Ux = y para encontrar la matriz x. En otras palabras es como que se pidiera resolver el siguiente sistema de ecuaciones, encontrando los valores de x1, x2 y x3:
La solución del sistema es:
Este es finalmente el valor de x1, x2 y x3; es decir, la respuesta del ejercicio utilizando la descomposición LU.
Página 10
2.- Encontrar los valores de x1, x2 y x3 para el siguiente sistema de ecuaciones:
SOLUCIÓN:
[A] =
11
-3
-2
5
-2
-8
4
-7
2
ITERACIÓN 1 factor 1 = (a21 / a11) = 5/11 = 0.4545454545 factor 2 = (a31 / a11) = 4/11 = 0.3636363636 Encontrando [U] fila 2 = - (factor 1) * (fila 1) + (fila 2) fila 3 = - (factor 2) * (fila 1) + (fila 3) a11 = a11 a12 = a12 a13 = a13 a21 = - (0.4545454545) * (11) + (5) = 0 a22 = - (0.4545454545) * (- 3) + (- 2) = - 0.6363636365 a23 = - (0.4545454545) + (- 2) + (- 8) = - 7.0909090919 a31 = - (0.3636363636) * (11) + (4) = 0 Página 11
18 [B] =
13 2
a32 = - (0.3636363636) * (- 3) + (- 7) = - 5.909090909 a33 = - (0.3636363636) * (- 2) + (2) = 2.7272727272
[U] =
11
-3
-2
0
- 0.6363636365
- 7.0909090919
0
- 5.909090909
2.7272727272
Encontrando [L] 1 [L] =
0
0
0.45454545 0
0
0.36363636 0
0
ITERACIÓN 2
factor 3 = (u32/u22) = - 5.909090909 / - 0.6363636365 = 9.285714284 Encontrando [U] fila 3 = - (factor 3) * (fila 2) + (fila 3) a31 = - (9.285714284) * (0) + (0) = 0 a32 = - (9.285714284) * (- 0.6363636365) + (- 5.909090909) = 0 a33 = - (9.285714284) * (- 7.0909090919) + (2.7272727272) = 68.57142857
Página 12
[U] =
11
-3
-2
0
- 0.6363636365
- 7.0909090919
0
0
68.57142857
Encontrando [L] 1
0
0
[L] = 0.4545454545 1
0
0.3636363636 9.285714284 1
Ahora ya se tiene la matriz [U] y la matriz [L]. El siguiente paso es resolver Ly = b para encontrar la matriz y. En pocas palabras es como que se pidiera resolver el siguiente sistema de ecuaciones, encontrando los valores de y1, y2 y y3:
Al resolver el sistema anterior, se obtienen los siguientes valores para y 1, y2 y y3:
El último paso es resolver Ux = y para encontrar la matriz x. En otras palabras es como que se pidiera resolver el siguiente sistema de ecuaciones, encontrando los valores de x1, x2 y x3: Página 13
La solución del sistema es:
Este es finalmente el valor de x1, x2 y x3; es decir, la respuesta del ejercicio utilizando la descomposición LU.
Página 14
PROGRAMA DE LU CODIGO FUENTE
#include
#include #include //Librería gsl para calcular lo que queremos using namespace std; //Instrucción para no tener que escribir std::cin int main (int argc, char **argv) { size_t i,j,m; int s; cout<<"\nIntroduce el orden de la matriz: "; cin>>m; gsl_matrix * A = gsl_matrix_alloc(m,m); gsl_permutation * p = gsl_permutation_alloc(m); gsl_matrix *invA = gsl_matrix_alloc(m,m); cout<<"\nIntroduce los elementos de la matriz:"<
}
gsl_matrix_free(A); gsl_matrix_free(invA); gsl_permutation_free(p); return(0);
Página 15
METODO DE CHOLESKY
En álgebra lineal, la descomposición de Cholesky o triángulo de Cholesky es una descomposición de Hermite, definida positiva de la matriz en el producto de una matriz triangular inferior y su transpuesta conjugada. Fue descubierto por André-Louis Cholesky para matrices reales. Cuando sea aplicable, la descomposición de Cholesky es aproximadamente dos veces tan eficiente como la descomposición LU para resolver sistemas de ecuaciones lineales. En un sentido amplio, metafórico, esto puede considerarse como el análogo de la matriz de la raíz cuadrada de un número. Si A tiene entradas reales y es simétrica (o más generalmente, es hermitiano) y definida positiva, entonces A se puede descomponer como:
Donde L es una matriz triangular inferior con elementos de la diagonal estrictamente positiva, y L * denota la transpuesta conjugada de L . Esta es la descomposición de Cholesky. La descomposición de Cholesky es único: dado un hermitiana, definida positiva matriz A, sólo hay una matriz triangular inferior L con entradas diagonales estrictamente positivos tales que A = LL * .Lo contrario es trivial: si A se puede escribir como LL * para algunos invertible L, triangular inferior o de otra manera, entonces A es hermitiana y definida positiva. El requisito de que L tiene entradas de la diagonal estrictamente positiva se puede quitar para ampliar la factorización en el caso positivo semidefinida. La declaración se lee: una matriz cuadrada Atiene una descomposición de Cholesky si y sólo si A es hermitiana y positiva semidefinida. Factorizaciones de Cholesky para matrices semidefinida positiva no son únicas en general. En el caso especial de que A es simétrico definida positiva de la matriz con entradas reales, L tiene entradas reales también. La descomposición de Cholesky se utiliza principalmente para la solución numérica de ecuaciones lineales Ax = b. Si A es simétrica y definida positiva, entonces podemos resolver Ax = b calculando primero la descomposición de Cholesky A = LL T, entonces la solución de Ly = b para y, y, finalmente, la solución de L T x = y para x. Sistemas de la forma Ax = b con A simétrica y positiva definida surgen muy a menudo en las aplicaciones. Por ejemplo, las ecuaciones normales en al menos lineal plazas problemas son de esta forma. También puede suceder que la Página 16
matriz A proviene de un funcional de energía que debe ser positivo a partir de consideraciones físicas, lo que ocurre con frecuencia en la solución numérica de ecuaciones diferenciales parciales. No lineal multivariado funciones pueden ser minimizados por sus parámetros utilizando las variantes del método de Newton se llama cuasi-Newton métodos. En cada iteración, la búsqueda da un pasos definido por la resolución de Hs = -g de s, donde s es el paso, g es el gradiente del vector de derivadas parciales primera de la función con respecto a los parámetros, y H es una aproximación a la matriz hessiana de segundas derivadas parciales formados por repetidas rango 1 actualiza en cada iteración. Dos conocidas fórmulas de actualización se llama Davidon-Fletcher-Powell (DFP) y Broydon-Fletcher-GoldfarbShanno (BFGS). La pérdida de la condición definida positiva a través de errores de redondeo se evita si en lugar de actualizar una aproximación a la inversa de la de Hesse, una actualización de la descomposición de Cholesky de una aproximación de la matriz hessiana en sí. La descomposición de Cholesky se usa comúnmente en el método de Monte Carlo para la simulación de sistemas con múltiples variables correlacionadas: La matriz de covarianza se descompone para dar el triangular inferior L. Aplicando esto a un vector de muestras correlacionadas, u, produce un vector de muestras Lu con las propiedades de la covarianza del sistema que se modela. Filtros de Kalman suelen utilizar la descomposición de Cholesky para elegir un conjunto de los llamados puntos de sigma. El filtro de Kalman sigue el estado medio de un sistema como un vector x de longitud N y covarianza como una N por N de la matriz P. La matriz P es siempre positiva semidefinida, y se puede descomponer en LL T. Las columnas de L pueden sumarse y restarse de la media x para formar un conjunto de 2 N vectores llamados los puntos sigma. Estos puntos sigma totalmente la captura de la media y la covarianza del estado del sistema. El algoritmo de Cholesky, que se utiliza para calcular la matriz de descomposición L , es una versión modificada de la eliminación de Gauss . El algoritmo recursivo comienza con i : = 1 y
En el paso i , la matriz A ( i ) tiene la siguiente forma:
Página 17
donde I i -1 denota la matriz identidad de dimensión i - 1. Si ahora definimos la matriz L i por
entonces podemos escribir A ( i ) como
donde
Tenga en cuenta que b i b i * es un producto externo, por lo tanto, este algoritmo se llama la versión del producto exterior en (Golub y Van Loan). Repetimos esto para i desde 1 hasta n. Después de n pasos, se obtiene un ( n +1) = I . Por lo tanto, la matriz triangular inferior L que estamos buscando se calcula como:
El Cholesky-Banachiewicz y algoritmos de Cholesky-Crout
Si escribimos la ecuación A = LL *
se obtiene la siguiente fórmula para las entradas de L :
Página 18
La expresión en la raíz cuadrada es siempre positiva, si A es real y definida positiva. Para el complejo de la matriz hermitiana, la siguiente fórmula:
Por lo tanto podemos calcular la ( i , j entrada) si sabemos las entradas a la izquierda y arriba. El cálculo es generalmente dispuesto en cualquiera de las siguientes órdenes:
El algoritmo de Cholesky-Banachiewicz se inicia desde la esquina superior izquierda de la matriz de L y se procede a calcular la matriz fila por fila. El algoritmo de Cholesky-Crout se inicia desde la esquina superior izquierda de la matriz de L y se procede a calcular la matriz columna por columna.
Una forma alternativa es la factorización
Esta forma elimina la necesidad de tomar raíces cuadradas. Cuando A es definida positiva los elementos de la matriz diagonal D son todos positivos. Sin embargo, esta factorización se pueden utilizar para la mayoría de las matrices simétricas invertible; un ejemplo de una matriz invertible cuya descomposición no está definido es aquel en la primera entrada es cero. Si A es real, las relaciones recursivas se aplicarán las siguientes entradas de D y L: Página 19
Para el complejo hermitiano matriz A , la siguiente fórmula:
El LDL T y LL T factorizaciones (tenga en cuenta que L es la diferencia entre los dos) puede ser fácilmente relacionada con:
La última expresión es el producto de una matriz triangular inferior y su transposición, como es el LL T factorización. Los algoritmos anteriores muestran que cada matriz definida positiva A tiene una descomposición de Cholesky. Este resultado se puede extender al positivo semidefinido caso por un argumento de la limitación. El argumento no es totalmente constructivo, es decir, no da algoritmos numéricos explícitos para el cálculo de factores de Cholesky. Si A es un n -por- n positiva semi-definida de la matriz, entonces la sucesión { A k } = { A + (1 / k ) I n } consiste en matrices definidas positivas. (Esto es una consecuencia inmediata de, por ejemplo, el teorema de espectro para el cálculo funcional polinomio.) Además, en norma del operador. Desde el caso definido positivo, cada una k tiene descomposición de Cholesky A k = L k L k *. Pr la propiedad de la norma de operador, por lo tanto { L k } es un conjunto acotado en el espacio de Banach de los operadores, por lo tanto, relativamente compacto (porque el espacio vectorial subyacente es de dimensión finita). En consecuencia, tiene una subsucesión convergente, también se denota por { L k }, con un límite L . Se puede comprobar fácilmente que L tiene las propiedades deseadas, es decir, A = LL * y L es triangular inferior con no negativo elementos de la diagonal: para todos los x , y y , Página 20
Por lo tanto A = LL * . Debido a que el espacio vectorial subyacente es de dimensión finita, todas las topologías en el espacio de los operadores son equivalentes. Por lo tanto L k tiende a L en la norma significa que L k tiende a Lentrywise. Esto a su vez implica que, dado que cada L k es triangular inferior con no negativo elementos de la diagonal, L también. La factorización de Cholesky se puede generalizar a (no necesariamente finito) matrices con entradas del operador. Vamos a ser una secuencia de espacios de Hilbert. Considere la matriz del operador
que actúa sobre la suma directa
donde cada
es un operador acotado . Si A es positiva (semidefinida) en el sentido de que para todos los finitos k y para cualquier
tenemos , entonces existe un operador de la matriz triangular inferior L tal que A = LL * . Uno también puede tomar los elementos diagonales de L a ser positivo.
Página 21
CODIGO FUENTE PROGRAMA DEL METODO DE CHOLESKY
#include > #include #include #include #include #define TINY 1.0e-20 static int cholesky1(double **A, int n) { int i, j, k; for (i = 0; i < n; ++i) { for (j = 0; j <= i; ++j) { double sum = A[i][j]; A[j][i] = 0; for (k = 0; k < j; ++k) sum -= A[i][k] * A[j][k]; /* computed previously */ if (i == j) { if (sum < 0) return (0); sum = sqrt(sum); sqrt(sum); if (fabs(sum) < TINY) return (0); A[i][j] = sum; } else A[i][j] = sum / A[j][j]; } } return (1); } #undef TINY static void cholbksb(double **A, int n, double *x, double *b) { int i, j; double sum; for (i = 0; i < n; i++) { sum = b[i]; Página 22
for (j = 0; j < i; ++j) sum -= x[j] * A[i][j]; x[i] = sum / A[i][i];
}
} for (i = n - 1; i >= 0; i--) { sum = x[i]; for (j = i + 1; j < n; j++) sum -= x[j] * A[j][i]; x[i] = sum / A[i][i]; }
Matrix *matrix_cholesky_decomp(Matrix * A) { if (A == NULL || A->n != A->m) return (NULL); A = matrix_cast_fill(A, double_v); if (!cholesky1 (A->el.double_v, A->n)) { matrix_free((Matrix *) A); return (NULL); } return (A); } Vector *matrix_cholesky_back_sub(Matrix * A, Vector * b) { int n; Vector *x = NULL; Vector *c = NULL; if (A == NULL || A->n != A->m) return (NULL);
}
n = A->n; if (b == NULL || b->n != n) return (NULL); A = matrix_cast_fill(A, double_v); c = vector_cast(b, double_v); x = vector_alloc(n, double_v); cholbksb (A->el.double_v, n, (double *) x->data, (double *) c->data); vector_free(c); matrix_free((Matrix *) A); return (x);
Vector *matrix_cholesky_sol(Matrix * A, Vector * b) Página 23
{
int n; int sucess; Vector *x = NULL; Vector *c = NULL; if (A == NULL || A->n != A->m) return (NULL); n = A->n; if (b == NULL || b->n != n) return (NULL); A = matrix_cast_fill(A, double_v); sucess = cholesky1(A->el.double_v, n); if (sucess) { c = vector_cast(b, double_v); x = vector_alloc(n, double_v); cholbksb(A->el.double_v, n, (double *) x->data, (double *) c->data); vector_free(c); } matrix_free((Matrix *) A); return (x);
} Vector *matrix_cholesky_weighted_least_square(Matrix * A, Matrix * W, Vector * b) { Matrix *At; Matrix *M; Matrix *temp; Vector *x; if (A == NULL || A->m != b->n || W == NULL || W->n != W->m || W->m != A->m) return (NULL); A = matrix_cast_fill(A, matrix_cast_fill(A, double_v); /* to maintain precision */ At = matrix_transp(A); x = matrix_vprod(W, b); b = matrix_vprod(At, x); vector_free(x); temp = matrix_prod(W, A); M = matrix_prod(At, temp); matrix_free((Matrix *) temp); x = matrix_cholesky_sol(M, b); vector_free(b); matrix_free((Matrix *) A); matrix_free((Matrix *) At); matrix_free((Matrix *) M); return (x); } Vector *matrix_cholesky_least_square(Matrix * A, Vector * b) { Página 24
}
Matrix *At; Matrix *M; Vector *x; if (A == NULL || b == NULL || A->m != b->n) return (NULL); A = matrix_cast_fill(A, double_v); /* to keep precision */ At = matrix_transp(A); b = matrix_vprod(At, b); M = matrix_prod(At, A); x = matrix_cholesky_sol(M, b); vector_free(b); matrix_free((Matrix *) A); matrix_free((Matrix *) At); matrix_free((Matrix *) M); return (x);
Página 25
METODO DE CHEBYSHEV
Para cada valor de n , ¿cuál será la mejor forma de elegir los nodos? Aquella para la cual el valor máximo de |ω(x)| sea lo más pequeño posible: min {max {| (x –x 1)(x- x 2) … (x- x n+1) j : a ≤ x ≤ b} : a ≤ x 1 < x 2 <… < x n+1 ≤ b} Para resolver esta cuestión podemos restringirnos, haciendo una traslación y un cambio de escala si es necesario, al intervalo [ - 1, 1]. La solución viene dada mediante los polinomios de Chebyshev {T n+1(x ) }, donde T n+1(x ) = cos((n +1) arccos(x )). Recordemos que la funciónTn +1(x ) es la solución polinomica de la ecuación diferencial (1 - x 2)y” – xy´ + ( n + 1)2y = 0 cuyo coeficiente líder es 2 n . La forma más fácil de construirlos es mediante la relación de recurrencia: T 0(x ) = 1; T 1(x ) = x; T n+1(x ) = 2xT n (x ) - T n-1(x ), (n = 1, 2, … ).
Los polinomios de Chebyshev tienen muchas propiedades interesantes. La que a nosotros nos concierne es que el valor absoluto máximo de 2 -n T n+1(x ) para x є[- 1; 1] es 2-n y que, además, es el polinomio con menor máximo absoluto de entre todos los polinomios con coeficiente líder igual a 1; en otras palabras, que 2 -n T n+1(x ) es la solución al problema anterior o, equivalentemente, que los mejores nodos posibles son las n + 1 ra³ces simples de T n+1:
() ()
que se llaman nodos de Chebyshev del intervalo [ - 1, 1]. Como ya hemos indicado, para trabajar en un intervalo cualquiera hay que hacer un cambio de variable lineal que nos lleve el intervalo [ - 1, 1] a nuestro intervalo de trabajo [ a; b ]:
[] [] () entonces los nodos de Chebyshev en [a; b ] son los puntos z (x i) (i = 1, 2, … , n +1). Error de interpolación
( ) ( )( )( ) () ()
Nótese que según cómo se elijan los x i el error podrá ser mayor o menor
Página 26
PROBLEMAS RESUELTOS 1.-Dada la función
() en el intervalo [-5,5] Si construimos el polinomio de interpolación p n(x) en este intervalo, entonces seguro que no hay convergencia en los puntos donde |x|>3.63. Resulta que para ciertas funciones, por ejemplo para f(x)=exp(x), si hay convergencia. Pero para otras no. El problema con estas últimas es que sus derivadas van creciendo demasiado en el intervalo considerado, esto es lo que sucede con esta función de Runge, que parecía al principio bastante inofensiva. En las figuras que siguen mostramos este comportamiento:
Aproximamos la función por el polinomio de interpolación que pasa por los puntos en asterisco. Se puede observar que la aproximación es mala cerca de los extremos del intervalo Ya que la aproximación es mala en los extremos del intervalo, una idea para mejorar la aproximación es la de olvidarnos de tomar nodos igualmente Página 27
espaciados y tomar nodos que se concentren más cerca de los extremos. De este modo al obligar al polinomio p n(x) a pasar por estos puntos quizás se mejore la aproximación. Su interés en la teoría de la aproximación nació porque Chebyshev que era profesor de matemáticas en San Petersburgo, fue profesor allí durante 35 años, viajó en 1852 a Paris y a Inglaterra donde estuvo muy interesado en la máquina de vapor. Poco antes, James Watt había mejorado mucho el mecanismo de los pistones de la máquina. El mecanismo funciona principalmente porque el pistón, que tiene un movimiento de vaivén rectilíneo, se tiene que transformar luego en circular, porque la máquina al final tiene que mover una rueda. Esto se consigue por una serie de varillas articuladas ( linkages ) que convierten movimiento rectilíneo en circular. Sin embargo, a pesar de que Watt había conseguido una buena aproximación, se presentaban numerosos problemas. Esto era principalmente porque si miramos el movimiento al revés, es decir, de circular al girar una rueda que se transforme mediante una serie de brazos articulados en rectilíneo, el movimiento que resulta nos es completamente lineal, sino que presenta oscilaciones (curva de Watt ), ello hacía que se resintiera la máquina.
Ahora tomamos los nodos de Chebyshev y la aproximación se hace en el intervalo [-1 1] (para mayor comodidad ya que los ceros de Chebyshev se definen en este intervalo y no hay que hacer cambio de variable). Obsérvese que los nodos están más concentrados hacia los extremos del intervalo.
Página 28
Mecanismo de Peaucellier, está formado por una serie de brazos articulados, donde C está fijo y los demás son pivotes movibles. Cuando el punto F' describe una línea, el otro punto F traza una circunferencia. El polinomio de Chebyshev de grado n se denota Tn (x), y se le da la fórmula explícita Tn (x) = cos (n arc cos x) Sin embargo se puede combinar con las identidades trigonométricas para producir expresiones explicitas para Tn(x): A (x) = 1 T1 (x) = X T2 (X) = 2x2 - 1 (5.6.2) T3 (x) = 4x3 - 3x T4 (x) = 8x4 - 8x2 + 1
Página 29
CODIGO FUENTE PROGRAMA DEL METODO DE CHEBYSHEV
# Include # Include # Define NMAX 51 # Define MEDIA 0.5 # Define UN 1,0 # Define CUARTO 0,25 # Define DOS 2.0 # Define CERO 0.0 # Define PI 4.0 * atan (1.0) doble función (double x); CHEBFT void (doble A, doble B, double * C, int N) { suma doble, F [NMAX]; doble BMA, BPA, FAC, Y; Int J, K; BMA = Media * (B-A); BPA = Media * (B + A); for (K = 1, K <= N, K + +) { Y = cos (PI * (K-MEDIA) / N); F [K] = FUNC (BMA Y * + BPA); } FAC = DOS / N; for (j = 1, J <= N, J + +) { = SUMA CERO; if (K = 1, K <= N, K + +) SUMA + = F [K] * cos ((PI * (J-1 ))*(( K-MEDIA) / N)); C [J] = SUMA * FAC; } } CHEBEV doble (doble A, doble B, double * C, M int, double x) { doble D, DD, SV, Y, Y2; Int J; if ((XA) * (RE)> ZERO) printf ("\ n X, no en el rango \ n \ n".) D = CERO; DD = CERO; Y = (X * DOS-A-B) / (B-A) Y2 = DOS * Y; for (j = M, J> 1, J -) { Página 30
SV = D; D = Y2 * D-DD + C [J]; DD = SV;
} return (Y * D-DD + MEDIA * C [1]);
} CHINT void (doble A, doble B, double * C, el doble * CINT, int N) { CON doble, FAC, SUM; Int J; CON CUARTO * = (B-A); = SUMA CERO; FAC = UNO; for (j = 2; J = 3) for (j = N-2, J> 0, J -) CDER [J] = CDER [J 2] + DOS * J * C [J +1]; CON DOS = / (B-A); for (j = 1, J <= N, J + +) CDER [J] *= CON; }
Página 31
METODO DE HERMITE
Definimos los polinomios de Hermite por:
son polinomios de grado n. Se tiene que:
es decir, Hn es par si n es par, e impar si n es impar. Los primeros polinomios de Hermite son: H0 (t) = 1 H1 (t) = 2t H2 (t) = 4t2 - 2 H3 (t) = 8t3 -12t H4 (t) = 16t4 - 48t2 + 12 Consideremos la función
Su relación en serie de Taylor será:
Como:
Se tiene:
Luego Página 32
Se dice que e2tx-x2 es la función generatriz de los polinomios de Hermite, vale decir, es aquella función de dos variables tal que su desarrollo de Taylor en una de las variables tiene como coeficientes precisamente los polinomios de Hermite. A partir de la ecuación anterior se pueden encontrar relaciones entre los polinomios de Hermite. La estrategia para hallarlas (para _esta o cualquier otra función generatriz de otros polinomios) es típica: derivar parcialmente respecto a alguna de las variables y luego comparar potencias de x en los desarrollos en Taylor resultantes. 1) Derivando respecto a t
Reordenando la suma en el lado izquierdo
Comparando los coeficientes de las potencias de x en cada serie encontramos:
La cual puede ser reescrito en la forma
Página 33
Observemos que, si bien solo tiene sentido considerar polinomios de Hermite con índice positivo, la expresión anterior puede ser extendida a n = 0, aunque ello haga aparecer un factor H-1. En general, las relaciones de recurrencia que obtendremos pueden considerarse validas para cualquier índice entero, adoptando la convención de que los polinomios con subíndices negativos tienen algún valor adecuado, por ejemplo, cero. La relación anterior expresa un polinomio de Hermite en términos de un operador (en este caso la derivada) aplicado sobre el polinomio de Hermite inmediatamente superior. Un operador que tiene tal propiedad se denomina operador de bajada. En este caso, el operador de bajada de los polinomios de Hermite es (2n) -1dt 2) Derivando respecto a x
Comparando potencias de x:
O bien
3) Podemos utiliza las dos relaciones de recurrencia anteriores para obtener una tercera
Hemos pues encontrado el operador de subida para los polinomios de Hermite, a saber, 2t - dt. Página 34
Derivando la tercera relación:
Con:
O sea,
Es decir el, los polinomios H n son una solución de la ecuación de Hermite
Consideremos la ecuación de Hermite, pero generalicémosla ligeramente:
Busquemos soluciones con un cierto desarrollo de Taylor:
Remplazando en la ecuación de Hermite
es decir, obtenemos una relación de recurrencia para los coeficientes de la serie:
Se desprenden las siguientes consecuencias: Página 35
a) Hay dos series independientes, la de coeficientes con índice par, que depende sólo de a0, y la de coeficientes con índice impar, que depende sólo de a1. Por tanto, hay dos coeficientes arbitrarios, a 0 y a1, y por ende dos soluciones linealmente independientes
b) Lo cual significa que el radio de convergencia de la serie es infinito. Vale decir, las soluciones no tienen singularidades en el plano c) La ecuación tiene por solución un polinomio sólo si ᵝ є N*. Si ᵝ es par, hay que tomar a0 ≠ 0 y a 1 = 0. Si ᵝ es impar, hay que tomar a 1 ≠ 0 y a0 = 0. d) Si ᵝ є N*, y si la solución es par o impar, entonces desde cierto V0 en adelante, de modo que los a v tiene todos el mismo signo para V > V 0. Esto es, la serie tiene un crecimiento rápido cuando t→∞
PROBLEMAS RESUELTOS
1.- Utilizar el método de Hermite para hallar un polinomio P(x) de grado 2 que satisfaga: p(1) = 0, p´(1) = 7, p(2) = 10 Como existe la derivada del polinomio p(x), quiere decir que el método a utilizar es el de hermite X 1 1 2
P(x) 0 0 10
7 P[1,2] = 10
P[1,1,2]= 3
P[1,2] = 10 Por lo tanto P(x)= 7(x - 1) + 3(x - 1)2 P(x)= 3x2 + x – 4 que satisface las condiciones dadas 2.- Sea P(x)= 3 + ½(x – 1) + 1/3(x – 1)(x – 1.5) – 2x(x – 1)(x – 1.5)x, Q(x) = 5/3 – 2/3(x – 2) – 5/3(x – 2)(x – 1.5)(x – 2). Obtener las tablas de diferencias divididas que dan origen a P(x) y Q(x) respectivamente Página 36
Lo que nos indica que los polinomios son iguales, entonces si evaluamos un punto P(x), será igual que si lo evaluamos en Q(x) P(xi) = Q(xi), que pertenece i,i = {0,1,1.5,2} X 1 1.5 0 A
P(x) 3 P(1.5) P(0) P(a)
1/2 P[1.5, 0] P[0, a]
1/3 P[1.5, 0, a]
-2
X 2 0 1.5 b
Q(X) 5/3 Q(0) Q(1.5) Q(b)
-2/3 Q[0, 1.5] Q[1.5, b]
-5/3 Q[0, 1.5, b]
-2
Si vemos las tablas, la de P(x), esta tiene los puntos {1, 1.5, 0, a} y la de Q(x) tiene {2, 0, 1.5, b} b= 1; Q(b=1)= 3 y a= 2; P(a=2)= 5/3 PROBLEMAS PROPUESTOS
1.- Se sabe que H4(X)= 4 + 3(x + 1) – 2(x + 1)2 + 3/2(x + 1)2(x – 1) – ½(x – 1)2 (x – 1)2 es el polinomio de interpolación de Hermite de cierta función f, basado en los datos f(-1), f´(-1), f(1), f´(1) y f‖(1) Sin evaluar H4(x) ni sus derivadas en -1 y 1, completar la tabla de diferencias divididas con repetición utilizada en la construcción de H 4(x) 2.- Probar que si g interpola a la función f en x 0, x1, …, xn-1 y h interpola a f x0, x1, …, xn-1, xn en entonces a la función g(x) + (x0 – x/ xn – x0)(g(x) – h(x)) interpola a f x0, x1, … , xn-1, xn en notar que h y g no necesitan ser polinomios).
Página 37
METODO DE INTERPOLACION DE NEWTON
El método de Newton de diferencias divididas es otra forma de obtener el polinomio interpolador. En este método el polinomio interpolador se escribe de la forma:
Y el algoritmo proporciona una regla para obtener los coeficientes a0, a1,..., an. Imponiendo que el polinomio interpolador pase por los puntos de interpolación obtenemos:
De estas ecuaciones, es obvio que a0 depende sólo de x0, a1 de x0 y x1 y así sucesivamente. Introducimos la nueva notación:
Restando las dos primeras ecuaciones obtenemos:
Restando la segunda y la tercera ecuación obtenemos:
Podemos proceder de igual modo para demostrar que:
Aunque la forma más cómoda es por inducción construimos un polinomio de grado n definido por:
Relacionando nuestro polinomio interpolador, la ecuación estaría dada por:
Página 38
Relación que es el origen del nombre de diferencias divididas para los coeficientes an. Podemos por lo tanto escribir el polinomio interpolador como:
Y la formula de aproximación de f(x) queda en la siguiente forma:
f (x0), como se ve de la primera ecuación. Restando las dos primeras ecuaciones obtenemos PROBLEMAS RESUELTOS 1. Se desea interpolar los siguientes datos mediante el polinomio de Newton para x = 5.2.
Primero, calculamos la tabla de diferencias divididas:
Ahora evaluamos el polinomio de Newton:
Página 39
PROBLEMAS PROPUESTOS 1.-Da los siguientes datos de las temperaturas (en ºC ) Entre las 8.00 am a 8.00 pm. el 10 de mayo de 2005, en Kanpur:
Tiempo 8 a.m. Temperatura 30 37
12 horas 43 38
16:00 20:00
Obtener polinomio de interpolación de Newton hacia atrás de grado para calcular la temperatura en Kanpur en ese día a las 5.00 pm. 2.-La velocidad de un tren, que corre entre dos estaciones se mide a diferentes distancias de la estación de partida. Si x es la distancia en km de la estación de partida, a continuación v(x), la velocidad (en km/hr ) Del tren en la distancia x está dada por la siguiente tabla:
x v (x)
0 0
50 60
100 80
150 110
200 90
250 0
Encontrar la velocidad aproximada del tren en el punto medio entre las dos estaciones.
Página 40
CODIGO FUENTE PROGRAMA DE METODO DE NEWTON
# include # include # include # include # include # include const int max_size=13; int n=0; int top=-1; int choice=0; long double xn[max_size]={0}; long double px[max_size]={0}; long double dd_table[max_size][max_size]={0}; char Non_linear_equation[100]={NULL}; char Stack[30][30]={NULL}; char Postfix_expression[30][30]={NULL}; void push(const char *); void convert_infix_expression_to_postfix_expression(const char *); const char* pop( ); const long double evaluate_postfix_expression(const long double); void show_screen( ); void clear_screen( ); void get_input( ); void construct_dd_table( ); void show_dd_table( ); void compute_newtons_dd_interpolation_polynomial( ); void estimate_pnx( ); int main( ) { clrscr( ); textmode(C4350); show_screen( ); get_input( ); construct_dd_table( ); show_dd_table( ); compute_newtons_dd_interpolation_polynomial( ); estimate_pnx( ); return 0; } Página 41
void show_screen( ) { textbackground(1); cprintf(" Newton's Divided Difference Interpolation Formula "); textbackground(8); for(int count=0;count<42;count++) cprintf("* *"); gotoxy(1,46); cprintf("********************************************************************************"); gotoxy(1,2); } void clear_screen( ) { for(int count=0;count<37;count++) { gotoxy(3,8+count); cout<<" } gotoxy(1,2); } void push(const char* Operand) { if(top==(max_size-1)) { cout<<"Error : Stack is full."<
";
} else {
}
getch( ); exit(0);
strcpy(Operand,Stack[top]); strset(Stack[top],NULL); top--;
} return Operand;
void convert_infix_expression_to_postfix_expression(const char* Expression) { char Infix_expression[100]={NULL}; char Symbol_scanned[30]={NULL}; push("("); strcpy(Infix_expression,Expression); strcat(Infix_expression,"+0)"); int flag=0; int count_1=0; int count_2=0; int equation_length=strlen(Infix_expression); if(Infix_expression[0]=='(') flag=1; do { strset(Symbol_scanned,NULL); if(flag==0) { int count_3=0; do { Symbol_scanned[count_3]=Infix_expression[count_1]; count_1++; count_3++; } while(count_1<=equation_length && Infix_expression[count_1]!='(' && Infix_expression[count_1]!='+' && Infix_expression[count_1]!='-' && Infix_expression[count_1]!='*' && Infix_expression[count_1]!='/' && Infix_expression[count_1]!='^' && Infix_expression[count_1]!=')'); flag=1; } Página 43
else if(flag==1) { Symbol_scanned[0]=Infix_expression[count_1]; count_1++; if(Infix_expression[count_1]!='(' && Infix_expression[count_1]!='^' && Infix_expression[count_1]!='*' && Infix_expression[count_1]!='/' && Infix_expression[count_1]!='+' && Infix_expression[count_1]!='-' && Infix_expression[count_1]!=')') flag=0; if(Infix_expression[count_1-1]=='(' && (Infix_expression[count_1]=='-' || Infix_expression[count_1]=='+')) flag=0; } if(strcmp(Symbol_scanned,"(")==0) push("("); else if(strcmp(Symbol_scanned,")")==0) { while(strcmp(Stack[top],"(")!=0) { strcpy(Postfix_expression[count_2],pop( )); count_2++; } pop( ); } else if(strcmp(Symbol_scanned,"^")==0 || strcmp(Symbol_scanned,"+")==0 || strcmp(Symbol_scanned,"-")==0 || strcmp(Symbol_scanned,"*")==0 || strcmp(Symbol_scanned,"/")==0) { if(strcmp(Symbol_scanned,"^")==0) { } else if(strcmp(Symbol_scanned,"*")==0 || strcmp(Symbol_scanned,"/")==0) { while(strcmp(Stack[top],"^")==0 || strcmp(Stack[top],"*")==0 || strcmp(Stack[top],"/")==0) { strcpy(Postfix_expression[count_2],pop( )); Página 44
}
}
count_2++;
else if(strcmp(Symbol_scanned,"+")==0 || strcmp(Symbol_scanned,"-")==0) { while(strcmp(Stack[top],"(")!=0) { strcpy(Postfix_expression[count_2],pop( ));
} } else {
}
}
count_2++;
push(Symbol_scanned);
strcat(Postfix_expression[count_2],Symbol_scanned); count_2++;
} while(strcmp(Stack[top],NULL)!=0);
}
strcat(Postfix_expression[count_2],"="); count_2++;
const long double evaluate_postfix_expression(const long double x) { long double function_value=0; int count_1=-1; char Symbol_scanned[30]={NULL}; do
{
count_1++; strcpy(Symbol_scanned,Postfix_expression[count_1]); if(strcmp(Symbol_scanned,"^")==0 || strcmp(Symbol_scanned,"*")==0 || Página 45
strcmp(Symbol_scanned,"/")==0 || strcmp(Symbol_scanned,"+")==0 || strcmp(Symbol_scanned,"-")==0) {
char Result[30]={NULL}; char Operand[2][30]={NULL}; strcpy(Operand[0],pop( )); strcpy(Operand[1],pop( )); long double operand[2]={0}; long double result=0; char *endptr; for(int count_2=0;count_2<2;count_2++) { int flag=0; if(Operand[count_2][0]=='-') { int length=strlen(Operand[count_2]); for(int count_3=0;count_3<(length-1);count_3++)
Operand[count_2][count_3]=Operand[count_2][(count_3+1)]; Operand[count_2][count_3]=NULL; }
flag=1;
if(strcmp(Operand[count_2],"x")==0) operand[count_2]=x; else if(strcmp(Operand[count_2],"e")==0) operand[count_2]=2.718282; else if(strcmp(Operand[count_2],"sinx")==0) operand[count_2]=sinl(x); else if(strcmp(Operand[count_2],"cosx")==0) operand[count_2]=cosl(x); else if(strcmp(Operand[count_2],"tanx")==0) operand[count_2]=tanl(x); Página 46
else if(strcmp(Operand[count_2],"lnx")==0) operand[count_2]=logl(x); else if(strcmp(Operand[count_2],"logx")==0) operand[count_2]=log10l(x); else operand[count_2]=strtod(Operand[count_2],&endptr);
}
if(flag) operand[count_2]*=-1;
switch(Symbol_scanned[0]) { case '^' : result=powl(operand[1],operand[0]); break; case '*' : result=operand[1]*operand[0]; break; case '/' : result=operand[1]/operand[0]; break; case '+' : result=operand[1]+operand[0]; break;
}
case '-' : result=operand[1]-operand[0]; break;
gcvt(result,25,Result); }
push(Result);
else if(strcmp(Symbol_scanned,"=")!=0) push(Symbol_scanned);
} while(strcmp(Symbol_scanned,"=")!=0); char Function_value[30]={NULL}; char *endptr;
strcpy(Function_value,pop( )); function_value=strtod(Function_value,&endptr); Página 47
}
return function_value;
void get_input( ) { do { clear_screen( ); gotoxy(6,9); cout<<"Number of Distinct Data Points :"; gotoxy(6,10); cout<<"ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ"; gotoxy(27,13); cout<<"[ min. n = 2 | max. n = 12 ]"; gotoxy(6,12); cout<<"Enter the max. number of distinct data points = n = "; cin>>n; if(n<2 || n>12) { gotoxy(12,25); cout<<"Error : Wrong Input. Press to exit or any other key"; gotoxy(12,26); cout<<" to try again."; n=int(getche( ));
}
if(n==27) exit(0);
} while(n<2 || n>12); gotoxy(6,19); cout<<"Input Mode :"; gotoxy(6,20); cout<<"ÍÍÍÍÍÍÍÍÍÍÍÍ"; gotoxy(8,23); Página 48
cout<<"Press : "; gotoxy(10,25); cout<<"- 'Y' or to enter function"; gotoxy(10,27); cout<<"- 'N' or to enter values of the function"; gotoxy(8,30); cout<<"Enter your choice : "; char Choice=NULL; Choice=getch( ); if(Choice=='y' || Choice=='Y' || int(Choice)==13) { choice=1;
} else { }
gotoxy(28,30); cout<<"Y";
gotoxy(28,30); cout<<"N";
gotoxy(25,43); cout<<"Press any key to continue..."; getch( ); if(choice) { clear_screen( ); gotoxy(6,11); cout<<"Non-Linear Function :"; gotoxy(6,12); cout<<"ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ";
(2*x)+3";
gotoxy(6,37); cout<<"Note : Write the function with proper Braces ( ) e.g; 2x+3 as Página 49
gotoxy(6,40); cout<<"Available Operators : ^ (raised to power) , * , / , + , -"; gotoxy(6,42); cout<<"Available Operands : x , e , sinx , cosx , tanx , lnx , logx ,"; gotoxy(6,44); cout<<"
n = any number";
gotoxy(6,14); cout<<"Enter the Function : f(x) = "; cin>>Non_linear_equation; }
convert_infix_expression_to_postfix_expression(Non_linear_equation);
clear_screen( ); gotoxy(6,9); cout<<"Data Points & Values of Function :"; gotoxy(6,10); cout<<"ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ"; gotoxy(25,12); cout<<"ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÂÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿"; gotoxy(25,13); cout<<"³ x
³
f(x)
³";
gotoxy(25,14); cout<<"ÃÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÅÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ´"; gotoxy(25,15); cout<<"³
³
³";
for(int count_1=0;count_1
}
gotoxy(25,(wherey( )+1)); cout<<"³ ³ ³";
Página 50
gotoxy(25,(wherey( )+1)); cout<<"ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÁÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ"; gotoxy(25,16); for(int count_2=0;count_2>xn[count_2]; if(choice) { dd_table[0][count_2]=evaluate_postfix_expression(xn[count_2]);
} else { }
gotoxy(43,(wherey( )-1)); cout<
gotoxy(43,(wherey( )-1)); cin>>dd_table[0][count_2];
if(choice) gotoxy(25,(wherey( )+2));
}
else gotoxy(25,(wherey( )+1));
gotoxy(25,43); cout<<"Press any key to continue..."; }
getch( );
void construct_dd_table( ) { for(int count_1=1;count_1
long double x2=xn[count_2];
}
}
}
dd_table[count_1][count_2]=((fx2-fx1)/(x2-x1));
void show_dd_table( ) { clear_screen( ); gotoxy(4,10); cout<<"Divided Difference Table :"; gotoxy(4,11); cout<<"ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ"; gotoxy(4,13); cout<<"ÚÄÄÄÄÄÄÄÄÂÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ"; gotoxy(4,14); cout<<"³ x ³
f(x)
";
gotoxy(4,15); cout<<"ÃÄÄÄÄÄÄÄÄÅÄÄÄÄÄÄÄÄÄÄÄÄÄÄ"; gotoxy(4,16); cout<<"³ ³
";
int x_cord=4; int y_cord=17; for(int count_1=0;count_1
}
y_cord+=2;
gotoxy(x_cord,y_cord); cout<<"ÀÄÄÄÄÄÄÄÄÁÄÄÄÄÄÄÄÄÄÄÄÄÄÄ"; x_cord=28; int count_2=0; for(int count_3=1;count_3
cout<<"³"; gotoxy(x_cord,(y_cord+1)); cout<<"³"; if(count_4>=count_3) { gotoxy((x_cord+2),(y_cord-count_3)); cout<
y_cord+=2;
x_cord+=16; count_2++; if((count_2%3)==0 && count_3<(n-1)) { gotoxy(x_cord,13); cout<<"ÂÄÄ"; gotoxy(x_cord,14); cout<<"³"; gotoxy(x_cord,15); cout<<"ÅÄÄ"; gotoxy(x_cord,16); cout<<"³"; gotoxy(x_cord,y_cord); cout<<"ÁÄÄ"; y_cord=17; for(int count_5=0;count_5
y_cord+=2;
Página 54
gotoxy(30,44); cout<<"Press any key to continue..."; getch( ); y_cord=13; x_cord=28; for(int count_6=0;count_6<=(n+2);count_6++) { gotoxy(x_cord,y_cord); cout<<" "; gotoxy(x_cord,(y_cord+1)); cout<<" }
}
}
y_cord+=2;
y_cord-=2; count_2=0; count_3--;
gotoxy(x_cord,13); cout<<"¿"; gotoxy(x_cord,14); cout<<"³"; gotoxy(x_cord,16); cout<<"³"; y_cord=17; for(int count_6=0;count_6
y_cord+=2;
Página 55
";
gotoxy(x_cord,15); cout<<"´"; gotoxy(x_cord,y_cord); cout<<"Ù"; gotoxy(30,44); cout<<"Press any key to continue..."; }
getch( );
void compute_newtons_dd_interpolation_polynomial( ) { for(int count_1=0;count_10;count_4--) { fx[count_4]*=xn[(count_3+1)]; fx[count_4]+=fx[(count_4-1)]; } }
fx[0]*=xn[(count_3+1)];
for(int count_5=0;count_5<(count_2+1);count_5++) fx[count_5]*=dd_table[count_2][count_2];
}
}
for(int count_6=0;count_6<(count_2+1);count_6++) px[count_6]+=fx[count_6];
Página 56
void estimate_pnx( ) { clear_screen( ); gotoxy(6,10); cout<<"Newton's Divided Difference Interpolation Formula :"; gotoxy(6,11); cout<<"ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ"; gotoxy(10,13); cout<<"P"<<(n-1)<<"(x) = "; for(int count_1=(n-1);count_1>=0;count_1--) { if(count_1==(n-1) && n>2) cout<1) cout<=70) gotoxy(30,(wherey( )+2)); if(count_1>0) { if(px[(count_1-1)]>0) cout<<" + ";
}
}
else cout<<" - ";
if(wherex( )>=70) gotoxy(30,(wherey( )+2));
gotoxy(6,19); cout<<"Estimation of Pn(x) :"; Página 57
gotoxy(6,20); cout<<"ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ"; char Choice=NULL; long double x=0; long double pnx=0; do
{
Choice=NULL; x=0; pnx=0; gotoxy(10,22); cout<<"Enter the value of x = "; cin>>x; pnx=px[0]; pnx+=(px[1]*x); for(int count_2=2;count_2
pnx+=(px[count_2]*temp);
gotoxy(10,26); cout<<"The estimated value of P"<<(n-1)<<"("<
} gotoxy(15,40); cout<<"Press to exit, 'V' to view D.D. Table again"; gotoxy(25,42); cout<<"or any other key to continue..."; Choice=getch( ); if(int(Choice)!=27 && Choice!='v' && Choice!='V') { gotoxy(10,22); cout<<" ";
}
gotoxy(10,26); cout<<"
";
gotoxy(10,28); cout<<"
";
gotoxy(10,30); cout<<"
";
gotoxy(15,40); cout<<"
";
gotoxy(25,42); cout<<"
";
else if(Choice=='v' || Choice=='V') { show_dd_table( ); estimate_pnx( ); } else if(int(Choice)==27) exit(0);
}
} while(1);
Página 59
CONCLUSIONES
Luego de haber estudiado a profundidad este tema o herramienta para resolver sistemas de ecuaciones, se concluye que para resolver estos sistemas de ecuaciones lineales existen diferentes métodos, pero dependerá del gusto de cada persona elegir uno en específico. Sin embargo, muchas veces la elección no será arbitraria, pues cada método tiene sus ventajas y sus desventajas. Algunos métodos son más exactos, otros más fáciles de programar, otros más cortos, etc. Para ser capaces de elegir un método apropiado, lo primero que se necesita es comprender cómo se desarrolla cada uno de estos procesos. El método de interpolación , es una herramienta muy utilizada en el ámbito laboral, debido a que gracias a sus procesos de aproximación, en muchos casos se pueden restablecer datos perdidos u extraviados de diversas bases de datos, tales como el control de pedidos en un almacén, o los ingresos y egresos que una sucursal de una panadería puede tener. Luego de la elaboración de este reporte, ya se tiene una buena base y el conocimiento del tema para poder comenzar a programar en la computadora estos procesos. Como se mencionó en la introducción, el método estudiado en este trabajo es ideal para programarlos por computadora, pues son iterativos y muy largos. Trabajar esto en papel podría resultar extremadamente largo y tedioso. Por ello son métodos ideales para trabajarlos en computadora. GLOSARIO DE TERMINOS Interpolación: En el subcampo matemático del análisis numérico, se denomina interpolación a la obtención de nuevos puntos partiendo del conocimiento de un conjunto discreto de puntos.
En ingeniería y algunas ciencias es frecuente disponer de un cierto número de puntos obtenidos por muestreo o a partir de un experimento y pretender construir una función que los ajuste. Interpolación: En el subcampo matemático del análisis numérico, se denomina interpolación a la obtención de nuevos puntos partiendo del conocimiento de un conjunto discreto de puntos. Matriz: En matemáticas, una matriz es una tabla bidimensional de números consistente en cantidades abstractas que pueden sumarse y multiplicarse. Las matrices se utilizan para describir sistemas de ecuaciones lineales, realizar un seguimiento de los coeficientes de una aplicación lineal y registrar los datos que Página 60
dependen de varios parámetros. Las matrices se describen en el campo de la teoría de matrices. Pueden sumarse, multiplicarse y descomponerse de varias formas, lo que también las hace un concepto clave en el campo del álgebra lineal. Ecuación: Una ecuación es una igualdad matemática entre dos expresiones algebraicas, denominadas miembros , en las que aparecen valores conocidos o datos , y desconocidos o incógnitas , relacionados mediante operaciones matemáticas. Los valores conocidos pueden ser números, coeficientes o constantes; y también variables cuya magnitud se haya establecido como resultado de otras operaciones. Las incógnitas, representadas generalmente por letras, constituyen los valores que se pretende hallar. Polinomio: En matemáticas , se le llama polinomio a la suma de varios monomios. Un monomio es un producto de un coeficiente y una variable elevado a un número natural, que se llama el exponente del monomio. Ortogonal: En matemáticas, el término ortogonalidad (del griego orthos —recto— y gonía —ángulo—) es una generalización de la noción geométrica de perpendicularidad. Cifra Significativa: Todas las cifras excepto el cero Ecuación cuadrática: Ecuación de segundo grado o cuadrática se expresa mediante la relación ax2 + bx + c = 0, donde a es distinto de 0. Ecuación cúbica: Las ecuaciones de tercer grado o cúbicas son del tipo ax3 + bx2 + cx +d = 0, donde a es distinto de 0. Ecuación cuártica: Las ecuaciones de cuarto grado o cuárticas, ax4 + bx3 + cx2 + dx + e = 0, para a distinto de 0. Ecuación Diferencial: Ecuación que contiene derivadas. Evaluar: Valorar una cosa
Página 61
