MANUAL DE REFERENCIA
 
 
Ejemplos Hechos por Vicente Fuentes Gea
2007
 
Indice
 

Banderas
Interpolación
Regresión lineal
Runge-Kutta

Sistema de ecuaciones algebraicas


 Runge-Kutta
 
//Mediante el método de Runge-Kutta de orden 4, este programa integra un
//conjunto de ecuaciones diferenciales y grafica los resultados.
//El programa utiliza solo arreglos unidimensionales porque el intérprete
//no conoce otro tipo de arreglos. El acceso a los elementos A[i][j] de un
//arreglo bidimensional declarado así: float A[X][Y], puede simularse
//mediante un arreglo unidimensional float B[X * Y], así: B[i * Y + j]. En
//donde i = 0... X - 1 y j = 0... Y - 1. En el caso de un arreglo
//tridimensional declarado así: float A[X][Y][Z], el acceso a los elementos
//A[i][j][k] se simula mediante el arreglo B[X * Y * Z] asi:
//B[i * Y * Z + j * Z + k]. En done  i = 0.. X - 1, j = 0... Y – 1 y
//k = 0... Z – 1
 
void Derivadas(float f[], float t, float y[])
{
 // Aquí se definen las ecuaciones diferenciales
 // Regresa las derivadas dy/dt = f(t,y) en f[]
 // Recibe en y[] las variables dependientes. En t la variable independiente.
 
 f[0] = -y[1];
 f[1] = y[0] - (1 / t) * y[1];
 f[2] = y[1] - (2 / t) * y[2];
 f[3] = y[2] - (3 / t) * y[3];
 
}
 
main()
{
 int i;
 int const np = 50, ne  = 4; // No de intervalos de integración y ecuaciones diferenciales
 float dt = 0.2, t0 = 1;     // Intervalo de integración y tiempo inicial
 float y[ne * np], t[np];  // Reciben los valores de las variables dependientes e independiente
 float y0[ne];               // Condiciones iniciales
 int color[ne] = {9, 10, 11, 12};
    
 for(i = 0; i < ne; i++)
   y0[i]  = besjn(i, 1);          // Asignación de condiciones iniciales usando funciones de Bessel
 
 RungeKuta4(t, y, np, ne, dt, t0, y0);
 
 ImprimeValores(t, y, np, ne);
 
 
 GraficaEcua(t, y, np, ne, color);
 
}
 
void RungeKuta4(float t[], float y[], int nite, int necua, float h, float t0, float y0[])
{
  //Regresa en los arreglos t[nite], y[necua * nite] los valores
  //calculados de las variable independiente y dependientes, respectivamente.
  //nite      = No de intervalos de integración
  //necua     = No de ecuaciones diferenciales
  //h         = Intervalo de integración
  //Derivadas = La función que entrega las derivadas en el arreglo f[] y recibe,
  //en ye[] las variables dependientes conteniendo los valores estimados
 
 int i, j;
 float h2 = h / 2;
 float f[necua], ye[necua], L1[necua], L2[necua], L3[necua], L4[necua];
 
 for( i = 0;i < necua;i++)            //  Asignacion de Condiciones iniciales
  ye[i] = y[i * nite] = y0[i];        //  ye --> valores estimados 
 
 t[0] = t0;
 
 for(j = 0; j < nite - 1; j++)
 {
  Derivadas(f, t[j], ye);             // paso 1 
  for(i = 0; i < necua; i++)
  {
   L1[i] = h * f[i];
   ye[i] = y[i * nite + j] + L1[i] / 2;
  }
 
  Derivadas(f, t[j] + h2, ye);        //  paso 2 
  for(i = 0; i < necua; i++)
  {
   L2[i] = h * f[i];
   ye[i] = y[i * nite + j] + L2[i] / 2;
  }
  
  Derivadas(f, t[j] + h2, ye);        //  paso 3 
  for(i = 0;i < necua; i++)
  {
   L3[i] = h * f[i];
   ye[i] = y[i * nite + j] + L3[i];
  }
  
  Derivadas(f, t[j] + h, ye);         //  paso 4 
  for(i = 0;i < necua; i++)
   L4[i] = h * f[i];
  
  for(i = 0;i < necua; i++)
    ye[i] = y[i * nite + j + 1] = y[i * nite + j] + (L1[i] + 2 * L2[i] + 2 * L3[i] + L4[i]) /6;
  
  t[j + 1] = t[j] + h;
   }
}
 
