Simpson & Romberg

Por: Arturo Fornés A01227071 y Miguel Montoya A01226045

Nota: Tanto en el diagrama de flujo, como en el código y los ejemplos, utilizaremos el método de Simpson 1/3 múltiple.

Antecedentes

Hay varios métodos para calcular integrales definidas donde se hace una aproximación por cada intervalo, como Riemann desde la izquierda y derecha, o por método de trapecios. Existen dos métodos más complejos que obtienen una mejor aproximación, estos son los métodos de Simpson (Trataremos especialmente el de 1/3 múltiple) y el e Romberg.

Simpson 1/3 está basado en la “interpolación cuadrática”, pues crea una función que pasa por todos los puntos dados e integra. La función es para dos intervalos continuos (Uniendo 3 puntos). La aplicación múltiple, es crear pares de intervalos para aplicar Simpson 1/3 a muchos puntos.

Romberg está basado en la integración por trapecios, y utiliza un método parecido a las diferencia divididas de Newton. Obtiene aproximaciones de integración con diferente número de trapecios para cada aproximación, y las mejora, hasta obtener un mejor resultado para la integral.

Justificación y propósito

Simpson 1/3 utiliza la siguiente función, para calcular la función perteneciente a los intervalos, antes de integrar.

{\displaystyle P_{2}(x)=f(a){\frac {(x-m)(x-b)}{(a-m)(a-b)}}+f(m){\frac {(x-a)(x-b)}{(m-a)(m-b)}}+f(b){\frac {(x-a)(x-m)}{(b-a)(b-m)}}.}

Debido a que siempre es la misma formula para obtener el polinomio, es fácil calcular la integral.

{\displaystyle I=\int _{a}^{b}f(x)\,dx}

Donde f(x) es el polinomio de Lagrange de grado dos. Podemos calcular una ecuación para cada intervalo, sin necesidad de calcular la integral de cada uno.{\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{6}}\left[f(a)+4f(m)+f(b)\right].}

Por otro lado, Romberg mejora el método de trapecios utilizando la siguiente formula

Screenshot from 2017-05-04 13-21-56

Donde se obtiene una integral mejorada, utilizando otras dos aproximaciones de trapecios, de menor nivel (Uno más exacto que el otro). Se denomina como aproximación más exacta a la que utilizo un mayor número de trapecios. Este método es recursivo, si obtienes x aproximaciones de un nivel mayor, todavía puedes obtener x-1 aproximaciones de nivel mayor.

Diagrama de flujo

Simpson 1/3 Múltiple.

Simpson1-3Multiple

Romberg.

Romberg(1)

Código en C++

Simpson 1/3 Múltiple: Github.

Romberg: Github.

Ejemplo comparativo

Para Simpson 1/3 múltiple, no se utiliza el método de trapecio para intervalos extra (En caso de que el número de intervalos sea impar, el último intervalo se integra con trapecios).


f (x) = x²

Se integra f(x) de 1 a 5:s_x2_integral

rs_x2_graph

Usando Simpson de 1/3 múltiple se obtiene una aproximación con una precisión de 0.00001:

s_x2_code

Usando Romberg de 3 aproximaciones iniciales, usando aproximaciones iniciales de 100, 150 y 200 trapecios:

r_x2_code

Romberg no fue preciso en este caso ya que este método fluctúa dependiendo de los valores de las aproximaciones por trapecios que se obtengan y qué tan diferentes sean estas entre sí.


f(x) = e^(x)

Se integra f(x) de 1 a 3:

rs_ex_int

rs_ex_graph

Usando Simpson de 1/3 múltiple se obtuvo una aproximación con una precisión de 0.1:

s_ex_code

Usando Romberg de 3 aproximaciones iniciales, usando aproximaciones iniciales de 100, 150 y 200 trapecios se obtuvo una aproximación a la integral definida con una precisión de 0.1:

r_ex_code

En este caso ambos métodos dieron valores certeros, con el mismo grado de precisión.

Anuncios

Interpolación de Newton y Lagrange

Antecedentes

