Matemáticas

Mínimos Cuadrados No Lineales

Posted by

7

Para ajustar un modelo m\left( x,t\right) de n parámetros a un conjunto de m datos \left( t_{i},y_{i}\right) , definimos f\left( x\right) =\frac{1}{2}R\left( x\right) ^{T}R\left( x\right) con R:\mathbb{R}^{n}\rightarrow \mathbb{R}^{m} y R\left( x\right) no lineal con r_{i}\left( x\right) =m\left( x,t_{i}\right) -y_{i} cómo i-ésima componente.

El problema finalmente es:

Minimizar f\left( x\right) =\frac{1}{2}R\left( x\right) ^{T}R\left( x\right) =\frac{1}{2}\sum\limits_{i=1}^{m}r_{i}^{2}\left( x\right)

Esto es equivalente a minimizar la suma de los cuadrados de las diferencias entre el modelo y los datos.

f\left( x\right) es una función semi-paraboidal en el sentido de que es la suma de parábolas y por tanto es minimizable.

Para hallar x tal que f\left( x\right) se minimiza, hallamos las soluciones del gradiente de f\left( x\right) , de este modo, definiendo J\left( x\right) \in \mathbb{R} ^{m\times n} como el jacobiano de R\left( x\right) ,J\left( x\right) _{ij}=\frac{\partial r_{i}\left( x\right) }{\partial x_{j}} tenemos las 2 primeras derivadas de f\left( x\right) :

\nabla f\left( x\right) =\sum\limits_{i=1}^{m}r_{i}\left( x\right) \cdot \nabla r_{i}\left( x\right) =J\left( x\right) ^{T}R\left( x\right)

\nabla ^{2}f\left( x\right) =\sum\limits_{i=1}^{m}\left( \nabla r_{i}\left( x\right) \cdot \nabla r_{i}\left( x\right) ^{T}+r_{i}\left( x\right) \cdot \nabla ^{2}r_{i}\left( x\right) \right) =J\left( x\right) ^{T}J\left( x\right) +S\left( x\right)

Siendo S\left( x\right) =\sum\limits_{i=1}^{m}r_{i}\left( x\right) \cdot \nabla ^{2}r_{i}\left( x\right)

Aplicando el método de Newton para resolver \nabla f\left( x\right) tenemos:

x_{n+1}=x_{n}-\left( \nabla ^{2}f\left( x_{n}\right) \right) ^{-1}\nabla f\left( x_{n}\right)

x_{n+1}=x_{n}-\left( J\left( x_{n}\right) ^{T}J\left( x_{n}\right) +S\left( x_{n}\right) \right) ^{-1}J\left( x_{n}\right) ^{T}R\left( x_{n}\right)

EJEMPLO

Ajustar el modelo y\left( t\right) =e^{tx} a los siguientes datos:

\begin{tabular}{|l|l|l|l|} \hline t & 1 & 2 & 3 \\ \hline y & 2 & 4 & 8 \\ \hline \end{tabular}

De aquí, m\left( x,t_{i}\right) =e^{t_{i}x} entonces:

R=r_{i}\left( x\right) =m\left( x,t_{i}\right) -y_{0}\left( t_{i}\right) =\left( \begin{array}{c} e^{x}-2 \\ e^{2x}-4 \\ e^{3x}-8\end{array}\right)

Ahora f\left( x\right) =\frac{1}{2}R\left( x\right) ^{T}R\left( x\right) =\frac{1}{2}\left[ \left( e^{x}-2\right) ^{2}+\left( e^{2x}-4\right) ^{2}+\left( e^{3x}-8\right) ^{2}\right]

De este modo:

\nabla f\left( x\right) =e^{x}\left( e^{x}-2\right) +2e^{2x}\left( e^{2x}-4\right) +3e^{3x}\left( e^{3x}-8\right)

\nabla ^{2}f\left( x\right) =\left( e^{x}\right) ^{2}+\left( 2e^{2x}\right) ^{2}+\left( 3e^{3x}\right) ^{2}+e^{x}\left( e^{x}-2\right) +4e^{2x}\left( e^{2x}-4\right) +9e^{3x}\left( e^{3x}-8\right)

Finalmente aplicando el método de Newton con:

x_{n+1}=x_{n}-\left( \nabla ^{2}f\left( x_{n}\right) \right) ^{-1}\nabla f\left( x_{n}\right)

Con x_{0}=1 y 8 decimales tenemos:

\begin{tabular}{|l|l|} \hline n & x \\ \hline 0 & 1.00000000 \\ \hline 1 & 0.87299182 \\ \hline 2 & 0.77360001 \\ \hline 3 & 0.71421136 \\ \hline 4 & 0.69491616 \\ \hline 5 & 0.69316064 \\ \hline 6 & 0.69314718 \\ \hline 7 & 0.69314718 \\ \hline \end{tabular}

Obtenemos una convergencia en 7 iteraciones con x_{7}=\ln 2\approx 0.69314718

De modo que y\left( t\right) =e^{t\ln 2}=2^{t} la cual se ajusta a los datos.

Post #10 – mycomplexsoul

Integrales Exponenciales