float Vmax(float a[], int n)
{
 int i;
 float ma = a[0];
 
 for(i = 1;i < n;i++)
  ma = max(ma, a[i]);
      
 return ma;
}
 
float Vmin(float a[], int n)
{
 int i;
 float mi = a[0];
 
 for(i = 1;i < n;i++)
  mi = min(mi, a[i]);
 
 return mi;
}
 
void GraficaEcua(float x[], float y[], int np, int ne, int color[])
{
 int i, j;
 float xmin, ymin, xmax, ymax;
 
 xmin = Vmin(x, np);
 ymin = Vmin(y, np * ne);
 xmax = Vmax(x, np);
 ymax = Vmax(y, np * ne);
 
escala(xmin, ymin, xmax, ymax);
for(i = 0; i < ne; i++)
 {
        setcolor(color[i]);
        moveto(x[0], y[i * np]);
        for(j = 1; j < np;j++)
               lineto(x[j], y[i * np + j]);
 }
}
 
void ImprimeValores(float x[], float y[], int np, int ne)
{
 int i, j;
 
 for(i = 0; i < ne; i++)
 {
  printf("\n\t\t    Ecuacion %d\n",i+1);
  printf("\n       t\t      y\t     solucion exacta \n");
  printf("   ---------------------------------------------\n");
  for(j = 0; j < np;j++)
   printf("   %10.4f  %10.4f  %10.4f\n", x[j], y[i * np + j], besjn(i, x[j]));
 }
}

Banderas
//      Función Ban. Permite manipular 32 banderas usando un número entero. Cada bit del
//      número entero se constituye en una bandera, permitiendo prescindir de un arreglo
 
 
const int TODAS = 0;
const int APAGA = 0;
const int PRENDE = 1;
const int PRUEBA = 2;
 
int banderas;  // Esta variable contiene las banderas, pero como
               // el interprete desconoce las variables locales estáticas
               // la tenemos que hacer global, lástima.
 
 
int Ban(int num, int comando)
{
 if(num > 32 ) return -2;      // solo 32 banderas, de la 1....32
 
 if(comando == PRUEBA) return (banderas & (1 << (num - 1))); // esta es una prueba
 
 if(!num && comando == PRENDE) // prende todas
 {
   banderas =  4294967295;     // UINT_MAX en binario 11111111111111111111111111111111
   return -1;
 }
 
 if(!num && comando == APAGA)  // apaga todas
 {
   banderas = 0;
   return -1;
 }
//                     la prende               o la         apaga según comando
 banderas = comando ? banderas | (1 << (num - 1)) : banderas & ~(1 << (num - 1));
 
 return -1;
}
 
int main()
{
 int i;
 
 Ban(TODAS, APAGA); // todas las banderas apagadas
 
 puts("Prendemos las banderas una a una\n");
 
 for(i = 1; i < 33; i++)
   Ban(i, PRENDE);
 
// Vamos a ver si es cierto
 for(i = 1; i < 33; i++)
 {
   if(Ban(i, PRUEBA))
     printf("bandera %d: prendida\n", i);
 }
 
 puts("\nAhora las apagamos también una a una\n");
 
 for(i = 1; i < 33; i++)
   Ban(i, APAGA);
 
// Vamos a ver si es cierto
 for(i = 1; i < 33; i++)
 {
   if(!Ban(i, PRUEBA))
     printf("bandera %d: apagada\n", i);
 }
 
 puts("\nAhora las prendemos al azar\n");
 
 srand(clock()); // sembramos un numero aleatorio con clock()
                 // clock() comenzó a contar desde que se corrió este interprete
 
 for(i = 1; i < 33; i++)
   Ban(i, rand() % 2);
 
 for(i = 1; i < 33; i++)
 {
   if(Ban(i, PRUEBA))
     printf("bandera %d: prendida\n", i);
   else
     printf("bandera %d: apagada\n", i);
 }
 
 Ban(TODAS, PRENDE);
 
 puts("\nAhora todas prendidas\n");
 for(i = 1; i < 33; i++)
 {
   if(Ban(i, PRUEBA))
     printf("bandera %d: prendida\n", i);
   else
     printf("bandera %d: apagada\n", i);
 }
 
 return 0;
}
 


