Departamento de Matem´ atica atica Aplicada
´ CALCULO COMPUTACIONAL. COMPUTACIONAL. Licenciatura en Qu´ Qu´ımica (Curso 2005-06) M´ınimos ıni mos Cuadrad Cua drados. os. Pr´ actica acti ca 7
1.
Intr Introdu oducc cci´ i´ on on
El M´ etodo etodo de ajuste por M´ınimos Cuadrados sirve para encontrar una funci´on on y = f ( f (x, α1 , α2 , . . . , αm), en la que habr´a que calcular los par´ametros ametros α1 , α2 , . . . , α m . Esta funci´on on debe ser la que se ajuste lo mejor posible a una tabla de valores que relaciona relaciona las dos variables variables x e y obtenida experimentalmente: x x1 y y1
x2 y2
. .. . . xn . . . yn
Para calcular los par´ametros ametros se impone la condici´on on de que sea s ea m´ınima ınima la funci´on on n
S (α1 , α2 , . . . , αm ) =
[f ( f (xi , α1 , α2 , . . . , αm ) − yi ]2
i=1
Como S (α1 , α2 , . . . , αm ) es una funci´on on de m variables, una condici´on on necesaria para que tenga un valor extremo en un punto punto es que sus derivadas derivadas parciales en ese punto sean todas nulas. De aqu´ aqu´ı obtenemos obtenemos un sistema sistema de m ecuaciones, con m inc´ ognitas: ognitas: ∂S = 0, 0, ∂α 1
∂S = 0, 0, ∂α 2
. .. .. ,
∂S =0 ∂α m
cuyas soluciones, los par´ametros ametros α1 , α2 , . . . , α m nos indican c´omo omo es la funci´on on que mejor se ajusta a los datos, es decir, f ( f (x, α1 , α2 , . . . , αm ). La funci´ on on f ( f (x, α1 , α2 , . . . , αm) puede ser de cualquier tipo, te´oricamente, oricamente, sin embargo, en la pr´actica, actica, con programas como MATLAB (o DERIVE) s´olo olo se pueden calcular (directamente) cuando la funci´on on es un polinomio. En otros casos, que veremos en los ejercicios, habr´a que modificar ligeramente el problema para que se pueda tratar con el ordenador. Para realizar estos c´alculos alculos MATLAB dispone del comando polyfit. Veamos en un ejemplo c´omo omo se puede ajustar una recta con M´ınimos Cuadrados:
2.
Ejem jemplo plo Dada la tabla de valores, x 1 2 3 4 y 2,1 4,3 6 7,8 Encontrar la funci´on on de la forma y = ax + b que mejor se ajuste a los datos. Como se trata de una funci´on on polin´ omica se puede hacer directamente. omica Introducimos primero la tabla de valores en dos variables: >>x=[1 2 3 4] >>y=[2 >>y=[2.1 .1 4.3 6 7.8] 7.8]
El comando a utilizar es polyfit(x,y,n), donde n es el grado del polinomio que queremos obtener. Por lo tanto para obtener una recta n = 1: >>c=polyfit(x,y,1)
Que nos da como resultado los coeficientes de la recta: c = 1.8800
0.3500
25
Es decir, que la recta que hemos encontrado es, y = 1,88x + 0,35 Para representar la informaci´on obtenida gr´aficamente: Primero dibujamos la tabla de valores, por ejemplo: >>plot(x,y,’*’)
De esta forma conseguimos que dibuje s´olo los puntos, con asteriscos o con cualquier otro formato. Para dibujar la recta, lo hacemos como para dibujar cualquier funci´on. Generamos una tabla, (que llamaremos con un nombre diferente de x para no borrar la tabla de los datos del problema): >>xp=linspace(1,4,20);
Para para calcular los valores de xp en la recta y = 1,88x + 0,35 podemos utilizar el comando polyval que eval´ ua el polinomio utilizando los coeficientes, que ten´ıamos en la variable c: >>yp=polyval(c,xp); >>hold on % para mantener el dibujo anterior
Y, por u ´ ltimo, plot(xp,yp)
Y nos aparece el dibujo de la figura 2. 8
7
6
5 y=1.88x+0.35
4
3
2
1
1.5
2
2.5
3
3.5
4
Figura 2: Tabla y Recta y = 1,88x + 0,35 El error que se comete al aproximar la funci´on emp´ırica (tabla inicial) por la funci´ on te´orica (recta) se puede cuantificar de varias formas. Una manera es precisamente calcular la suma de las desviaciones en cada punto de la tabla al cuadrado, es decir, el valor de la funci´on: S (α1 , α2 , . . . , αm ): Primero calculamos los valores de f (xi , α1 , α2 , . . . , α m ): >>fxi=polyval(c,x) fxi = 2.2300
4.1100
5.9900
7.8700
Y ahora sustituimos en la f´ormula: >>error=sum((fxi-y).^2) error = 0.0580
26
Otra forma consiste en calcular lo que se llama Error Cuadr´ atico Medio, que viene dado por la expresi´on:
1 [f (x , α , α , . . . , α ε= n
n
i
1
2
2
m
) − yi ]
i=1
Para calcularlo: >>n=length(x) % calcula el n’umero de t’erminos n = 4 >>errorcm=sqrt(sum((fxi-y).^2)/n) errorcm = 0.1204
Por u ´ ltimo, se pueden utilizar los c´odigos que llevan incorporados los comandos polyfit y polyval para obtener un intervalo de confianza al 95 % en el que se encuentra cada yi. Veamos en qu´e consiste: Primero volvemos a utilizar el comando polyfit, de una forma diferente: >>[c,s]=polyfit(x,y,1) c = 1.8800 0.3500 s = R: [2x2 double] df: 2 normr: 0.2408 s es la estructura con la que se va a calcular el error. Ahora volvemos a utilizar polyval, pero tambi´en
modificado: >>[fxi,delta]=polyval(c,x,s) fxi = 2.2300 4.1100 5.9900 delta = 0.2220 0.1942 0.1942
7.8700 0.2220
Hemos obtenido, en cada xi , un intervalo de confianza al 95 %, [fxi-2*delta,fxi+2*delta], dentro del cual se encuentra el valor yi. Si representamos la gr´afica siguiente se entender´an de manera gr´afica los c´alculos que acabamos de realizar. >>plot(x,y,’*’,x,fxi,’g-’,x,fxi+2*delta,’r:’,x,fxi-2*delta,’r:’)
3.
Ejercicios
Ejercicio 1 Mediante el M´ etodo de M´ınimos Cuadrados, elegir una funci´ on cuadr´ atica del tipo f (x) = ax2 + bx + c para los valores de x e y dados por la siguiente tabla x y
7 8 9 10 11 12 13 7,4 8,4 9,1 9,4 9,5 9,5 9,4
Calcular la funci´ on f (x). El error, el error cuadr´ atico medio y representar la funci´ on obtenida junto a los datos de la tabla con una banda alrededor de la gr´ afica de f (x), dentro de la cual se encuentren los puntos de la tabla a un nivel del 95 %.
27
