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);
}