Regresión lineal
 // Este programa efectúa una regresión lineal por el método de los mínimos cuadrados 
 
int main()
{
    float mb[2], r;
    float minx, miny, maxx, maxy;
    char tmp[100];
    const int n = 25;
    // datos utilizados en el programa :
    float x[n] = {450,545,518,487,477,474,479,520,497,479,475,473,478,
                  509,487,475,472,470,477,568,542,505,485,477,481};
    float y[n] = {437,598,520,493,482,479,477,507,505,498,483,480,478,
                  497,500,490,482,478,476,585,517,492,482,478,477};
 
    r = RegLin(x, y, n, mb);
    
    settextstyle(0, HORIZ_DIR, 1.5);
    sprintf(tmp, "Coeficiente de correlación = %2.2f", r);
    setcolor(WHITE);
    outtextxy(getmaxx() / 4 , 20, tmp);
    minx = Vmin(x, n);
    miny = Vmin(y, n);
    maxx = Vmax(x, n);
    maxy = Vmax(y, n);
    
    escala(minx, miny, maxx, maxy);
    
    setcolor(LIGHTGREEN);
    line(0, mb[1], maxx, mb[0] * maxx + mb[1]);
    GraficaCirculos(x, y, n, (maxy - miny) / 90, LIGHTRED);
    
    escala(); // destruimos la escala
    sprintf(tmp, "y = %2.2f*x%+2.2f", mb[0], mb[1]);
    outtextxy(getmaxx() / 3 , 70, tmp);
    
    return 0;
}
 
 
float RegLin(float x[], float y[], int n, float mb[])
{
// Ajusta por mínimos cuadrados los datos contenidos en los arreglos x[], y[]
// a la recta y = mx + b, retornando el coeficiente de correlación r de la
// manera usual. Los valores de m, b se regresan usando el arreglo mb[]
// (lástima, pero el intérprete no cuenta con apuntadores).
 
    int i ;
    float xy, x2, y2, sx, sy, sx2, sy2, sxy ;
    float r ; 
    
    sx = sy = sx2 = sy2 = sxy = 0 ;
    
    for(i = 0; i < n; i++)
    {
        sx  += x[i] ;
        sy  += y[i] ;
        sx2 += x[i] * x[i] ;
        sy2 += y[i] * y[i] ;
        sxy += x[i] * y[i] ;
    }
    
    xy = sxy - sx * sy / n ;
    x2 = sx2 - sx * sx / n ;
    y2 = sy2 - sy * sy / n ;
    
    mb[0] = xy / x2 ;                    // pendiente 
    mb[1] = sy / n - mb[0] * sx / n ;    // ordenada al origen 
    r = sqrt(xy * xy / (x2 * y2)) ;      // coeficiente de correlación
    return r;
}
 
float Vmax(float a[], int n)
{
    int i;
    float ma = a[0];
    
    for(i = 1;i < n;i++)
        ma = max(ma, a[i]);
    
    return ma;
}
 
float Vmin(float a[],int n)
{
    int i;
    float mi = a[0];
    
    for(i = 1;i < n;i++)
        mi = min(mi, a[i]);
    
    return mi;
}
 
void GraficaCirculos(float x[], float y[], int n, float r, int c)
{ 
    int i; 
    int ca = getcolor();
    
    setcolor(15);
    setfillstyle(0, c);
    for(i = 0; i < n;i++)
    {
        fillellipse(x[i], y[i], r, r);  
        circle(x[i], y[i], r);
    }
    
    setcolor(ca);
}
 


