Método de Jacobi programado en C/C++ para la resolución de sistemas lineares de ecuaciones:
El método Jacobi es el método iterativo para resolver sistemas de ecuaciones lineales más simple y se aplica sólo a sistemas cuadrados, es decir a sistemas con tantas incógnitas como ecuaciones, del tipo:
Ax=b
Idea:
1-Separar la matriz A en su diagonal (D) y en su matrices estrictamente superior (U) e inferior (L):
A=D+L+U=D+R
2-Reescribir el sistema con este desarrollo de la matriz A:
b-(D+R)x=0
b-Rx=Dx
D^(-1)(b-Rx)=x
3- A partir de esta ultima expresion se itera tomando un valor inicial de x denominado x0:
D^(-1)(b-Rx(i))=x(i+1)
donde x(i) hace referencia al valor de x en la iteración i.
Hay que notar que la complejidad de calcular la inversa de la matriz diagonal es nula:
D^(-1)(i,i)=1/D^(-1)(i,i)
Nota: El método de Jacobi siempre converge si la matriz A es estrictamente diagonal dominante y puede converger incluso si esta condición no se satisface.
Idea:
1-Separar la matriz A en su diagonal (D) y en su matrices estrictamente superior (U) e inferior (L):
A=D+L+U=D+R
2-Reescribir el sistema con este desarrollo de la matriz A:
b-(D+R)x=0
b-Rx=Dx
D^(-1)(b-Rx)=x
3- A partir de esta ultima expresion se itera tomando un valor inicial de x denominado x0:
D^(-1)(b-Rx(i))=x(i+1)
donde x(i) hace referencia al valor de x en la iteración i.
Hay que notar que la complejidad de calcular la inversa de la matriz diagonal es nula:
D^(-1)(i,i)=1/D^(-1)(i,i)
Nota: El método de Jacobi siempre converge si la matriz A es estrictamente diagonal dominante y puede converger incluso si esta condición no se satisface.
Dejo aquí los enlaces a otros metodos de resolucion de sistemas programados en C y en PASCAL:
JACOBI en PASCAL: http://cypascal.blogspot.fr/2012/03/resolucion-de-sistemas-mediante-jacobi.html
GAUSS en C: http://cypascal.blogspot.fr/2012/10/metodo-de-gauss-en-cc.html
GAUSS en PASCAL: http://cypascal.blogspot.fr/2012/10/metodo-de-gauss-en-cc.html
LU en PASCAL:http://cypascal.blogspot.fr/2012/02/factorizacion-lu.html
Código:
#include <stdio.h>
#include <stdlib.h>
int dim;
float norma(float vector1[],float vector2[]);
float suma_jacobi(float Matriz[], float vector[], int componente);
int main(){
int i,j,iteraciones=0;
float error,epsilon;
printf("\n METODO DE JACOBI DE RESOLUCION DE SISTEMAS Ax=b \n");
printf("Dimension de la matriz A: ");
scanf("%d",&dim);
float A[dim][dim],b[dim],x[dim],x_prev[dim],aux[dim];
printf("\n Elementos de la matriz A: \n");
for(i=0;i<dim;i++) for(j=0;j<dim;j++){
printf("A(%d,%d)=",i,j); scanf("%f",&A[i][j]);
}
printf("\n Elementos del vector b: \n");
for(i=0;i<dim;i++){
printf("b(%d)=",i); scanf("%f",&b[i]);
}
printf("\n Error de parada: \n");
printf("E=",i); scanf("%f",&epsilon);
error=epsilon+1;
//cominezo algoritmo de Jacobi
//Error se mide como la norma del vector diferenceia entre la iteracion i e i+1
printf("\n Valor inicail de la iteracion: \n");
for(i=0;i<dim;i++){
printf("x0(%d)=",i); scanf("%f",&x_prev[i]);
}
while (error>epsilon){
for(i=0;i<dim;i++){
for(j=0;j<dim;j++) aux[j]=A[i][j];
x[i]=(1/A[i][i])*(b[i]-suma_jacobi(aux,x_prev,i));
}
error=norma(x,x_prev);
printf("\n\n Iteracion %d: \n",iteraciones);
for(i=0;i<dim;i++){
x_prev[i]=x[i];
printf("X(%d)=%f \n",i,x[i]);
}
iteraciones++;
if (iteraciones==10) error=epsilon-1;
}
printf("Solucion del sistema\n");
printf("Numero de iteraciones: %d \n", iteraciones);
for(i=0;i<dim;i++){
printf("x(%d)=%f\n",i,x[i]);
}
return 1;
}
float norma(float vector1[],float vector2[]){
float aux=0;
int i;
for(i=0;i<dim;i++){
aux=aux+(vector1[i]-vector2[i])*(vector1[i]-vector2[i]);
}
return aux;
}
float suma_jacobi(float Matriz[], float vector[], int componente)
{
float aux=0;
int i;
for(i=0;i<dim;i++){
if (componente!=i){
aux=aux+Matriz[i]*vector[i];
}
}
return aux;
}
Salu10
#include <stdlib.h>
int dim;
float norma(float vector1[],float vector2[]);
float suma_jacobi(float Matriz[], float vector[], int componente);
int main(){
int i,j,iteraciones=0;
float error,epsilon;
printf("\n METODO DE JACOBI DE RESOLUCION DE SISTEMAS Ax=b \n");
printf("Dimension de la matriz A: ");
scanf("%d",&dim);
float A[dim][dim],b[dim],x[dim],x_prev[dim],aux[dim];
printf("\n Elementos de la matriz A: \n");
for(i=0;i<dim;i++) for(j=0;j<dim;j++){
printf("A(%d,%d)=",i,j); scanf("%f",&A[i][j]);
}
printf("\n Elementos del vector b: \n");
for(i=0;i<dim;i++){
printf("b(%d)=",i); scanf("%f",&b[i]);
}
printf("\n Error de parada: \n");
printf("E=",i); scanf("%f",&epsilon);
error=epsilon+1;
//cominezo algoritmo de Jacobi
//Error se mide como la norma del vector diferenceia entre la iteracion i e i+1
printf("\n Valor inicail de la iteracion: \n");
for(i=0;i<dim;i++){
printf("x0(%d)=",i); scanf("%f",&x_prev[i]);
}
while (error>epsilon){
for(i=0;i<dim;i++){
for(j=0;j<dim;j++) aux[j]=A[i][j];
x[i]=(1/A[i][i])*(b[i]-suma_jacobi(aux,x_prev,i));
}
error=norma(x,x_prev);
printf("\n\n Iteracion %d: \n",iteraciones);
for(i=0;i<dim;i++){
x_prev[i]=x[i];
printf("X(%d)=%f \n",i,x[i]);
}
iteraciones++;
if (iteraciones==10) error=epsilon-1;
}
printf("Solucion del sistema\n");
printf("Numero de iteraciones: %d \n", iteraciones);
for(i=0;i<dim;i++){
printf("x(%d)=%f\n",i,x[i]);
}
return 1;
}
float norma(float vector1[],float vector2[]){
float aux=0;
int i;
for(i=0;i<dim;i++){
aux=aux+(vector1[i]-vector2[i])*(vector1[i]-vector2[i]);
}
return aux;
}
float suma_jacobi(float Matriz[], float vector[], int componente)
{
float aux=0;
int i;
for(i=0;i<dim;i++){
if (componente!=i){
aux=aux+Matriz[i]*vector[i];
}
}
return aux;
}
Salu10
