ALA C, 2006-07, vecka 2, Studio 2b

[Kurssidan, vecka 1, vecka 2, Linjär algebra i R^n: Minsta kvadrat-metoden. 1.  Least squares methods. Define the matrix B and vector u as in the previous exercise (Studio 1). We recall that the vector u does not belong to R(B), so that the equation Bx=u has no solution.   It is therefore impossible to make the residual Bx-u equal to 0, and we try instead to minimize the norm of the residual: find x in R2 such that  || Bx-u ||2 is minimized. This is called "the least squares method".  We know that the distance  || Bx-u ||  is minimal when Bx is the projection of u onto V=R(B)=span(a,b), which is exactly what we computed in exercise 1 of Studio 1a.  So we have already computed the least squares solution x.    Now we will do it again in a different way. Pu is determined by the equation (u-Pu,v)=0  for all v in V=R(B) Here we have Pu=Bx and v=By for some x,y in R2.  Hence, we want to find x in R2  such that (u-Bx,By)=0  for all y in R2. Here  (u-Bx,By)=(BTu-BTBx,y)  so that (BTu-BTBx,y)=0 for all y in R2.  This implies that (*)   BTBx=BTu           (the normal equations) Here we used the fact that (w,y)=0 for all y, implies w=0.  Prove this!   Hint:  we can take y=w because the equality should hold for all y. Prove that BTB is symmetric. Set up the normal equations with the given matrices B, u and solve the equations using Matlab's backslash (x = (BTB)\(BTu)). Compare with what you got in Studio 1a. What happens when you do only x=B\u? In Matlab the backslash operator:  x=B\u  gives the least squares solution if the equation Bx=u has no solution. 2. Example Suppose that the variables y and x are related by y=kx+m. In order to determine the coefficients k and m we make measurements of y and x: Table 1 x 5 6 7 8 9 10 y 19.5888 23.4043 25.5754 29.1231 31.9575 35.8116 This leads to an overdetermined system of the form kx1+m=y1 . . . kx6+m=y6 or, in matrix form Av=y, where the unknown vector is v=[k, m]T. Solve this system by the least squares method in Matlab. Hint: set up the column vectors x, y and the matrix A=[x ones(size(x))], then form the matrices B=A'*A and g=A'*y. Solve the system Bv=g by the command v=B\g. Plot the data points (xi,yi) and the fitted function y=kx+m in the same figure. The following commands are useful: plot(x,y,'or'), fplot('ykxm',[x(1) x(6)]). Here ykxm.m is a function file that implements the function y=kx+m. Don't forget to declare global k m both inside the function file and in the main program. Actually, Matlab's backslash command v=A\y automatically uses the least squares method when the system Av=y is overdetermined. Try this also! 3.  The tank reactor.(Advanced.) The Arrhenius rate law is  k=k0 exp(-E/(RT)).   The task is to determine the rate constant k0 and the activation energy E from the given information in Table 1.   Form the logarithm of this equation to get at linear relation between y=log(k) and x=1/T, y=b+ax.  Use the given information about the reaction rate to generate several data points   Ti, ki.   (where T is the absolute temperature [K].)  Then set up a linear system of equations and solve them by the least squares method.    The idea is that we determine the coefficients a, b by finding a straight line y=b+ax that fits the given data points as well as possible, i.e., find a,b so that the norm of the residual b+axi-yi is minimized. Table 2 T [K] 343 353 363 373 383 393 403 k [s-1] 2.8e-5 5.6e-5 11.2e-5 22.4e-5 44.8e-5 89.6e-5 179.2e-5 Hint:   B=[ones(size(T)) -1./T].   The result should be somewhere around  k0=1e8, E=8e4.

Editor: Nils Svanstedt Last modified: Thu Jan 25 2007