Sistema de ecuaciones algebraicas
// Resuelve n ecuaciones algebraicas lineales usando el método de
//              Gauss - Jordan con pivoteo parcial
 
 
#include <math.h>
#include <stdio.h>
 
 
//#define ne 14
const int ne = 14;
 
void ImprimeEcuaciones(double coef[], int filas);
void ImprimeIncognitas(double coef[], int filas);
int ResuelveEcuaciones(double coef[], int filas, double tolerancia);
 
 
int main()
{
    double coef[ne * (ne + 1)] =
    {
        1,    1,    1,  0,   0,  0,   0,   0,   0,   0,  0,   0,   0,   0,   43.93,
        1.17, 0,    0,  -1,  0,  0,   0,   0,   0,   0,  0,   0,   0,   0,   0,
        0,    0,    0,  0,   1,  0,   0,   0,   0,   0,  0,   0,   0,   0,   95.798,
        0,    0,    1,  0,   1, -1,  -1,  -8,   0,   0,  0,   0,   1,   0,   99.1,
        0,    0,    0,  0,   0,  1,   1,   1,   1,  -1, -1,   0,   0,   0,   -8.4,
        0,    0,    0,  1,   0,  0,   0,   0,   0,   0,  0,   0,  -1,   0,   24.2,
        1,    0,    0, -1,   0,  0,   0,   0,   0,   1,  0,   0,   0,   1,   189.14,
        0,    0,    0,  0,   0,  0,   0,   0,   0,4.594, 0,   0,   0, .11,   146.55,
        0,    0,    0,  0,   0,  0,   0,   0,   1,   0,  0,   0,   0,   0,   10.56,
        0,    1,    0,  0,   0,  0,   0,   0,   0,   0,  0,   0,   0,   0,   2.9056,
        0,    0,    0,  0,   0,  1,   0,   0,   0,   0,  0,   0,   0,-.0147, 0,
        0,    0,    1,  0,   0,  0,   0,   0,   0,   0,  0, -.07,  0,  0,    0,
        0,    0,    0,  0,   0,  0,   1,   0,   0,   0,  0,   0,   0,  0,  14.618,
        0,    0,    0,  0,   0,  0,   0,   0,   0,   1,  0,  -1,   0,  1,  -97.9
        
    };
    
    ImprimeEcuaciones(coef, ne);
    
    if(ResuelveEcuaciones(coef, ne, 1e-10))
        ImprimeIncognitas(coef, ne);
    else
        puts("\n\aMatriz singular");
 
    return 0;
}
 
int ResuelveEcuaciones(double coef[], int filas, double tolerancia)
{
    int i, j, k, l, m;
    int col = filas + 1;
    double temp;
    
    for (i = 0; i < filas; i++)
    {
        // pivoteo parcial, coef[i][i] = pivote 
        for(m = i; m < filas; m++)
        {
            if(fabs(coef[i * col + i]) < fabs(coef[m * col + i]))
            {
                for(l = 0; l < filas; l++)
                {
                    temp = coef[m * col + l];
                    coef[m * col + l] = coef[i * col + l];
                    coef[i * col + l] =  temp;
                }
                temp = coef[m * col + col - 1];
                coef[m * col + col - 1] = coef[i * col + col - 1];
                coef[i * col + col - 1] = temp;
            }
        }
        
        if(fabs(coef[i * col + i]) < tolerancia) return 0;       // matriz singular 
        
        for (j = col - 1; j >= 0; j--)
            coef[i * col + j] /= coef[i * col + i];
        
        for (j = i + 1; j < filas; j++)
            for (k = col - 1;k >= i;k--)
                coef[j * col + k] -=  coef[j * col + i]  *  coef[i * col + k];
    }
    
    for (i = filas - 2; i >= 0; i--)
        for (j = col - 2; j > i; j--)
        {
            coef[i * col + col - 1] -= coef[i * col + j] * coef[j * col + col - 1];
            coef[i * col + j] = 0;
        }
        
    return 1;
}
 