Posted by

2

La integral exponencial simple \int x^{n}e^{ax}dx con n\geq 0 y a
constante tiene por resultado:

\int x^{n}e^{ax}dx=e^{ax}\sum\limits_{i=0}^{n}\left( -1\right) ^{i}\frac{1}{a^{i+1}}\frac{d^{i}}{dx^{i}}\left( x^{n}\right)

Donde \frac{d^{0}}{dx^{0}}\left( x^{n}\right) =x^{n}.

Para demostrar este resultado usamos inducción sobre n.

DEMOSTRACIÓN

Para n=0.

\int x^{0}e^{ax}dx=e^{ax}\sum\limits_{i=0}^{0}\left( -1\right) ^{i}\frac{1}{a^{i+1}}\frac{d^{i}}{dx^{i}}\left( x^{0}\right) =e^{ax}\left( \frac{1}{a}\right)

Suponemos para n=k.

\int x^{k}e^{ax}dx=e^{ax}\sum\limits_{i=0}^{k}\left( -1\right) ^{i}\frac{1}{a^{i+1}}\frac{d^{i}}{dx^{i}}\left( x^{k}\right)

Probamos para n=k+1. Usando integración por partes.

\int x^{k+1}e^{ax}dx=\frac{1}{a}x^{k+1}e^{ax}-\frac{\left( k+1\right) }{a}\int x^{k}e^{ax}dx

Sustituyendo

\int x^{k+1}e^{ax}dx=\frac{1}{a}x^{k+1}e^{ax}-\frac{\left( k+1\right) }{a}e^{ax}\sum\limits_{i=0}^{k}\left( -1\right) ^{i}\frac{1}{a^{i+1}}\frac{d^{i}}{dx^{i}}\left( x^{k}\right)

Reescribiendo el índice de la sumatoria

\int x^{k+1}e^{ax}dx=e^{ax}\left( \frac{1}{a}x^{k+1}-\frac{\left( k+1\right) }{a}\sum\limits_{j=1}^{k+1}\left( -1\right) ^{j-1}\frac{1}{a^{j}}\frac{d^{j-1}}{dx^{j-1}}\left( x^{k}\right) \right)

Reescribiendo

\int x^{k+1}e^{ax}dx=e^{ax}\left( \frac{1}{a}x^{k+1}+\sum\limits_{j=1}^{k+1}\left( -1\right) ^{j}\frac{1}{a^{j+1}}\frac{d^{j-1}}{dx^{j-1}}\left( \left( k+1\right) x^{k}\right) \right)

\int x^{k+1}e^{ax}dx=e^{ax}\left( \frac{1}{a}x^{k+1}+\sum\limits_{j=1}^{k+1}\left( -1\right) ^{j}\frac{1}{a^{j+1}}\frac{d^{j-1}}{dx^{j-1}}\left( \frac{d}{dx}\left( x^{k+1}\right) \right) \right)

\int x^{k+1}e^{ax}dx=e^{ax}\left( \frac{1}{a}x^{k+1}+\sum\limits_{j=1}^{k+1}\left( -1\right) ^{j}\frac{1}{a^{j+1}}\frac{d^{j}}{dx^{j}}\left( x^{k+1}\right) \right)

Introduciendo el termino restante con j=0.

\int x^{k+1}e^{ax}dx=e^{ax}\sum\limits_{j=0}^{k+1}\left( -1\right) ^{j}\frac{1}{a^{j+1}}\frac{d^{j}}{dx^{j}}\left( x^{k+1}\right)

Demostrado \forall n\in \mathbb{N}.

Post #02 – mycomplexsoul

Raíz Cúbica

Posted by

9

El algoritmo para calcular la raíz cúbica a mano es idéntico al de la raíz cuadrada salvo las diferencias en la expansión binomial \left( a+b\right) ^{3} respecto a \left( a+b\right) ^{2}.

En el algoritmo de raíz cuadrada, se usa la expansión \left(a+b\right) ^{2}=\allowbreak a^{2}+2ab+b^{2} de tal forma que al separar los digitos de a+b en pares podamos aproximar el primer numero a dejando el resto de la expresión para repeticiones sucesivas del algoritmo. De este modo cuando en el siguiente paso hacemos «el doble del numero que elegimos agregando otro tal que multiplicado por el numero que agregamos sea menor que el resto dentro de la casita» es exactamente lo mismo que hacer \left(2a\ast 10+b\right) \ast b=2ab+b^{2} que es el complemento de la expansión binomial (con la ligera diferencia de que $a$ es un multiplo de 10). Para la raíz cúbica se hará lo mismo en cada paso, entonces:

ALGORITMO

1.- Separamos los digitos de a+b en ternas de derecha a izquierda a partir del punto decimal y después del punto, de izquierda a derecha.

2.- Se aproxima por abajo la raíz cúbica del numero (en bloques de comas) más a la izquierda y se anota en el resultado calculando el resto y agregando la siguiente terna.

3.- Usando el resultado actual multiplicado por 10 como a y agregando un b de tal forma que la expresión 3a^{2}b+3ab^{2}+b^{3} aproxime por abajo al resto actual. Calculamos el resto, agregamos la siguiente terna (si la hay).