Ambas interpolaciones son métodos que permiten la creación de un polinomio de grado n-1, donde n es el número de datos que se tienen. Ambos métodos asumen que no existe ruido en sus mediciones de datos, es decir, que el polinomio generado por el método pasara por todas las coordenadas insertadas al método.

A menos que se conozca la función original (Lo cual es prácticamente imposible al aplicarlos en la vida real), no hay manera de calcular algún error.

Justificación y propósito

El propósito de ambos métodos es poder generar una función para la cual se puedan introducir todos los datos originales y obtener 0 error, pues la curva se va modelando punto a punto. Al obtener una función, se puede crear aproximaciones y estimaciones.
(En nuestros métodos, no proporcionamos al usuario la función, solamente el valor f(x) para la x deseada)

El método de Newton utiliza la formula de polinomios de Newton:

N(x)=[y_{0}]+[y_{0},y_{1}](x-x_{0})+\cdots +[y_{0},\ldots ,y_{k}](x-x_{0})(x-x_{1})\cdots (x-x_{{k-1}}).
Donde las “y” entre corchetes se refieren al calculo de diferencias dividas, el cual es un algoritmo usado para computar tablas de funciones logarítmicas y trigonométricas, usando división recursiva.

El método de Larange obtiene una misma función, pero elimina por completo la necesidad de usar diferencias divididas:

L(x):=\sum _{j=0}^{k}y_{j}\ell _{j}(x)

donde

\ell _{j}(x):=\prod _{\begin{smallmatrix}0\leq m\leq k\\m\neq j\end{smallmatrix}}{\frac {x-x_{m}}{x_{j}-x_{m}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{j-1})}{(x_{j}-x_{j-1})}}{\frac {(x-x_{j+1})}{(x_{j}-x_{j+1})}}\cdots {\frac {(x-x_{k})}{(x_{j}-x_{k})}},

Los polinomios de lagrange son más fáciles de computar, pues elimina la necesidad de recurrir a métodos de recursión,

Explicación gráfica de la interpolación

 

Estas gráficas representan dado un conjunto de puntos la interpolación de grado 1 a 4 de una funciócon la intención de obtener el valor de f(10), y cómo se ve reflejado esto en una gráfica.

Para la gráfica a) se hace una interpolación de grado 4, por lo que se utilizan los 5 puntos dados y se obtiene una función que cruce por estos puntos, luego se evalúa en x = 10.

Para la gráfica b) se hace una interpolación cúbica, por lo que sólo se utilizan los 4 puntos más cercanos a 10, se obtiene una función que cruce por estos puntos y se evalúa en x = 10.

Para la gráfica c) se hace una interpolación cuadrática, por lo que sólo se utilizan los 3 puntos más cercanos a 10, se obtiene una función cuadrática que cruce por estos puntos y se evalúa esa función en x = 10.

Para la gráfica d) se hace una interpolación lineal, por lo que sólo se utilizan los 2 puntos más cercanos a 10, se obtiene una función lineal que cruce por estos puntos y se evalúa esa función en x = 10.

Por estas gráficas se puede observar, también, la influencia que tiene el tamaño de la muestra para el resultado.

Diagrama de flujo

Newton

NewtonInterpolation

Lagrange

Lagrange Intepolation

Código en C++

Ejemplos resueltos

A continuación se muestra los resultados de varias muestras resueltas por ambos métodos:

Con los siguientes puntos:

x, f(x) = (1, 0), (4, 1.386294), (6, 1.791759)

Se interpola usando estos datos para obtener una aproximación de f(5). Estos datos son muestras de la función f(x) = ln(x), por lo que ya conocemos el valor de f(5) = 1.61.

int_graph1

Con Polinomios de Interpolación de Newton se obtienen los valores para b0 … bn y aproximación siguientes:

int_1newton

Para Polinomios de Interpolación de Lagrange se obtienen los coeficientes de Lagrange y aproximación siguientes:

int_1lagrange


Con los siguientes puntos:

x, f(x) = (1, 1), (2, 0.5), (2.5, 0.4), (4, 0.25), (7, 0.14)

Se interpola usando estos datos para obtener una aproximación de f(3). Estos datos son muestras de la función f(x) = 1/x, por lo que ya conocemos el valor de f(3) = 0.33.

