0

He intentado resolver un sistema de ecuaciones por el método de Jacobi (Previamente se hizo el cálculo por el método LU) Tengo que obtener el número de iteraciones y guardarlo en un archivo externo. Para ello creo una variable llamada k que me guardaría el número de iteraciones. Sin embargo, he intentado emplear la k como un valor unidimensional de tipo int pero al llamarlo en mi funcion main esta queda como valor 0. Sin embargo, si declaro ahora k como un array unidimensional k[1] si que he conseguido obtener el número de iteraciones. Es como si el void solo me guarda los valores tipo array pero si no lo declaro así no consigo resolverlo. Alguna idea?

#include <iostream>
#include <cmath>
#include <fstream>
#include<iomanip>

using namespace std;
const int n=5;
void lu(double A[n][n], double b[n], double x[n], const int n);
void Jacobi(double A[n][n], double B[n], double x[n], double xl[n], int n, double tol,int k[]);
float norma(double x[n],int n);
int main()
{
    //Vamos a cambiar las posiciones de la matriz que nos dan y lo vamos a poner de tal forma que en la diagonal esté los números con mayor valor
    //Nuestra matriz no es estritamente dominante por lo que es posible que no converja. El cambio que vamos a hacer es cambiar la tercera fila con la cuarta
    double A[n][n]={{12, -2, 1, 4, -0.5},{-1, 15, 4, 6, 7},{9, 0.5, 12, 3, 1},{-3, -6, 9, 24 ,1},{1, -1, 3, 2, 8}};
    double b[n]={-1, 21, 43, 32, 55};
    double x[n]; //array que guarda las soluciones para el método LU
    lu(A,b,x,n);
    cout<< "La solucion obtenida por el metodo LU es: "<< endl;
    //hemos reordenado la matriz luego la x3 obtenida corresponde a la x4 original y viceversa
    for(int i=0;i<n;i++){
            cout << "x("<< i << ") es: "<< x[i]<< endl;
    }
    double xj[n]; //array que guarda las soluciones del método de Jacobi
    double epsilon=1e-3;
    int k[1];
    for(int i=0; i<10; i++){
        epsilon= epsilon*0.1;
        Jacobi(A,b,xj,x,n,epsilon,k);
    cout<< "El numero de iteraciones realizado por Jacobi es: "<< k[0]<< "para epsilon= "<< epsilon << endl;
    cout<< "La solucion hallada por Jacobi es: "
    for(int i=0;i<n;i++){
            cout << "x("<< i << ") es: "<< xj[i]<< endl;
    }
}
    return 0;
}

void lu(double A[n][n], double b[n], double x[n], const int n){    //Declaramos la funcion que nos calcula la matriz Lower Upper
    // A[n][n]: Matriz principal b[n] términos independientes x[n] solución n [número de incógnitas]
    double L[n][n], U[n][n], z[n]; //Matriz Lower-Upper y vector z
    int i, j, k;
    for(i=0;i<n;i++){
        for(j=0; j<n; j++){
            if(i<j){
                double sumu= 0; //Sumatoria para la matriz upper
                L[i][j]=0;
                for(k=0; k<i; k++){
                    sumu=sumu+L[i][k]*U[k][j];
                }
            U[i][j]= A[i][j]-sumu;
            }
        else if(i==j){
        double sumu= 0;
        L[i][i]=1.0;
            for(k=0; k<i; k++){
                    sumu=sumu+L[i][k]*U[k][j];
                }
                U[i][i]= A[i][i]-sumu;
        }
        else{
            U[i][j]=0.0;
            double suml=0.0;
            for(k=0; k<j; k++){
                    suml=suml+L[i][k]*U[k][j];
                }
                L[i][j]=(A[i][j]-suml)/(U[j][j]);

        }
        }
    }
    cout<< endl;
    cout << "La matriz Lower es: "<< endl;
    for (int i=0;i<=n-1;i++){
            for (int k=0;k<=n-1;k++){
                cout << L[i][k] << " " ;}
    cout << endl;}
    cout<< endl;
    cout << "La matriz Upper es: "<< endl;
        for (int i=0;i<=n-1;i++){
            for (int k=0;k<=n-1;k++){
                cout << U[i][k] << " " ;}
    cout << endl;}
    cout<< endl;
    //Vamos a obtener z
    cout<< "El valor de la z es: "<<endl;
    for(int k=0; k<n; k++){
            double sumz=0;
        for (int l=0; l<k; l++){
            sumz=sumz+L[k][l]*z[l];
        }
        z[k]=b[k]-sumz;
        cout << z[k] << endl;
    }
    for(int i=n-1;i>=0;i--){
        double sum=0.0;
        for(j=i+1;j<n;j++){
            sum=sum+U[i][j]*x[j];
        }
        x[i]=(z[i]-sum)/U[i][i];
    }

}
float norma(double x[n],int n){
    float s=0;
        for(int i=0;i<n;i++){
            s += x[i]*x[i];
        }
        double sq= sqrt(s);
        return sq;
}

void Jacobi(double A[n][n], double B[n], double x[n], double xl[n], int n, double tol, int k[]){
    //Escogemos unos valores iniciales para x[i], para ello vamos a buscar solución teniendo en cuenta solo términos de la diagonal pricipal
    for(int l=0;l<n;l++){
 x[l]=B[l]/A[l][l];
 cout << endl << "La primera aproximacion es x("<<l<<") :";
 cout<< x[l]<< endl;
    }
    k[0]=0;
    double xold[n];
    do{
        for (int i=0; i<n; i++){
xold[i]=x[i];
 }
for(int i=0;i<n;i++){
    double sigma=0;
    for(int j=0;j<n;j++){
        if(j != i){
        sigma=sigma+A[i][j]*xold[j];
        }
    }
    x[i]=(B[i]-sigma)/A[i][i];
}
k[0]=k[0]+1;
}while((fabs(norma(x,n)-norma(xl,n)))>tol);

}
2

0

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

Examina otras preguntas con la etiqueta o formula tu propia pregunta.