4.- Se repite el paso 3 hasta terminar con las ternas o hasta la precisión deseada.

EJEMPLO \sqrt[3]{15363967256}=2486

Paso 1: Separamos los digitos en ternas.

\sqrt[3]{15,363,967,256}

Paso 2: Aproximamos \sqrt[3]{15} por debajo, esto es 2. De modo que 15-2^{3}=7, lo anotamos y agregamos la siguiente terna de arriba 363.

\begin{array}{ll}\sqrt[3]{15,363,967,256} & 2 \\\text{ \ \ \ }7,363 &\end{array}

Paso 3: Hacemos a=2\ast 10=20 y pensamos un b tal que 3a^{2}b+3ab^{2}+b^{3}\leq 7,363. Reexpresando tenemos: 3\left( 20\right)^{2}b+3\left( 20\right) b^{2}+b^{3}=\allowbreak 1200b+60b^{2}+b^{3}, donde el termino «lider» es 1200b seleccionando b=4 obtenemos 1200\left( 4\right) +60\left( 4\right) ^{2}+4^{3}=\allowbreak 5,824. Anotamos el 4 (el valor de b que elegimos), calculamos el resto 7363-5824=1539, lo anotamos debajo y bajamos la siguiente terna. (Note que con b=5 nos hubieramos pasado del resto obteniendo 1200\left( 5\right) +60\left(5\right) ^{2}+5^{3}=\allowbreak 7625).

\begin{array}{ll}\sqrt[3]{15,363,967,256} & 24 \\\text{ \ \ \ }7,363 & 3\left( 20\right) ^{2}\left( 4\right) +3\left(20\right) \left( 4\right) ^{2}+\left( 4\right) ^{3}=\allowbreak 5,824 \\\text{ \ \ \ }1,539,967 &\end{array}

Paso 3 (Segunda iteración): Ahora hacemos a=24\ast 10=240 y pensamos un b tal que 3a^{2}b+3ab^{2}+b^{3}\leq 1,539,967. Reexpresando: 3\left( 240\right) ^{2}b+3\left( 240\right) b^{2}+b^{3}=\allowbreak 172800b+720b^{2}+b^{3}, donde el termino «lider» es 172800b seleccionando b=8 obtenemos 172800\left( 8\right) +720\left( 8\right)^{2}+8^{3}=\allowbreak 1428\,992. Anotamos el 8 (el valor de b que elegimos), calculamos el resto 1,539,967-1,428,992=110,975, lo anotamos debajo y bajamos la siguiente terna.

\begin{array}{ll}\sqrt[3]{15,363,967,256} & 248 \\\text{ \ \ \ }7,363 & 3\left( 20\right) ^{2}\left( 4\right) +3\left(20\right) \left( 4\right) ^{2}+\left( 4\right) ^{3}=\allowbreak 5,824 \\\text{ \ \ \ }1,539,967 & 3\left( 240\right) ^{2}\left( 8\right) +3\left(240\right) \left( 8\right) ^{2}+\left( 8\right) ^{3}=1,428\,,992 \\\text{ \ \ \ \ \ \ }110,975,256 &\end{array}

Paso 3 (Tercera iteración): Ahora hacemos a=248\ast 10=2480 y pensamos un b tal que 3a^{2}b+3ab^{2}+b^{3}\leq 110,975,256. Reexpresando tenemos: 3\left( 2480\right) ^{2}b+3\left( 2480\right) b^{2}+b^{3}, donde el termino «lider» es 3\left( 2480\right) ^{2}b seleccionando b=6 obtenemos 3\left( 2480\right) ^{2}\left( 6\right) +3\left( 2480\right)\left( 6\right) ^{2}+\left( 6\right) ^{3}=\allowbreak 110\,975\,256. Anotamos el 6 (el valor de $b$ que elegimos), calculamos el resto 110,975,256-11,097,256=0, lo anotamos debajo y como el resto es cero, la raíz es exacta.

\begin{array}{ll}\sqrt[3]{15,363,967,256} & 2486 \\\text{ \ \ \ }7,363 & 3\left( 20\right) ^{2}\left( 4\right) +3\left(20\right) \left( 4\right) ^{2}+\left( 4\right) ^{3}=\allowbreak 5,824 \\\text{ \ \ \ }1,539,967 & 3\left( 240\right) ^{2}\left( 8\right) +3\left(240\right) \left( 8\right) ^{2}+\left( 8\right) ^{3}=1,428\,,992 \\\text{ \ \ \ \ \ \ }110,975,256 & 3\left( 2480\right) ^{2}\left( 6\right)+3\left( 2480\right) \left( 6\right) ^{2}+\left( 6\right) ^{3}=\allowbreak 110\,975\,256 \\\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }0 &\end{array}

De este modo:

\sqrt[3]{15,363,967,256}=2486

Como este procedimiento es general, para cualquier entero n>3 basta con sustituir las referencias cúbicas por n-ésimas y la expansión binomial cúbica por la n-ésima \left( a+b\right) ^{n}.

Post #01 – mycomplexsoul