int_2grapg

Con Polinomios de Interpolación de Newton se obtienen los valores para b0 … bn y aproximación siguientes:

int_2newton

Para Polinomios de Interpolación de Lagrange se obtienen los coeficientes de Lagrange y aproximación siguientes:

int_2lagrange

Regresión Polinomial

Antecedentes

Como ya fue explicado en el anterior post. La regresión (De cualquier tipo) es utilizada para generar una recta o curva que se ajusta a un diagrama de dispersión. Este diagrama contiene datos previamente capturados.

En el caso de regresión polinomial, el objetivo es crear una función que mejor se adapte a los datos introducidos. El grado de la función es especificado por el usuario.

Justificación y propósito

El propósito de todas las regresiones lineales es encontrar la relación de la variable dependiente con la independiente. Así realizar hipótesis, analizar tendencias y hacer estimaciones. En este caso se introducen términos polinomiales para realizar esta aproximación: a0 + a1*x + a2*x^2 + a3*x^3 + … + am*x^m

Al igual que anteriores métodos de regresión lineal, el error de regresión polinomial se mide a través del error estándar (Sy/x), coeficiente de determinación (r^2)  y el coeficiente de correlación (r).

Diagrama de flujo

Regresión Polinomial

Código en C++

Link a Github.

rp_codigo1rp_codigo2rp_codigo3rp_codigo4rp_codigo5

 

Ejemplos resueltos con gráficos y cuantificación del error

ejemplo1_grado2

ejemplo1_grado3

ejemplo1_grafica

 


 

ejemplo2_grado2

ejemplo2_grado3.png

ejemplo2_grafica

 


 

ejemplo3_grado2

ejemplo3_grafica

Regresión lineal y linealización

Miguel Angel Montoya Zaragoza – A01226045

Arturo Fornés Arvayo – A01227071

*Se abordara en este post la regresión lineal por mínimos cuadrados, y la linealización de modelos no lineales (Potencial y exponencial)

Antecedentes

La regresión lineal es un método que, a partir de de datos, y la creación de un diagrama de dispersión de los mismos, determina la dependencia funcional de dos variables. Es decir, determina una función que mejor permita modelar el diagrama de dispersión (Con el menor error posible).

Se denomina regresión lineal, si el objetivo de la determinación de la función es encontrar los valores m y b de f(x) = mx + b, es decir, que la función que modelará los datos es lineal.

Si por otro lado, el objetivo es encontrar los valores α y β de f(x) = αx^β, sería una linealización de un modelo no lineal, donde la dispersión se ajusta más a una función potencial, que a una lineal.

O si quisiéramos encontrar los valores α y β de f(x) = αe^(βx), también es una linealización, donde la dispersión se ajusta a una función exponencial.

Más adelante de explicara como es que usando el método de mínimos cuadrados, podemos (Con mínimas modificaciones en código), determinar la recta de regresión lineal o las funciones exponenciales o potenciales .

Justificación y propósito

El propósito de las regresiones lineales o la linealización (De cualquier tipo) es encontrar la influencia de una variable sobre otra, es decir, la dependencia funcional. Crear una función (Recta, potencial o exponencial), para poder predecir valores para x que no se hallan en la dispersión, analizar tendencias y realizar pruebas de hipótesis.

Existen 3 medidas para cuantificar el error de los 3 modelos.

Error estándar (S): Esta medida se podría considerar una desviación estándar, es decir, un promedio de las diferencias (La desviación estándar trabaja alrededor de la media, y el error estándar alrededor de la recta de regresión).

Syx = sqrt(Sr/(n – 2))
Donde Syx = Error estándar, Sr = Suma de diferencias al cuadrado (y real – y medida) alrededor de la recta de regresión, n = número de datos.

Coeficiente de determinación (r^2): Describe, con valores normalizados, que tanto se redujo el error al modelar la dispersión con una recta de regresión, en vez de la media.

r^2 = (St – Sr)/St

Donde r^2 = Coeficiente de determinación, St = Suma total de cuadrados alrededor de la media, Sr = Suma de diferencias al cuadrado (y real – y medida) alrededor de la recta de regresión.

