miércoles, 17 de octubre de 2012

MINIMOS CUADRADOS EN C/C++

El siguiente programa trata del problema de minimos cuadrados. Dados unos datos experimentales, se busca una correlacion que relacione ambas variables medainte un polinomio de grado n. Este metodo nos permite ajustar datos experimentales a rectas, parabolas, polinomios de grado superior... obteniendo asi el menor error cudratico. El programa usa funciones como ResolverGauss del programa (Ver alguno de los dos enlaces):    http://cypascal.blogspot.com.es/2012/10/metodo-de-gauss-en-cc.html   o la función potencia, que calcula la potencia de un numero en coma flotante elevado a un entero.
La ejecución del programa es la siguiente:






El codigo del programa de ajuste por minimos cuadrados es este:

#include <stdio.h>

void PideDatos(int *Dat, int *Ord, float Val[][102]);
float Potencia(int n, float Num);
void PreparaSistema(int Ord, int Dat, float Sist[][102], float Val[][102]);
void ResuelveGauss(int Dim, float Sist[][102]);
void EscribeDatos(int Dim, float Sist[][102]);


int main(void)
{
    int Datos,Orden,C;
   
    float Valores[2][102],Sistema[102][102];
    PideDatos(&Datos,&Orden,Valores);
    PreparaSistema(Orden,Datos,Sistema,Valores);
    printf("\n\n El sistema a resolver es el siguiente:");
    EscribeDatos(Orden,Sistema);
    ResuelveGauss(Orden,Sistema);
    printf("\n\n El sistema resuelto:");
    EscribeDatos(Orden,Sistema);
    printf("\n\n La Ecuacion del Polinomio ajustado por minimos Cuadrados\n\n:");
    for(C=1;C<=Orden;C++)
        printf(" + (%f)X^%d",Sistema[C][Orden+1],C-1);
    return(0);
}



void PideDatos(int *Dat, int *Ord,float Val[][102])
{
    int A,B;
   
    printf("\n\n\n METODO DE MINIMOS CUADRADOS.\n\n");
    printf("\n Introduce el numero de datos (Puntos): ");scanf("%d",&*Dat);
    printf("\n\n\n Introce los valores de cada punto\n");
   
    for(A=1;A<=*Dat;A++)
    {
        printf(" -Valores del Punto %d:\n",A);
        printf("   X%d: ",A); scanf("%f",&Val[0][A]);
        printf("   Y%d: ",A); scanf("%f",&Val[1][A]);
    }
    printf("\n\n\n Introduce el orden del polinomio: "); scanf("%d",&B);
    *Ord=B+1;
}


float Potencia(int n, float Num)
{
    int A;
    float res;
   
    res=1;
    for(A=1;A<=n;A++) res=res*Num;
    return res;
}



void PreparaSistema(int Ord, int Dat, float Sist[][102], float Val[][102])
{
    int A,B,C,Exp;
    float suma,termino;
   
    for(A=1;A<=Ord;A++)    for(B=1;B<=Ord;B++)
    {
        suma=0;
        Exp=A+B-2;
      
        for(C=1;C<=Dat;C++)
        {
            termino=Val[0][C];
            suma=suma+Potencia(Exp,termino);
        }
        Sist[A][B]=suma;
    }
    for(A=1;A<=Ord;A++)
    {
        suma=0;
        Exp=A-1;
      
        for(C=1;C<=Dat;C++)
        {
            termino=Val[0][C];
            suma=suma+Val[1][C]*Potencia(Exp,termino);
        }
        Sist[A][Ord+1]=suma;
    }
}


void ResuelveGauss(int Dim, float Sist[][102])
{
    int NoCero,Col,C1,C2,A;
    float Pivote,V1;
   
    for(Col=1;Col<=Dim;Col++){
        NoCero=0;A=Col;
        while(NoCero==0){
            if(Sist[A][Col]!=0){
                NoCero=1;}
            else A++;}
        Pivote=Sist[A][Col];
        for(C1=1;C1<=(Dim+1);C1++){
            V1=Sist[A][C1];
            Sist[A][C1]=Sist[Col][C1];
            Sist[Col][C1]=V1/Pivote;}
        for(C2=Col+1;C2<=Dim;C2++){
            V1=Sist[C2][Col];
            for(C1=Col;C1<=(Dim+1);C1++){
                Sist[C2][C1]=Sist[C2][C1]-V1*Sist[Col][C1];}
        }}
   
    for(Col=Dim;Col>=1;Col--) for(C1=(Col-1);C1>=1;C1--){
        Sist[C1][Dim+1]=Sist[C1][Dim+1]-Sist[C1][Col]*Sist[Col][Dim+1];
        Sist[C1][Col]=0;
    }
}

void EscribeDatos(int Dim, float Sist[][102])
{
    int A,B;
    printf("\n\n");
    for(A=1;A<=Dim;A++){
        for(B=1;B<=(Dim+1);B++){
            printf("%7.2f",Sist[A][B]);
            if(B==Dim) printf("   |");}
        printf("\n");
    }
    printf("\n\n");
}










12 comentarios:

  1. Fantastico amigo, no tendras el algoritmo para diferencias divididas con 6 datos en matlab

    ResponderEliminar
    Respuestas
    1. Lo siento Alvaro, pero no tengo ningún programa interpolación. En el caso de que lo hiciera ya te comentaría por aquí.
      Gracias por la visita.

      Eliminar
  2. amigo puedes compartir un poco de la teoría que utilizaste para este programa. Gracias

    ResponderEliminar
  3. hola; me podrías ayudar con este algoritmo en pseint, lo tengo en grado 1 pero no se como hacerlo en grado 4

    ResponderEliminar
  4. Hola Adriana, pues la verdad es que noconozco nada a cerca de Pseint. Lo que sí te puedo decir es que si quieres implementar el método de orden cuatro, deberás ser capaza de resolver un sistema de orden cuatro. Te recomiendo si te interesa que te mires el programa del blog en el que se explica como resolver mediante Gauss-Jordan un sistema lineal: http://cypascal.blogspot.fr/2013/05/metodo-de-resolucion-de-ecuaciones-en-c.html
    La verdad que basta con que te mires la función resuleveGauss.
    Un saludo

    ResponderEliminar
  5. Hola, si me podrías explicar a que se refiere con orden del polinomio?, por favor!, me sería de gran ayuda, gracias.

    ResponderEliminar
    Respuestas
    1. Hola, se refiere al orden del polinomio con el que quieres aproximar los datos. Si eliges de orden dos, te aproximara los datos introducidos a una funcion del tipo a*x^2+b*x+c, de manera que se minimice el error cuadratico medio
      Un saludo

      Eliminar
  6. Amigo ayudame con algo. Estoy colocando el codigo en c++ pero cuando tiene que mostrar el resultado la consola se cierra muy rapido y no me deja verlo, ya probe con system(PAUSE); y no pasa nada. Espero que me puedas ayudar

    ResponderEliminar
  7. Gracias amigo me sirvió mucho... saludos

    ResponderEliminar
  8. no tendrás un ejemplo de un programa sobre el método de interpolación lineal?

    ResponderEliminar