# Homework #5 > # Name: Key # # Collaborators: # # NOTE: You may work on the problems together, but you must formulate # the phrasing of your answers independantly. In particular you must # type your own submission. > with(linalg): > grade := 0; grade := 0 # Problem #1 (Least Squares). # # You are given the following data sampling points in the form [x,y,z], # and accurate to about 2 decimal places from a surface defined by z = # F(x,y), where the function F is unknown. > > values := [[3.8000, -.60000, 141.38], [1.8000, 4.3000, 140.17], > [3.3000, -10., 116.50], [4., 6., 158.40], [-.60000, 9.8000, 135.42], > [-1.7000, 4.8000, 118.57], [2.7000, -3.3000, 128.71], [.40000, 1.1000, > 124.64], [-2.8000, -3.7000, 96.843], [1.7000, -2.9000, 123.92]]; values := [[3.8000, -.60000, 141.38], [1.8000, 4.3000, 140.17], [3.3000, -10., 116.50], [4., 6., 158.40], [-.60000, 9.8000, 135.42], [-1.7000, 4.8000, 118.57], [2.7000, -3.3000, 128.71], [.40000, 1.1000, 124.64], [-2.8000, -3.7000, 96.843], [1.7000, -2.9000, 123.92]] # a) If we guess that the surface is in fact plane, then the equation # it must satisfy is of the form z = a*x + b*y + c*1. Find a matrix A # and a vector u so that a solution to A*v = u would give the # coefficients of the equation (i.e. v = [a, b, c]) if all the # measurements were prefectly accurate. > A := matrix(10,3); u:= vector(10); A := array(1 .. 10, 1 .. 3, []) u := array(1 .. 10, []) > for i from 1 to 10 do A[i,1] := values[i][1]; A[i,2] := values[i][2]; > A[i,3] := 1; u[i] := values[i][3]; od: print(A); print(u); [3.8000 -.60000 1] [ ] [1.8000 4.3000 1] [ ] [3.3000 -10. 1] [ ] [ 4. 6. 1] [ ] [-.60000 9.8000 1] [ ] [-1.7000 4.8000 1] [ ] [2.7000 -3.3000 1] [ ] [.40000 1.1000 1] [ ] [-2.8000 -3.7000 1] [ ] [1.7000 -2.9000 1] [141.38, 140.17, 116.50, 158.40, 135.42, 118.57, 128.71, 124.64, 96.843, 123.92] > # 3/3 > grade := grade + 3; grade := 3 > # b) Verify that the system does not have a solution. Find the least # squares solution to the system. > > linsolve(A,u); # No solution. > vhat := linalg[linsolve](transpose(A)&*A,transpose(A)&*u); # least > squares solution vhat := [6.02814772, 2.142869568, 119.6812556] > # 3/3 > grade := grade +3; grade := 6 # c) How far off is your this solution? (i.e. what is the norm of the # difference A * v - u). > norm(u - A &* vhat,2); 3.736800931 > # Pretty far off. > # 3/3 > grade := grade + 3; grade := 9 # d) It may be that this set of points did not lie in a plane. Perhaps # they lie on a quadratic surface, i.e. their equation is one of the # form: # z = a*x^2 + b*y^2 + c*x*y + d* x + e*y + f. Repeat steps (a)-(c) # this time finding the least squares solution v = [a,b,c,d,e,f]. How # does the norm of A*v - u, compare to the norm in part (c)? > A := matrix(10,6); b:= vector(10); for i from 1 to 10 do A[i,4] := > values[i][1]; A[i,5] := values[i][2]; A[i,6] := 1; A[i,1] := A[i,4]^2; > A[i,2] := A[i,5]^2; A[i,3] := A[i,4] * A[i,5]; b[i] := values[i][3]; > od: print(A); print(b); A := array(1 .. 10, 1 .. 6, []) b := array(1 .. 10, []) [14.44000000 , .3600000000 , -2.280000000 , 3.8000 , -.60000 , 1] [3.24000000 , 18.49000000 , 7.74000000 , 1.8000 , 4.3000 , 1] [10.89000000 , 100. , -33.0000 , 3.3000 , -10. , 1] [16. , 36. , 24. , 4. , 6. , 1] [.3600000000 , 96.04000000 , -5.880000000 , -.60000 , 9.8000 , 1] [2.89000000 , 23.04000000 , -8.16000000 , -1.7000 , 4.8000 , 1] [7.29000000 , 10.89000000 , -8.91000000 , 2.7000 , -3.3000 , 1] [.1600000000 , 1.21000000 , .440000000 , .40000 , 1.1000 , 1] [7.84000000 , 13.69000000 , 10.36000000 , -2.8000 , -3.7000 , 1] [2.89000000 , 8.41000000 , -4.93000000 , 1.7000 , -2.9000 , 1 ] [141.38, 140.17, 116.50, 158.40, 135.42, 118.57, 128.71, 124.64, 96.843, 123.92] > vhat2 := linalg[linsolve](transpose(A)&*A,transpose(A)&*b); vhat2 := [-.000349227, .000074482, .100314044, 6.000742758, 1.999212234, 120.0016548] > norm(b-A &* vhat2,2); .02145214926 > # Quite a bit better than the other fit. > > # 6/6 > grade := grade + 6; grade := 15 > > # NOTE: the numbers for the problem were obtained by evaluating the > function F(x,y) = 6*x + 2*y + .1 * x*y + 120 + > rand(-100..100)()/10000. > # Problem #2 (Inner Products) # # Let x = (x1,x2,x3) and y = (y1,y2,y3). Determine whether the # following are inner products on R^3: > # a) = x1*y3 + x2*y2 + x3*y1 > # Symmetric? Yes. # = y1*x3 - x2*y2 + y3*x1 = # Bilinear? Yes. # = a * + # Positive Definite? No. # = 2*x1*x3 + x2^2 =0 if x1 = -1 , x2 = 1, x3 = 1/2 # b) = y1*x2+y2*x3+x3*y1 # Symmetric? No # - = ( x1*y2 + x2*y3 + y3*x1) - (y1*x2+y2*x3+x3*y1) != 0 # in general > > > # c) = x1^2 + y2^2 + x3*y3 > # Bilinear? No. > > > > # 6/6 > > grade := grade + 6; grade := 21 # Problem #3 (Orthogonalization) # > # a) Use the Gram-Schmidt process to find an orthogonal basis for the # span of the vectors: # u = (1, 3, 4, 6) # v = (2, 7, 7, 12) # w = (4, 16, 12, 24) # x = (3, 10, 11, 18) > u := vector([1, 3, 4, 6]); > v := vector([2, 7, 7, 12]); > w := vector([4, 16, 12, 24]); > x := vector([3, 10, 11, 18]); u := [1, 3, 4, 6] v := [2, 7, 7, 12] w := [4, 16, 12, 24] x := [3, 10, 11, 18] > base := GramSchmidt([u,v,w,x]); 65 -29 base := [[1, 3, 4, 6], [1/62, --, ---, 3/31]] 62 31 > > # 3/3 > grade := grade +3; grade := 24 # b) Use the orthogonal basis to compute the projection of b = (1, 1, # 1, 1) onto span(u,v,w,x). > > b := vector([1, 1, 1, 1]); b := [1, 1, 1, 1] > b1:=evalm(dotprod(b,base[1]) * base[1]/norm(base[1],2)^2 + > dotprod(b,base[2]) * base[2]/norm(base[2],2)^2); [28 98 98 56] b1 := [---, ---, ---, --] [123 123 123 41] > # 3/3 > grade := grade +3; grade := 27 # c) Compare your answer in (b) to what you get from computing the least # squares solution to A y = b, where A the columns of A are the taken # from a maximal linearly independant subset of {u,v,w,x}, and b is as # in part (b) (i.e find the b_hat in the column space of A which is # closest to b). > # Note: u and v are linearly independant > A := transpose(matrix([[1, 3, 4, 6],[2, 7, 7, 12]])); [1 2] [ ] [3 7] A := [ ] [4 7] [ ] [6 12] > xhat := linsolve(transpose(A)&*A,transpose(A)&*b); [ 14 ] xhat := [0, ---] [ 123] > b2 := evalm(A&*xhat); [28 98 98 56] b2 := [---, ---, ---, --] [123 123 123 41] > # They are the same! > # 3/3 > grade := grade + 3; grade := 30 >