El objetivo de los ensayos realizados ha sido caracterizar el comportamiento de un material viscoelástico. Para ello, se ha realizado un ensayo de compresión y un ensayo de cortante, aplicando una carga cíclica en cada uno de los casos. A partir de los valores registrados del ensayo, esto es, carga aplicada y desplazamiento experimentado por la probeta, se obtendrán los parámetros que caracterizan y describen el módulo de Young del material.
Los métodos iterativos o aproximados proveen una alternativa aceptable a los métodos de pivoteo como el método de Gauss y Gauss-Jordan. En el caso que les expongo en esta ocasión, el método de Gauss-Seidel emplea un enfoque similar al utilizado por los métodos para obtener las raíces de ecuaciones simples lineales o no lineales.
Dichos planteamientos se fundamentan en suponer un valor inicial de la raíz y luego utilizar un método sistemático para obtener una aproximación más precisa. En este caso extenderemos este concepto a sistemas de ecuaciones lineales.
El método de Gauss-Seidel, es el método iterativo más comúnmente utilizado para obtener la solución de sistemas de ecuaciones lineales, Supongamos que se nos da en conjunto de $n$ ecuaciones con $n$ incógnitas:
Ahora es posible obtener la solución de los sistemas de ecuaciones presentados en \eqref{ec2}, \eqref{ec1} y \eqref{ec3} a partir de un valor inicial $x$, el cual representa un vector con los valores iniciales para $x_1$, $x_2$ y $x_3$. Lo que se acostumbra a hacer para fijar estos valores iniciales es asumirlos todos iguales a $0$, de esta forma se remplaza $x_2=0$ y $x_3=0$ en la ecuación \eqref{ec2}, obteniéndose de esta manera el primer valor para $x_1$. Luego se reemplaza este nuevo valor de $x_1$ y el valor $x_3=0$ en la ecuación \eqref{ec2} obteniéndose así el nuevo valor de $x_2$. Y por último los nuevos valores de $x_1$ y $x_2$ se reemplazan en la ecuación \eqref{ec3} y se obtiene el nuevo valor estimado de $x_3$. El método iterativo se repite hasta alcanzar una aproximación lo suficientemente precisa.
Como criterio para verificar que el método converge a una solución se utiliza el error relativo, esto es:
\begin{equation}
|\epsilon_{r,i}|=\left | \frac{x_i^j-x_i^{j-1}}{x_i^j}\right |
\end{equation}
Para todas las $i$, en donde $j$ y $j-1$ son las iteraciones actuales y previas respectivamente.
Ahora volvamos al caso general en el que se cuenten con $n$ ecuaciones y $n$ incógnitas. ¿Cómo podemos determinar una expresión general para evaluar los valores de $x$, para de esta forma realizar un algoritmo que luego nos permita generar un programa?. Veamos.
Podemos escribir un conjunto de expresiones como las siguientes:
\begin{eqnarray}
x_i&=&\frac{b_i-\sum_{j=i+1}^n a_{ij}x_j}{a_{ii}}& i=1 \label{ec4}\\
x_i&=&\frac{b_i-\sum_{j=i+1}^n a_{ij}x_j-\sum_{j=1}^{i-1} a_{ij}x_j}{a_{ii}}& i=2,3,..,n-1 \label{ec5}\\
x_i&=&\frac{b_i-\sum_{j=1}^{n-1} a_{ij}x_j}{a_{ii}}& i=n \label{ec6}
\end{eqnarray}
De esta forma se nos es fácil implementar el método en un lenguaje de programación. verifiquemos que sucede con las ecuaciones \eqref{ec4}, \eqref{ec5} y \eqref{ec6} para $n=3$.
De la ecuación \eqref{ec4} tenemos:
\begin{eqnarray}
x_1&=&\frac{b_1-\sum_{j=2}^3 a_{1j}x_j}{a_{11}} \nonumber\\
x_1&=&\frac{b_1-a_{12}x_{2}+a_{13}x_3}{a_{11}} \nonumber
\end{eqnarray}
De la ecuación \eqref{ec5}:
\begin{eqnarray}
x_2&=&\frac{b_2-\sum_{j=3}^3 a_{2j}x_j-\sum_{j=1}^1 a_{2j}x_j}{a_{22}}\nonumber\\
x_2&=&\frac{b_2-a_{23}x_3-a_{21}x_1}{a_{22}}\nonumber
\end{eqnarray}
y de la ecuación \eqref{ec6}:
\begin{eqnarray}
x_3&=&\frac{b_3-\sum_{j=1}^2 a_{3j}x_j}{a_{33}}\nonumber\\
x_3&=&\frac{b_3-a_{31}x_1-a_{32}x_2}{a_{33}}
\end{eqnarray}
Adicionalmente recordemos que los elementos de la diagonal principal deben ser distintos de cero para evitar una división por cero. De modo que se debe tener esto en cuenta a la hora de especificar la matriz de coeficientes $A$.
Cabe mencionar que al ser un método numérico, pudiera darse el caso que el mismo no convergiera. De tal forma que para garantizar una convergencia la matriz $A$ debe cumplir lo siguiente:
Lo cual quiere decir que el elemento perteneciente a la diagonal principal en cada ecuación debe ser mayor que la suma del resto de elementos de la misma ecuación. De esta forma la convergencia se garantiza si se cumple la condición, aunque pueden existir algunos casos en que el método converja sin cumplir con ella. Los sistemas que cumplen con la condición antes mencionada son conocidos como diagonalmente dominantes.
Iteración de Jacobi
Como en el caso del método de Gauss-Seidel cada valor de $x_i$ que se calcula es utilizado inmediatamente para calcular el valor siguiente, de tal forma que si la solución es convergente se utilizará la mejor combinación posible de los valores de $x_i$. Pero también existe otra forma de hacer esto, y es primero reemplazar cada uno de los valores de $x_i$ en las ecuaciones \eqref{ec4}, \eqref{ec5} y \eqref{ec6}, obteniendo otro nuevo conjunto de valores $x_i$, este es el planteamiento conocido como Iteración de Jacobi
Veamos un ejemplo utilizando el método de Gauss-Seidel. Obtener la solución del siguiente sistema de ecuaciones:
\begin{eqnarray}
x_1&+&10x_2&+&4x_3&-&2x_4&=&2 \nonumber\\
8x_1&&&-&10x_3&+&2x_4&=&-7 \nonumber\\
8x_1&+&3x_2&+&x_3&+&17x_4 &=&8\nonumber\\
11x_1&+&7x_2&-&3x_3&+&2x_4&=&12 \nonumber
\end{eqnarray}
Si creamos una función llamada en este caso gseidel.m en matlab, tendríamos lo siguiente:
De esta forma podemos ver que el método converge. A continuación les dejo la función implementada en Matlab Función Gauss-Seidel
Función Gauss-Seidel
%Desarrollada por: Jorge De La Cruz%Para un sistema de la forma 'Ax=b'%imax: Número máxima de iteraciones en caso de%que el método no converja%tol: tolerancia deseada
function x=gseidel(A,b,tol,imax)
n=max(size(A)); % Obtenemos de la matriz A% el número de ecuaiones% y por ende de incognitas.
c=diag(A); % Almacenamos en 'c' los elementos% de la diagonal principal.
d=1;
i=1;
while (d~=0)&&(i<n) % revisamos si en la diagonal principal 'c'
d=c(i,1)*c(i+1,1); % existe algún elemento igual a 0.
i=i+1; % En caso de que lo hayaend% salir del bucle y d=0.while d==0 % Si se ha detectado algún cero en la
B(1,:)=A(n,:); % diagonal principal se entra en este bucle
U(1,1)=b(n,1); % para así cambiar las filas. Se alternaránfor i=2:n % las filas hasta que todos los elementos
B(i,:)=A(i-1,:);% de la diagonal principal
U(i,1)=b(i-1,1);% sean distintos de cero.end% Todo esto para evitar una división
A=B; % por cero. Aunque se recomienda
b=U; % siempre ordenar las ecuaciones apriori.
c=diag(A);
d=1;i=1;
while (d~=0)&&(i<n)
d=c(i,1)*c(i+1,1);
i=i+1;
endend% Se selecciona un vector X(n,1) sólo para evaluar% el error en la primera iteración.for i=1:n
X(i,1)=b(i,1)/A(i,i);
end% se sespecifica el vector x(n,1) con valores iguales a cero% para iniciar las iteraciones.
x=zeros(n,1);
err=tol+1;
iter=1;
% Se inician las iteracioneswhile (tol<err)&&(iter<imax)
for i=1:n
if i==1
sum=0;
for j=i+1:n
sum=A(i,j)*x(j,1)+sum;
end
x(i,1)=(b(i,1)-sum)/A(i,i);
endif i==n
sum=0;
for j=i-1:-1:1
sum=A(i,j)*x(j,1)+sum;
end
x(i,1)=(b(i,1)-sum)/A(i,i);
endif 2<=i<=n-1
sum1=0;
for j=i+1:n
sum1=A(i,j)*x(j,1)+sum1;
end
sum2=0;
for j=i-1:-1:1
sum2=A(i,j)*x(j,1)+sum2;
end
x(i,1)=(b(i,1)-sum1-sum2)/A(i,i);
endend
err=abs((x-X)./x);
err=max(err);
X=x;
iter=iter+1;
end
Uno de los problemas de la implementación del método de Newton-Raphson es el de la evaluación de la derivada. Puede que esto no sea un inconveniente para los polinomios y algunas otras funciones, pero existen algunas otras funciones para las cuales la derivada puede ser difícil de evaluar. Además también pudiera darse el caso en que en realidad no contemos con una función como tal sino que solo dispongamos de un conjunto de puntos $(x_i,f(x_i))$ de tal manera que tampoco sería posible obtener de forma analítica la derivada de nuestra función objetivo.Para casos como estos, la derivada su puede aproximar como una diferencia dividida finita regresiva. La derivada se puede aproximar entonces de la siguiente forma:
De esta forma la ecuación \eqref{ec3} es la formula para el Método de la Secante. Nótese que en este caso a diferencia del método de Newton-Raphson es necesario dos puntos para poder iniciar el método.
Figura: Ilustración de cómo el Método de la Secante se aproxima a la raiz.
Método de la Secante
Método de la Secante
%///////////////////////////////////////////////////%// Función Secante ///%// Desarrollada por Jorge De La Cruz ///%/ ///%//////////////////////////////////////////////////
function xr=secante(fun,x0,x1,tol,imax)
x=[x0 x1];
f=eval(fun);
err=tol+1;
n=1;
while (err>tol)&&(n<imax)
xr=x(1,2)-f(1,2)*(x(1,1)-x(1,2))/(f(1,1)-f(1,2));
err=abs((xr-x(1,2))/xr);
error=err*100;
x(1,1)=x(1,2);
x(1,2)=xr;
f=eval(fun);
fprintf('Iteración: %d, Error: %f, Raiz: %f\n',n,error,xr)
n=n+1;
end
Entradas
