|
Tratamiento de datos |
Supongamos que hemos medido un conjunto de pares de datos (xi, yi) en una experiencia, por ejemplo, la posición de un móvil en ciertos instantes de tiempo.
Queremos obtener una función y=f(x) que se ajuste lo mejor posible a los valores experimentales. Se pueden ensayar muchas funciones, rectas, polinomios, funciones potenciales o logarítmicas.
Una vez establecido la función a ajustar se determinan sus parámetros, en el caso de un polinomio, serán los coeficientes del polinomio de modo que los datos experimentales se desvíen lo menos posible de la fórmula empírica.
La función más sencilla es la función lineal y=ax+b que hemos tratado en la página anterior. El procedimiento de ajustar los datos experimentales a una línea recta se denomina regresión lineal
Queremos aproximar un polinomio de grado n, a un conjunto de m pares de datos (xi, yi) de modo que n£ m.
Sea el polinomio
P(x)=a0+a1x+a2x2+...anxn
Se calcula la cantidad
Para obtener los valores de los coeficientes del polinomio aproximador se tienen que determinar los valores de los coeficientes a0, a1, a2, ...an de forma que la cantidad S tome un valor mínimo.
Hagamos las derivadas parciales de S respecto de a0, a1, a2, ...an iguales a cero
(1)
Obtenemos un sistema de n+1 ecuaciones con n+1 incógnitas, a0, a1, a2, ...an
Ejemplo:
Supongamos que tenemos 4 pares de datos y que queremos ajustarlos al polinomio de segundo grado y=a0+a1x+a2x2
| x | x0 | x1 | x2 | x3 |
| y | y0 | y1 | y2 | y3 |
Las ecuaciones (1) se escribirán
agrupando términos
Volvamos al sistema de n+1 ecuaciones con n+1 incógnitas. Introduzcamos las expresiones
(2)
Se obtiene el siguiente sistema de n+1 ecuaciones con n+1 incógnitas
(3)
Si todos los puntos son distintos, el sistema de ecuaciones tiene una solución única.
Para resolver el sistema de ecuaciones se puede emplear alguno de los varios procedimientos existentes. El método empleado en este programa es el método de iteración. Ya que la matriz del sistema de ecuaciones es positiva, para este sistema se emplea el procedimiento de Seidel.
La clase PolRegresion tiene los siguientes miembros dato
public class PolRegresion {
private double[] x; //datos
private double[] y;
private int nDatos;
double[][] m; //matriz de los coeficientes
double[] t; //términos independientes
public double[] a; //polinomio a[0]+a[1]·x+a[2]·x2+...
public int grado; //grado del polinomio
El constructor inicializa los miembros, array x, array y que guardan los resultados experimentales y el número nDatos de dichos datos. El constructor crea la matriz m, el vector t, y reserva espacio para el array que va a guardar los coeficientes del polinomio a, cuyo grado se guarda en el miembro grado.
public PolRegresion(double[] x, double[] y, int grado) {
this.x=x;
this.y=y;
nDatos=x.length;
this.grado=grado;
t=new double[grado+1];
m=new double[grado+1][grado+1];
a=new double[grado+1];
}
A partir de datos (xi, yi) se tiene que obtener la matriz m cuadrada de los coeficientes (3) de dimensión n (grado+1) y el vector t de dimensión n usando las fórmulas (2). La función miembro coeficientes calcula los elementos del vector t y de la matriz m.
private void coeficientes(){
double[] s=new double[2*grado+1];
double suma;
for(int k=0; k<=2*grado; k++){
suma=0.0;
for(int i=0; i<nDatos; i++){
suma+=potencia(x[i], k);
}
s[k]=suma;
}
for(int k=0; k<=grado; k++){
suma=0.0;
for(int i=0; i<nDatos; i++){
suma+=potencia(x[i], k)*y[i];
}
t[k]=suma;
}
for(int i=0; i<=grado; i++){
for(int j=0; j<=grado; j++){
m[i][j]=s[i+j];
}
}
}
La función coeficientes llama a la función auxiliar potencia para calcular la potencia de un número elevado a un exponente entero
private double potencia(double base, int exp){
double producto=1.0;
for(int i=0; i<exp; i++){
producto*=base;
}
return producto;
}
Finalmente, nos queda la tarea de resolver el sistema de n ecuaciones con n incógnitas (n es el grado del polinomio). El procedimiento empleando es el de Seidel. Explicaremos este procedimiento mediante un ejemplo. Sea el sistema lineal de tres ecuaciones
Se despeja x1 en la primera ecuación, x2 en la segunda y x3 en la tercera. Reduciéndose el sistema a una forma conveniente para aplicar el procedimiento de iteración o de aproximaciones sucesivas.
(4)
Como podemos apreciar, las diagonales de la matriz de los coeficientes de las x son ceros.
Tomemos como aproximación inicial
La primera aproximación se escribe
Se introduce la primera aproximación de x0 en la segunda ecuación. Se introduce la primera aproximación de x0 y de x1 en la tercera ecuación, y así sucesivamente. De este modo el método de Seidel ofrece una mejor convergencia que la iteración simple.
La segunda aproximación es
En general el método se Seidel lo podemos escribir del siguiente modo.
En todo proceso iterativo, es preciso establecer la condición de finalización, no puede correr el proceso de forma indefinida. El error relativo entre dos iteraciones consecutivos k y k+1 de los valores que toman cada una de las incógnitas xi deberá ser menor que un cierto valor (por ejemplo, una milésima).
Veamos ahora el código
En primer lugar, transformamos la matriz m de los coeficientes del sistema de ecuaciones (3) y el vector t de los términos independientes, para expresarlo de la forma (4) adecuada para la iteración. Para ello, se divide los elementos de cada fila i por el elemento correspondiente de la diagonal principal m[i][i] y se cambia de signo. Los elementos de la diagonal m[i][i] se hacen cero.
for(int i=0; i<=grado; i++){
aux=m[i][i];
for(int j=0; j<=grado; j++){
m[i][j]=-m[i][j]/aux;
}
t[i]=t[i]/aux;
m[i][i]=0.0;
}
Se elige la aproximación inicial, y se guardan en el array p de las incógnitas
p[0]=t[0];
for(int i=1; i<=grado; i++){
p[i]=0.0;
}
Se calcula siguiente aproximación y se guarda en el array a de las incógnitas
for(int i=0; i<=grado; i++){
a[i]=t[i];
for(int j=0; j<i; j++){
a[i]+=m[i][j]*a[j];
}
for(int j=i+1; j<=grado; j++){
a[i]+=m[i][j]*p[j];
}
}
La iteración actual a es la iteración anterior p en el siguiente cálculo
for(int i=0; i<=grado; i++){
p[i]=a[i];
}
El proceso se repite, hasta que se cumple la condición de finalización.Para ello, calculamos el error relativo de cada una de las incógnitas. Guardamos en la variable maximo el mayor de los errores relativos.
maximo=0.0;
for(int i=0; i<=grado; i++){
error=Math.abs((a[i]-p[i])/a[i]);
if (error>maximo) maximo=error;
}
Todo el procedimiento lo ponemos en el interior de un bucle do...while que se ejecuta hasta que se cumple que el valor de la variable maximo sea menor o igual que una milésima.
double error=0.0, maximo=0.0;
do{
error=0.0;
maximo=0.0;
for(int i=0; i<=grado; i++){
a[i]=t[i];
for(int j=0; j<i; j++){
a[i]+=m[i][j]*a[j];
}
for(int j=i+1; j<=grado; j++){
a[i]+=m[i][j]*p[j];
}
error=Math.abs((a[i]-p[i])/a[i]);
if (error>maximo) maximo=error;
}
for(int i=0; i<=grado; i++){
p[i]=a[i];
}
}while(maximo>0.001);
El código completo de la clase PolRegresion es el siguiente
public class PolRegresion {
private double[] x; //datos
private double[] y;
private int nDatos;
double[][] m; //matriz de los coeficientes
double[] t; //términos independientes
public double[] a; //polinomio a[0]+a[1]·x+a[2]·x2+...
public int grado; //grado del polinomio
public PolRegresion(double[] x, double[] y, int grado) {
this.x=x;
this.y=y;
nDatos=x.length;
this.grado=grado;
t=new double[grado+1];
m=new double[grado+1][grado+1];
a=new double[grado+1];
}
private void coeficientes(){
double[] s=new double[2*grado+1];
double suma;
for(int k=0; k<=2*grado; k++){
suma=0.0;
for(int i=0; i<nDatos; i++){
suma+=potencia(x[i], k);
}
s[k]=suma;
}
for(int k=0; k<=grado; k++){
suma=0.0;
for(int i=0; i<nDatos; i++){
suma+=potencia(x[i], k)*y[i];
}
t[k]=suma;
}
for(int i=0; i<=grado; i++){
for(int j=0; j<=grado; j++){
m[i][j]=s[i+j];
}
}
}
private double potencia(double base, int exp){
double producto=1.0;
for(int i=0; i<exp; i++){
producto*=base;
}
return producto;
}
//procedimiento de Seidel
public void calculaPolinomio(){
coeficientes();
//matriz
double aux;
for(int i=0; i<=grado; i++){
aux=m[i][i];
for(int j=0; j<=grado; j++){
m[i][j]=-m[i][j]/aux;
}
t[i]=t[i]/aux;
m[i][i]=0.0;
}
// aproximación inicial
double[] p=new double[grado+1];
p[0]=t[0];
for(int i=1; i<=grado; i++){
p[i]=0.0;
}
//aproximaciones sucesivas
double error=0.0, maximo=0.0;
do{
maximo=0.0;
for(int i=0; i<=grado; i++){
a[i]=t[i];
for(int j=0; j<i; j++){
a[i]+=m[i][j]*a[j];
}
for(int j=i+1; j<=grado; j++){
a[i]+=m[i][j]*p[j];
}
error=Math.abs((a[i]-p[i])/a[i]);
if (error>maximo) maximo=error;
}
for(int i=0; i<=grado; i++){
p[i]=a[i];
}
}while(maximo>0.001);
}
}
|
Determinar el polinomio aproximador de segundo grado para la siguientes tablas de datos, tomadas de una experiencia con fluidos
| x | 0.05 | 0.1 | 0.15 | 0.2 | 0.25 | 0.3 | 0.35 | 0.4 |
| y | 0.375 | 0.625 | 0.750 | 0.815 | 0.875 | 1.00 | 1.065 | 1.125 |
| x | 0.05 | 0.1 | 0.15 | 0.2 | 0.25 | 0.3 | 0.35 | 0.4 |
| y | 0.25 | 0.42 | 0.58 | 0.72 | 0.85 | 0.98 | 1.10 | 1.12 |
| x | 0.05 | 0.1 | 0.15 | 0.2 | 0.25 | 0.3 | 0.35 | 0.4 |
| y | 0.28 | 0.50 | 0.64 | 0.75 | 0.82 | 0.90 | 0.95 | 1.02 |
El siguiente applet permite introducir los datos experimentales y ajustarlos mediante un polinomio de grado menor o igual al número de datos. Los datos que aparecen cuando se inicia el applet en los controles de edición corresponden a la primera tabla.
Se elige el grado del polinomio pulsando en el botón titulado Grado polinomio>> y a continuación se pulsa en el botón titulado Calcular.
Se observa a la derecha la representación de los datos, del polinomio de ajuste y los valores de los coeficientes del polinomio.
P(x)=a0+a1x+a2x2+...anxn
Para ensayar con otros datos se pulsa el botón titulado Borrar, y a continuación se introducen los pares de datos en los controles de edición respectivos.
PolRegresion.java, RegresionApplet2.java, MiCanvas.java
B.P. Demidowitsch, Maron, Schuwalowa. Métodos numéricos de análisis. Editorial Paraninfo (1978), págs 21-23.
Demidovich, Maron. Cálculo numérico fundamental. Editorial Paraninfo (1977), págs 339-340.