void ImprimeIncognitas(double coef[], int filas)
{
    int i, j;
    int col = filas + 1;
    
    puts("\nSolucion del sistema de ecuaciones algebraicas:");
    
    for (i = 0;i < filas;i++)
    {
        printf("\n");
        for(j = 0;j < filas; j++)
            printf("%5.1f",coef[i * col + j]);
        
        printf("      x%d    = %6.4g",i + 1, coef[i * col + j]);
    }
}
 
void ImprimeEcuaciones(double coef[], int filas)
{
    int i, j;
    int col = filas + 1;
    
    puts("Sistema de ecuaciones algebraicas:\n");
    
    for (i = 0;i < filas;i++)
    {
        for(j = 0;j < filas - 1; j++)
            printf("%8.4fx%d +", coef[i * col + j],j + 1);
 
        printf("%8.4fx%d ", coef[i * col + j], j + 1 );
        printf("    = %8.2f \n", coef[i * col + j + 1]);
    }
}


Interpolación
// Interpola una serie de puntos mediante la técnica de interpolación coseno y de Hermite
 
// Correr el programa varias veces para obtener curvas diferentes
 
int main()
{
    int n0, n;
    double xx;
    double x[100], y[100], xi[1000], yi[1000];
 
    srand(clock()); // queremos generar curvas diferentes en cada corrida
 
    setbkcolor(WHITE);
    cleardevice();
 
    setcolor(BLACK);
    settextstyle( 0, HORIZ_DIR, 2);
    outtextxy(getmaxx() / 20, getmaxy() - 75, "Interpolacion coseno");
    outtextxy(3 * getmaxx() / 6, getmaxy() - 75, "Interpolacion de Hermite");
 
    n0 = 0;
    for(xx = 0;xx <= 115;xx += 15) 
    {
        x[n0] = xx;
        y[n0] = rand() % 116;
        n0++;
    }
 
    setcolor(DARKGRAY);
    rectangle(25, 100, getmaxx() / 2 - 25, getmaxy() - 100);
    setviewport(25, 100, getmaxx() / 2 - 25, getmaxy() - 100, 1);
 
    n = CreaCurva2(x, y, n0, xi, yi);  // interpolación coseno
    Grafica(x, y, n0, LIGHTGRAY, 1);   // los puntos sin alizar
    Grafica(xi, yi, n, LIGHTRED, 0);   // la curva ya alizada
 
    escala();   // eliminamos la escala y la región de la gráfica anterior 
    setviewport();  
 
    // creamos una nueva región:              
    rectangle(getmaxx() / 2 + 25, 100, getmaxx() - 25, getmaxy() - 100);
    setviewport(getmaxx() / 2 + 25, 100, getmaxx() - 25, getmaxy() - 100, 1);
    
    n = CreaCurva4(x, y, n0, xi, yi);  // interpolación de Hermite
    Grafica(x, y, n0, LIGHTGRAY, 1);   // los puntos sin alizar
    Grafica(xi, yi, n, LIGHTBLUE, 0);  // la curva ya alizada
 
    return 0;
}
 
void Interpola2(double x1, double y1, double x2, double y2, double xi[], double yi[], int n, int i0)
{
    int i;
    double m, dm, x, dx;
 
    dm = 1. / n;
    dx = (x2 - x1) / n;
 
    i = i0;
    m = 0;
    for(x = x1; x <= x2; x += dx)
    {
        xi[i] = x;
        yi[i] = IntCoseno(y1, y2, m);
        i++;
        m += dm;
    }
}
 
int CreaCurva2(const double x[], const double y[], const int Np, double xi[], double yi[])
{
    // llena xi[], yi[] con puntos interpolados usando la interpolación coseno
    // y regresa el número de puntos interpolados. 
    const int nsi = 10; // dividir los intervalos en 10 subintervalos
    int i;
    int Nint = Np * nsi, i0 = 0;
 
    for(i = 0 ;i < Np - 1;i++)
    {
        Interpola2(x[i], y[i], x[i + 1], y[i + 1], xi, yi, nsi, i0);
        i0 += nsi;
    }
 
    return Nint - nsi + 1;
}
 