Coeficiente de relación (r): Determina el grado de dependencia entre las dos variables. Entre más se acerca r a 1, indica que la relación es perfecta directa; a -1, perfecta inversa; o 0, inexistente, valores independientes.

r = sqrt((St – Sr)/St)

Donde r = Coeficiente de relación, St = Suma total de cuadrados alrededor de la media, Sr = Suma de diferencias al cuadrado (y real – y medida) alrededor de la recta de regresión.

Tipos de coeficiente de relación:

Correlation_coefficient

Fundamento matemático con gráficos

El método consiste en obtener m(a1) y b(a0) de la ecuación de una recta.
Las formulas son las siguientes:

p

Adaptación formula potencial

f(x) = αx^β

Lo cual se puede representar como log(f(x)) = βlog(x) + log(α) -> F(x) = βX + log(α)
Considerando a0 = log(α) y a1 = β se puede utilizar el método de regresión lineal.

Adaptación formula exponencial

f(x) = αe^βx

Lo cual se puede representar como ln(f(x)) = βx + ln(α) -> F(x) = βx + ln(α)
Considerando a0 = ln(α) y a1 = β se puede utilizar el método de regresión lineal.

Diagrama de flujo

Regresión(1)

Anotación *1: Para linealización potencial, a la lectura de las x y y, se les aplicara log10() antes de usarlas dentro de cualquier proceso. Para linealización exponencial, a la lectura de las y, se les aplicara ln() antes de usarlas dentro de cualquier proceso.

Anotación *2: Para linealización potencial, el calculo de a0 se utiliza como exponente para una base 10, y así obtener α; por otro lado a1 = β. Para linealización exponencial, el calculo de a0 se utiliza como exponente de e, y así obtener α; por otro lado a1 = β.

Anotación *3: Para linealización potencial, el calculo de yMes (y medida) = a0*xReal^a1. Para linealización exponencial, el calculo de yMes (y medida) = a0*e^a1*xReal.

Código en C++

Regresión lineal (Github)

Linealización modelo potencial (Github)

Linealización modelo exponencial (Github)

Ejemplos resueltos con gráficos y cuantificación del error

Ejemplo 1

Gráfica de distribución

lin

Usando regresión linear, se obtiene la recta f(x) = 0.40042x  –  0.203835

p

Ni el modelo potencial o exponencial arroja resultados validos para esta dispersión. Esto significa que la regresión lineal es el único método para modelar la dispersión.

Gráfica y recta

lin

Ejemplo 2

Gráfica de dispersión

Sin título2

Regresión lineal: f(x) = -0.436x + 1.598
Captura
Modelo potencial: f(x) = 0.784511x^(-0.100102)
pot
Modelo exponencial: f(x) = 1.7953e^(-0.496167x)
exp

Gráfica y modelos

Regresión lineal: Rojo
Modelo potencial: Negro
Modelo exponencial: Azul

Sin título

 

Método de Jacobi & Método de Gauss-Seidel

Introducción y antecedentes

Ambos son métodos semejantes, usados para calcular la solución a un sistema de ecuaciones lineales (De forma Ax = b).

Requisitos de los métodos

Para ambos:

  • El sistema de ecuaciones se tiene que poder representar en una matriz cuadrada
  • La diagonal no debe tener elementos nulos
  • La matriz tiene que ser dominante (Los elementos en la diagonal tienen que tener el coeficiente de mayor magnitud de la ecuación)

Diferencias

En Jacobi, para cada iteración, se usan los valores previamente estimados de las variables (O si es la primer iteración, los valores de las variables definidas por el usuario o las default). En cambio en Gauss-Seidel, el método va reemplazando las variables calculadas en la iteración anterior, por las de la iteración actual inmediatamente. Así el método reduce el número de iteraciones totales.

Diagramas de flujo

Se realizaron los diagramas de flujo para resolver un sistema de 3 ecuaciones y 3 incógnitas

Jacobi

Jacobi

Gauss-Seidel

Gauss-Seidel

¿Es un método iterativo?

