ALA C, 2003-04, vecka 1, Studio 1b

[Kurssidan, vecka 1, vecka 2, vecka 3, vecka 4, vecka 5, vecka 6, vecka 7, Matlab] Linjär algebra i R^n: Minsta kvadrat-metoden, dynamiska system 1.  Least squares methods, Ch 42.45. Define the matrix B and vector u as in the previous exercise (Studio 1a). 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. 4.  Dynamical systems and linearization Do problems 91.4-5 d) from 91. Linearization and stability (ps). Use plots with the solution plotted against time as well as phase portraits (if you have forgotten about phase portraits, check studio 5.1 from ALA-B).

Editor: Stig Larsson Last modified: Sun Jan 18 21:30:05 MET 2004