void Interpola4(double y1, double x2, double y2, double x3, double y3, double y4,
                    double xi[], double yi[], int n, int i0)
{
    int i;
    double m, dm, x, dx;
 
    dm = 1. / n;
    dx = (x3 - x2) / n;
 
    i = i0;
    m = 0;
    for(x = x2; x <= x3; x += dx)
    {
        xi[i] = x;
        yi[i] = IntHermite(y1, y2, y3, y4, m, .5, 0);
        i++;
        m += dm;
    }
}
 
int CreaCurva4(const double x[], const double y[], const int Np, double xi[], double yi[])
{
    // llena xi[], yi[] con puntos interpolados usando la interpolación de Hermite
    // y regresa el número de puntos interpolados.
    const int nsi = 10; // dividir los intervalos en 10 subintervalos
    int i;
    int Nint = Np * nsi, i0;
 
    Interpola2(x[0], y[0], x[1], y[1], xi, yi, nsi, 0); // el primer segmento con el coseno
 
    i0 = nsi;
    for(i = 1;i < Np - 2;i++)
    {
        Interpola4(y[i - 1], x[i], y[i], x[i + 1], y[i + 1], y[i + 2], xi, yi, nsi, i0);
    i0 += nsi;
    }
                                    //el último segmento también con el coseno
    Interpola2(x[Np - 2], y[Np - 2], x[Np - 1], y[Np - 1], xi, yi, nsi, i0); 
 
    return Nint - nsi + 1;
}
 
double IntCoseno(double y1, double y2, double mu)
{
    // 0 <= mu <= 1;   mu = 0 en y1, mu = 1 en y2. Regresa el valor de y en mu
   double mu2;
 
   mu2 = (1. - cos(mu * PI)) / 2.;
   return (y1 * (1. - mu2) + y2 * mu2);
}
 
double IntHermite(double y0, double y1, double y2, double y3, double mu,
                    double tension, double sesgo)
{
    // 0 <= mu <= 1;   mu = 0 en y1, mu = 1 en y2. Regresa el valor de y en mu
    // Tensión: 1 alta, 0 normal, -1 baja
    // Sesgo: 0 es neutral, + hacia el primer segmento, - hacia el otro
 
    double m0, m1, mu2, mu3;
    double a0, a1, a2, a3;
    double ft = (1. - tension) / 2.;
    
    mu2 = mu * mu;
    mu3 = mu2 * mu;
    m0  = (y1 - y0) * (1. + sesgo) * ft;
    m0 += (y2 - y1) * (1. - sesgo) * ft;
    m1  = (y2 - y1) * (1. + sesgo) * ft;
    m1 += (y3 - y2) * (1. - sesgo) * ft;
    a0 = 2. * mu3 - 3. * mu2 + 1.;
    a1 = mu3 - 2. * mu2 + mu;
    a2 = mu3 -  mu2;
    a3 = -2. * mu3 + 3. * mu2;
 
    return(a0 * y1 + a1 * m0 + a2 * m1 + a3 * y2);
}
 
double Vmax(double a[], int n)
{
    int i;
    double ma = a[0];
 
    for(i = 1;i < n;i++)
        ma = max(ma, a[i]);
      
    return ma;
}
 
double Vmin(double a[],int n)
{
    int i;
    double mi = a[0];
 
    for(i = 1;i < n;i++)
        mi = min(mi, a[i]);
 
    return mi;
}
 
void Grafica(double x[], double y[], int n, int c, int ban)
{
    int i;
    double xmin, ymin, xmax, ymax;
    int ca = getcolor();
 
    if (ban)       //si la bandera esta puesta calcula la escala
    {
        xmin = Vmin(x, n);
        ymin = Vmin(y, n);
        xmax = Vmax(x, n);
        ymax = Vmax(y, n);
        escala(xmin, ymin, xmax, ymax);
    }
 
    setcolor(c);
    moveto(x[0], y[0]);
    for(i = 1; i < n;i++)
            lineto(x[i], y[i]);
 
    setcolor(ca);
}