Sí. Al final de cada iteración, el método ya tiene una aproximación a la solución. Durante la siguiente iteración, se utilizan los resultados previos obtenidos para calcular las nuevas soluciones.

Código

Jacobi: GitHub

Gauss-Seidel: GitHub

Criterio de convergencia del método

Para garantizar la convergencia de ambos métodos, Jacobi y– Gauss-Seidel, se tiene que cumplir con alguna de las siguientes condiciones:

  • Ser diagonal dominante.
  • Ser simétrica y definida positva.

En los siguientes ejemplos nos aseguraremos de cumplir con la propiedad de la diagonal dominante, esto se refiere a que el elemento de la diagonal de cada fila debe ser el elemento de mayor magnitud en esa fila.

Relajación (qué es y para qué sirve?)

Relajación es un método iterativo para resolver sistemas de ecuaciones. Consiste en actualizar los elementos del vector solución, uno a uno, asumiendo que los valores previos están correctos.

Pruebas y resultados

3 incógnitas

3x-2y+z=4

x+5y+2z=6

-2x+3y+6z=1

La solución del sistema es x = 1.6857143, y = 0.7142857, z = 0.3714286.

Usando scilab se encontró la solución:

01-scilab

El sistema cumple con la propiedad de la diagonal dominante, por lo que ambos métodos deberían encontrar la solución.

gs-01-gs.png
Usando Gauss-Seidel
gs-01-j
Usando Jacobi

Ambos métodos encontraron la solución con un error absoluto de 10^-6.


4 incógnitas

8x-3y+z-w=20

-x+20y-3z+8w=3

x+y+3z-w=3

6x-4y+11z+15w=9

La solución del sistema es x = 2.621404, y = 0.3945017, z = – 0.1140893, w = – 0.3263623.

Usando scilab se encontró la solución:

gs-02-matriz

Dado que el sistema cumple con la propiedad de la diagonal dominante, ambos métodos deberían encontrar la solución.

gs-02-gs.png
Usando Gauss-Seidel
gs-02-j
Usando Jacobi

Ambos métodos encontraron la solución con un error absoluto de 10^-6.

Casos de falla

Ambos métodos fallan cuando no se cumple con alguna de las condiciones mencionadas en la sección Criterio de convergencia del método, ambos métodos divergen, por ejemplo, si no se cumple con la propiedad de la diagonal dominante.

Conclusiones

Es clara la superioridad del método de Gauss-Seidel sobre Jacobi. Pues en menos iteraciones, se obtienen las soluciones al sistema. Solamente es necesario aplicar los cambios por cada variable del sistema, más que por cada iteración,

Eliminación de Gauss y Gauss-Jordan

Arturo Fornés – A01227071
Miguel Montoya – A01226045


Introducción y antecedentes

Gauss y Gauss-Jordan son ambos algoritmos del Álgebra Lineal que sirven para encontrar la solución a un sistema lineal (sistema de ecuaciones). Un sistema lineal puede tener una solución única, una infinidad de soluciones o ninguna solución. La solución única de un sistema de ecuaciones se encuentra despejando las incógnitas y encontrando su valor; se tiene una infinidad de soluciones cuando al despejar alguna de las ecuaciones quedan en 0=0; no se tiene solución cuando alguna de las ecuaciones resulta en 0=a, a != 0.

Requisitos de los métodos

  • Una matriz de dimensión n x n.
  • Un sistema de ecuaciones con una solución única.

Diferencias entre uno y otro método

Usando el método de Gauss se obtiene una matriz triangular superior; la última incógnita tiene una solución y se sustituye en la penúltima y se sigue este procedimiento para las siguientes filas en este orden.

Usando Gauss-Jordan se obtiene una matriz en forma escalonada reducida, es decir donde el pivote de cada fila es 1 y todo elemento que no sea el pivote es 0, se forma una diagonal de 1’s. Al llegar a la forma escalonada reducida de la matriz, ya se tiene la solución en el vector con el cual se expandió la matriz de coeficientes para formar el sistema lineal.

Diagramas de flujo

Gauss

Gauss

Gauss-Jordan

Gauss-Jordan

¿Es un método iterativo?

Ninguno de los dos métodos, Gauss ni Gauss-Jordan, es iterativo, hay iteraciones dentro de ellos pero no se aproxima de manera iterativa al resultado. Un método iterativo es definido por realizar una aproximación a la solución en cada una de sus iteraciones; Gauss y Gauss-Jordan calculan una solución una única vez, en ambos casos un vector solución.

Archivo de script en Scilab

Gauss: GitHub.

Screenshot from 2017-03-11 13-27-05

Gauss-Jordan: GitHub.

Screenshot from 2017-03-09 07-40-50

Pivoteo, y escalamiento.

El pivote en una matriz es definido como el primer elemento que sea diferente de 0 en una fila. La forma escalonada se obtiene al seguir un patrón de escalones entre los pivotes, donde el pivote de una fila se encuentre a la derecha del pivote de la fila anterior.

Pruebas y resultados con sistemas de más de 3 incógnitas

Caso 1:

3x + 3y + 12z +21w = 18
2x + 4y + 10z +14w = 16
-2x + -8z -18w = 10
3x + 11z +28w = 11

Gauss

Screenshot from 2017-03-11 13-24-48

Gauss-Jordan

Screenshot from 2017-03-11 13-25-15.png

Ambos métodos obtienen los resultados de manera exitosa

Caso 2:

2x +3y -2z = 7
x -2y +w +4v = -2
10y +4z +5v = 43
7w +7v = 0
3x +5y -5z +11w +v = -2

Gauss

Screenshot from 2017-03-11 14-11-47.png

Gauss-Jordan

Screenshot from 2017-03-11 14-12-18

Ambos métodos resuelven el sistema de 5 variables correctamente

Casos de falla

Ambos métodos, Gauss y Gauss-Jordan, fallarán para sistemas que no tengan una solución única o en los que el pivote de una fila se encuentre a la izquierda del pivote anterior. Lo último se puede evitar reacomodando manualmente las filas (intercambiándolas entre sí).

Conclusiones

El método de Gauss requiere menos pivoteo, pero al mismo tiempo necesita realizar el proceso de sustitución. mientras tanto, Gauss-Jordan elimina por completo la sustitución, obteniendo un valor para cada variable de manera independiente, aunque necesita de mayor pivoteo.

Método de Cramer

Arturo Fornés – A01227071

Miguel Montoya – A01226045


Introducción y antecedentes

La regla de Cramer es un método númerico para obtener las soluciones a un sistema de ecuaciones dado. Es costoso en recursos, sobretodo en matrices grandes.

La regla es: xj = det(Aj)/det(A)

Donde cada Aj es la matriz A donde la columna j es reemplazada por el vector solución

Requisitos del método

  • Para el funcionamiento la matriz tiene que ser cuadrada.
  • El vector solución no deberá ser 0, pues causaría que det(Aj) fuera 0, y por lo tanto la solución. Cuando puede no ser el caso.

Qué es un determinante

Es un termino complejo de definir. El determinante de una matriz es un número, si es 0 es que es un sistema singular, y si es cercano a 0 es que es un sistema mal condicionado. Un sistema singular en un sistema de ecuaciones ocurre cuando varias de las ecuaciones tienen la misma pendiente.

Diagrama de flujo

Cramer

Control de iteraciones del método

Archivo de script en Scilab (GitHub)

Screenshot from 2017-03-08 21-15-05

Pruebas con 2 y 3 incógnitas

2 incógnitas

[  1   1] = [ 2]
[-1   1]    [ 3]

Screenshot from 2017-03-09 20-22-39

El método funciona de manera exitosa

3 incógnitas

[  1   2    3]    [ -5]
[  3   1  -3] = [  4]
[-3   4   7]    [ -7]

Screenshot from 2017-03-09 20-25-06

El método funciona de manera exitosa

Conclusiones

Es un método mucho más fácil de implementar que el método de eliminación  Gaussiana y sustitución. El problema comienza al momento de tener que calcular determinantes de matrices muy grandes. Aunque nunca hay que olvidar que esto permite evitar el pivoteo de Gauss.