(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2"; MacintoshStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times"; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times"; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times"; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times"; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times"; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-3, 12, "Courier"; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-3, 12, "Courier"; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R32768, L-3, 12, "Courier"; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-3, 12, "Courier"; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B32768, L-3, 12, "Courier"; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Geneva"; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = leftheader, inactive, L2, 12, "Times"; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; fontset = leftfooter, inactive, L2, 12, "Times"; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; paletteColors = 128; currentKernel; ] :[font = title; inactive; preserveAspect; startGroup] Linear Methods of Applied Mathematics Projections of vectors and functions ;[s] 3:0,1;38,2;39,1;76,-1; 3:0,25,18,Times,1,24,0,0,0;2,24,17,Times,1,23,0,0,0;1,12,8,Times,1,11,0,0,0; :[font = text; inactive; preserveAspect; plain; bold] (c) Copyright 1994-1997 by Evans M. Harrell II and James V. Herod. All rights reserved. :[font = text; inactive; Cclosed; preserveAspect; plain; bold; startGroup] Notes for the instructor. ;[s] 2:0,1;25,0;26,-1; 2:1,13,9,Times,1,12,0,0,0;1,11,8,Times,1,9,0,0,0; :[font = text; inactive; preserveAspect; endGroup] This contains calculations and examples which correlate with chapter 2 of the WWW text by Harrell and Herod. Students can be encouraged to cut and paste from this notebook to do homework. :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Instructions :[font = text; inactive; preserveAspect; endGroup] This notebook uses Mathematica to perform calculations for Harrell and Herod's hypertext book, Linear Methods of Applied Mathematics. The student needs only a basic knowledge of Mathematica to use the notebook, which is designed both to show how to work problems in the text and to provide a template for using Mathematica to work other problems of the student's own design. Calculations will be performed when the reader presses enter in a given calculation cell (bold print). It is best to activate the cells in order, so that Mathematica will be able to call on operators and functions defined in earlier cells. Red color coding is used to warn the reader when a given calculation relies on an earlier one. ;[s] 13:0,0;19,1;30,0;96,1;133,0;181,1;192,0;315,1;326,0;538,1;549,0;626,2;642,0;721,-1; 3:7,13,9,Times,0,12,0,0,0;5,13,9,Times,2,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = section; inactive; Cclosed; preserveAspect; startGroup] Projections in ordinary 3-D space :[font = text; inactive; preserveAspect; fontSize = 14] Let us begin by defining the projection of the vector v into the direction of w, for ordinary 3-D vectors v,w. :[font = input; preserveAspect] Proj[v_,w_] := w (v.w)/(w.w) :[font = text; inactive; preserveAspect; fontSize = 14] Notice that the projection has the direction of w but the length given by |v| cos(q), where q is the angle between v and w. You may wish to try out the command Proj for various pairs of vectors before proceeding. ;[s] 15:0,0;48,1;49,0;75,1;76,0;82,3;83,0;92,3;93,0;115,1;116,0;121,1;122,0;162,2;167,0;215,-1; 4:8,16,12,Times,0,14,0,0,0;4,16,12,Times,1,14,0,0,0;1,15,11,Courier,0,14,0,0,0;2,16,12,Symbol,0,14,0,0,0; :[font = text; inactive; preserveAspect; fontSize = 14] We do not always wish to project a vector onto another vector. For instance, we may wish to find the point in a plane which is closest to a given vector. Let us assume that a plane is specified by two orthogonal vectors (the orthogonality is important), such as: :[font = input; preserveAspect] P1 := {{1,1,1},{1,-1,0}} :[font = text; inactive; preserveAspect; fontSize = 14] The projection onto the plane is the sum of the projections onto the two orthogonal basis vectors: :[font = input; preserveAspect] Projp[v_,P_] := Proj[v,P[[1]]] + Proj[v,P[[2]]] ;[s] 5:0,0;16,1;20,0;33,1;37,0;48,-1; 2:3,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect; fontSize = 14] For example, as in the text: :[font = input; Cclosed; preserveAspect; startGroup] Projp[{3,4,5},P1] ;[s] 4:0,1;5,0;14,1;16,0;18,-1; 2:2,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {7/2, 9/2, 4} ;[o] 7 9 {-, -, 4} 2 2 :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] Model Problem II.8. We wish to find a) the vector in the plane spanned by {1,1,1} and {1,2,3}, which is the closest to v ={3, 4, 5}; and b) the vector in the same plane closest to {1,-1,0}. ;[s] 3:0,0;23,1;209,0;212,-1; 2:2,16,12,Times,1,14,0,0,0;1,16,12,Times,0,14,0,0,0; :[font = text; inactive; preserveAspect; fontSize = 14] Since the vectors spanning the plane are not given as orthogonal, we replace them with another pair: :[font = input; preserveAspect] v1 := {1,1,1} :[font = input; Cclosed; preserveAspect; startGroup] Solve[({1,2,3} - alpha {1,1,1}).{1,1,1} == 0, alpha] :[font = output; output; inactive; preserveAspect; endGroup] {{alpha -> 2}} ;[o] {{alpha -> 2}} :[font = input; Cclosed; preserveAspect; startGroup] v2 = {1,2,3} - alpha {1,1,1} /. % ;[s] 2:0,0;32,1;34,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {{-1, 0, 1}} ;[o] {{-1, 0, 1}} :[font = input; preserveAspect] P2 := {{1,1,1},{-1,0,1}} :[font = input; Cclosed; preserveAspect; startGroup] Projp[{3,4,5},P2] ;[s] 4:0,1;5,0;14,1;16,0;18,-1; 2:2,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {3, 4, 5} ;[o] {3, 4, 5} :[font = text; inactive; preserveAspect; fontSize = 14] Explanation: This vector is already in the plane, so the projection to the plane has no effect. Not so for the second vector: :[font = input; Cclosed; preserveAspect; startGroup] Projp[{1,-1,0},P2] ;[s] 4:0,1;5,0;15,1;17,0;19,-1; 2:2,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {1/2, 0, -1/2} ;[o] 1 1 {-, 0, -(-)} 2 2 :[font = text; inactive; preserveAspect; fontSize = 14] Now that it is in the plane, however, another projection won't move it: :[font = input; Cclosed; preserveAspect; startGroup] Projp[{1/2,0,-1/2},P2] ;[s] 4:0,1;5,0;19,1;21,0;23,-1; 2:2,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup; endGroup; endGroup] {1/2, 0, -1/2} ;[o] 1 1 {-, 0, -(-)} 2 2 :[font = section; inactive; Cclosed; preserveAspect; startGroup] Projections in function space :[font = text; inactive; preserveAspect; fontSize = 14] The following calculations do similar things in function space. The point here is to replace our original function with some other, convenient function, which is as close as possible to it. What do we mean by "close"? The notion of distance is given by the norm of a function: ||f - g||, which you may recognize from your courses on laboratory statistics as the r.m.s. deviation between f and g. The best approximation of a function by a nicer function will mean the specific nicer function with the smallest r.m.s. deviation. Let us begin by exploring this notion of distance: :[font = subsection; inactive; preserveAspect; startGroup] The distance between functions :[font = text; inactive; preserveAspect; fontSize = 14] We have learned that the r.m.s. distance between two real-valued functions defined on an interval a ˛ x ˛ b is: :[font = input; preserveAspect] RMSDist[f_,g_,{x_,a_,b_}] \ := Sqrt[Integrate[(f - g /. x -> intvar)^2, {intvar,a,b}]] :[font = text; inactive; preserveAspect; fontSize = 14] Warning! This formula does not work for complex-valued functions. The RMS norm of a single function f could be obtained with the command RMSDist[f,0,{x,a,b}]. ;[s] 3:0,1;8,0;151,2;173,-1; 3:1,16,12,Times,0,14,0,0,0;1,16,12,Times,1,14,0,0,0;1,15,11,Courier,0,14,0,0,0; :[font = input; preserveAspect] RMSNorm[f_,{x_,a_,b_}] := RMSDist[f,0,{x,a,b}] ;[s] 3:0,0;26,1;33,0;47,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect; fontSize = 14] What is the constant which is most like the cosine? 1/2? 1/4? Let's first check some values? :[font = input; Cclosed; preserveAspect; startGroup] {N[RMSDist[1/2,Cos[Pi x/2], {x,-1,1}]], N[RMSDist[1/4,Cos[Pi x/2], {x,-1,1}]], N[RMSDist[-1,Cos[Pi x/2], {x,-1,1}]]} ;[s] 7:0,0;3,1;10,0;43,1;50,0;84,1;91,0;120,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {0.4761937161122952, 0.6988420620085904, 2.355096407680655} ;[o] {0.476194, 0.698842, 2.3551} :[font = input; preserveAspect; endGroup] Plot[{Cos[Pi x/2], 1,1/2,-1}, {x, -1,1}] :[font = text; inactive; preserveAspect; fontSize = 14] The closest constant to Cos[Pi x/2] would be, geometrically speaking, the projection of the sine onto the direction of the constant functions. Intuitively, this should be the average of the function. As we learned in the text, in order to calculate the best approximation to a function, we perform a projection in function space, using virtually the same formula we use to project a vector onto another vector or onto a plane: ;[s] 3:0,0;4,1;11,0;430,-1; 2:2,16,12,Times,0,14,0,0,0;1,16,12,Times,1,14,0,0,0; :[font = input; preserveAspect] Projf[f_,g_, {x_,a_,b_}] := g Integrate[(f g /. x -> intvar), {intvar,a,b}] \ / Integrate[(g^2 /. x -> intvar), {intvar,a,b}] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Explanation of the syntax :[font = text; inactive; preserveAspect; endGroup] Ľ The syntax is analogous to the Mathematica command Integrate. You specify two functions (with the independent variables explicitly written), the variable on which they depend, and the range of that variable. That is, a < x < b. Ľ In the definition of the command the integration variable was named intvar because that is an unlikely variable to arise elsewhere in the formula. With this definition, the projected function and the original functions will depend on the same independent variable, which can be called anything (other than intvar in rare cases or reserved Mathematica words). ;[s] 11:0,0;34,2;45,0;54,1;64,0;307,1;313,0;546,1;552,0;579,2;590,0;600,-1; 3:6,13,9,Times,0,12,0,0,0;3,13,10,Courier,0,12,0,0,0;2,13,9,Times,2,12,0,0,0; :[font = text; inactive; preserveAspect; fontSize = 14] Indeed, we find that the projection of a function f[x] onto the direction of the constant function 1 is the average. For our example, :[font = input; Cclosed; preserveAspect; startGroup] Projf[Cos[Pi x/2],1, {x,-1,1}] ;[s] 2:0,1;5,0;31,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] 2/Pi ;[o] 2 -- Pi :[font = input; Cclosed; preserveAspect; startGroup] N[RMSDist[2/Pi,Cos[Pi x/2], {x,-1,1}]] ;[s] 3:0,0;2,1;9,0;39,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] 0.4352361782541725 ;[o] 0.435236 :[font = input; preserveAspect] Plot[{2/Pi,Cos[Pi x/2]}, {x,-1,1}] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] Some more examples of projection of functions :[font = text; inactive; preserveAspect; fontSize = 13] The projection of x onto the direction of x2: ;[s] 4:0,1;43,2;44,1;45,0;46,-1; 3:1,14,10,Times,0,13,0,0,0;2,16,12,Times,0,14,0,0,0;1,24,16,Times,32,14,0,0,0; :[font = input; Cclosed; preserveAspect; startGroup] Projf[x, x^2, {x, 0, 1}] :[font = output; output; inactive; preserveAspect; endGroup] (5*x^2)/4 ;[o] 2 5 x ---- 4 :[font = input; preserveAspect] Plot[{x, 5 x^2/4}, {x, 0, 1}] :[font = text; inactive; preserveAspect; fontSize = 14] The projection in the other direction. This time let's call the independent variable y: :[font = input; Cclosed; preserveAspect; startGroup] Projf[y^2, y, {y, 0, 1}] :[font = output; output; inactive; preserveAspect; endGroup] (3*y)/4 ;[o] 3 y --- 4 :[font = input; preserveAspect] Plot[{y^2, 3 y/4}, {y, 0, 1}] :[font = text; inactive; preserveAspect; endGroup] Notice that the coefficients are different when we project in different directions between the two functions. A question to ponder: What would happen if we projected x into the cirection of x2, then projected x2 back in the direction ox x, then repeated the process many times? ;[s] 6:0,1;195,2;196,1;214,2;215,1;283,0;284,-1; 3:1,13,9,Times,0,12,0,0,0;3,16,12,Times,0,14,0,0,0;2,24,16,Times,32,14,0,0,0; :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] Model Problem II.9 :[font = text; inactive; preserveAspect] Model Problem II.9. Consider the set of functions on the interval -1 ˛ x ˛ 1. We wish to find the function in the span of 1, x, and x2, which is the closest to f(x) = cos(šx/2). In other words, replace the cosine by a simple quadratic. ;[s] 5:0,1;18,2;135,3;136,2;240,0;241,-1; 4:1,13,9,Times,0,12,0,0,0;1,16,12,Times,1,14,0,0,0;2,16,12,Times,0,14,0,0,0;1,24,16,Times,32,14,0,0,0; :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Solution :[font = text; inactive; preserveAspect] Solution. First let us ask whether the three functions are orthogonal. Actually, no. The function x is orthogonal to each of the other two, but 1 and x2 are not orthogonal. We should replace x2 by x2 minus its projection onto the direction of 1, that is, x2 minus its average, which is easily found to be 1/3. Then we can project cos(šx/2) onto the span of the three functions 1, x, and x2 - 1/3: ;[s] 11:0,0;154,1;155,0;196,1;197,0;202,1;203,0;260,1;261,0;393,1;394,0;402,-1; 2:6,13,9,Times,0,12,0,0,0;5,21,13,Times,32,12,0,0,0; :[font = input; Cclosed; preserveAspect; startGroup] Projf[Cos[Pi x/2], 1, {x, -1,1}] ;[s] 2:0,1;5,0;33,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] 2/Pi ;[o] 2 -- Pi :[font = input; Cclosed; preserveAspect; startGroup] Projf[Cos[Pi x/2], x, {x, -1,1}] ;[s] 2:0,1;5,0;33,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] 0 ;[o] 0 :[font = text; inactive; preserveAspect] (This is also obvious without a calculation, since the cosine is even and x is odd.) :[font = input; Cclosed; preserveAspect; startGroup] Projf[Cos[Pi x/2], x^2 - 1/3, {x, -1,1}] ;[s] 2:0,1;5,0;41,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (45*((-2*(24 - 2*Pi^2))/(3*Pi^3) + (2*(-24 + 2*Pi^2))/(3*Pi^3))*(-1/3 + x^2))/8 ;[o] 2 2 -2 (24 - 2 Pi ) 2 (-24 + 2 Pi ) 1 2 45 (--------------- + ---------------) (-(-) + x ) 3 3 3 3 Pi 3 Pi -------------------------------------------------- 8 :[font = text; inactive; preserveAspect] The best quadratic approximation to a function on this interval would be: :[font = input; preserveAspect] BestApprox[f_] := Projf[f,1, {x, -1,1}] + \ Projf[f,x, {x, -1,1}] + \ Projf[f,x^2 - 1/3, {x, -1,1}] ;[s] 7:0,0;18,1;23,0;62,1;67,0;106,1;111,0;136,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect] We have the sum of three projections because there are three orthogonal functions in the set we are using to approximate the function f. Otherwise, the formula is just like the one for the plane. :[font = input; Cclosed; preserveAspect; startGroup] BestApprox[Cos[Pi x/2]] ;[s] 2:0,1;10,0;24,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] 2/Pi + (45*((-2*(24 - 2*Pi^2))/(3*Pi^3) + (2*(-24 + 2*Pi^2))/(3*Pi^3))*(-1/3 + x^2))/8 ;[o] 2 2 -2 (24 - 2 Pi ) 2 (-24 + 2 Pi ) 1 2 45 (--------------- + ---------------) (-(-) + x ) 3 3 3 2 3 Pi 3 Pi -- + -------------------------------------------------- Pi 8 :[font = input; Cclosed; preserveAspect; startGroup] Simplify[%] ;[s] 3:0,0;9,1;10,0;12,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (3*(20 - Pi^2 - 60*x^2 + 5*Pi^2*x^2))/Pi^3 ;[o] 2 2 2 2 3 (20 - Pi - 60 x + 5 Pi x ) ------------------------------- 3 Pi :[font = input; preserveAspect] Plot[{%,Cos[Pi x/2]},{x,-1,1}] ;[s] 3:0,0;6,1;7,0;31,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect] Let's compare with the Taylor approximation, :[font = input; Cclosed; preserveAspect; startGroup] Series[Cos[Pi x/2], {x, 0, 2}] :[font = output; output; inactive; preserveAspect] SeriesData[x, 0, {1, 0, -Pi^2/8}, 0, 3, 1] ;[o] 2 2 Pi x 3 1 - ------ + O[x] 8 :[font = input; preserveAspect; endGroup] SeriesData[x, 0, {1, 0, -Pi^2/8}, 0, 3, 1] ;[o] 2 2 Pi x 3 1 - ------ + O[x] 8 :[font = input; preserveAspect; endGroup; endGroup] Plot[{2/Pi + (45*((-2*(24 - 2*Pi^2))/(3*Pi^3) + (2*(-24 + 2*Pi^2))/(3*Pi^3))*(-1/3 + x^2))/8, \ Cos[Pi x/2], \ 1 - Pi^2 x^2 /8}, {x, -1,1}] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] A Variant of Model Problem II.9 :[font = text; inactive; preserveAspect] Model Problem II.9'. Consider the set of functions on the interval 0 ˛ x ˛ š. We wish to find the function in the span of 1, sin(2x) and sin(3x), which is the closest to f(x) = x2. In other words, find the function of the form a + b sin(2x) + c sin(3x) which is closest in the r.m.s. sense to x2. ;[s] 5:0,0;181,1;182,0;301,1;302,0;304,-1; 2:3,13,9,Times,0,12,0,0,0;2,21,13,Times,32,12,0,0,0; :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Solution :[font = text; inactive; preserveAspect] Solution. First let us ask whether the three functions are orthogonal. Actually, no. The function sin(2 x) is orthogonal to each of the other two, but 1 and sin(3x) are not orthogonal. We should replace sin(3x) by sin(3 x) minus its projection onto the direction of 1, that is, sin(3 x) minus its average. :[font = input; Cclosed; preserveAspect; startGroup] f3[x_] = Sin[3 x] - Projf[Sin[3 x],1,{x,0,Pi}] ;[s] 3:0,0;20,1;25,0;47,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] -2/(3*Pi) + Sin[3*x] ;[o] -2 ---- + Sin[3 x] 3 Pi :[font = input; preserveAspect] BestApprox[f_] := Projf[f,1,{x,0,Pi}] + \ Projf[f,Sin[2 x],{x,0,Pi}] + \ Projf[f,f3[x],0,Pi] ;[s] 7:0,0;18,1;23,0;60,1;65,0;109,1;114,0;129,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect] We have the sum of three projections because there are three orthogonal functions in the set we are using to approximate the function f. Otherwise, the formula is just like the one for the plane. :[font = input; Cclosed; preserveAspect; startGroup] BestApprox[x^2] ;[s] 2:0,1;10,0;16,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] Pi^2/3 + Projf[x^2, -2/(3*Pi) + Sin[3*x], 0, Pi] + (2*(-1/4 + (1 - 2*Pi^2)/4)*Sin[2*x])/Pi ;[o] 2 Pi 2 -2 --- + Projf[x , ---- + Sin[3 x], 0, Pi] + 3 3 Pi 2 1 1 - 2 Pi 2 (-(-) + ---------) Sin[2 x] 4 4 ----------------------------- Pi :[font = input; Cclosed; preserveAspect; startGroup] Simplify[%] ;[s] 3:0,0;9,1;10,0;12,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] Pi^2/3 + Projf[x^2, -2/(3*Pi) + Sin[3*x], 0, Pi] - Pi*Sin[2*x] ;[o] 2 Pi 2 -2 --- + Projf[x , ---- + Sin[3 x], 0, Pi] - Pi Sin[2 x] 3 3 Pi :[font = input; preserveAspect] Plot[{%,x^2},{x,0,Pi}] ;[s] 3:0,0;6,1;7,0;23,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect; endGroup; endGroup; endGroup] It is not a really close match (wait until the next section for something more impressive). On the other hand, there is definitely a correlation between the two graphs. :[font = section; inactive; Cclosed; preserveAspect; startGroup] The best polynomial (Legendre series) :[font = text; inactive; preserveAspect; fontSize = 14] We ought to do a better job of approximating Cos[Pi x/2] with a polynomial of 3rd or 4th degree. You may initially think of the Taylor series at x = 0, and wish regard Sin[x] as close to 0 + x + 0 - x^3/6. Or you might try the Taylor series at x= Pi, which is built into Mathematica in the Series command. ;[s] 3:0,0;273,1;284,0;308,-1; 2:2,16,12,Times,0,14,0,0,0;1,16,12,Times,2,14,0,0,0; :[font = input; Cclosed; preserveAspect; fontSize = 14; startGroup] Series[Cos[Pi x/2], {x, 0, 3}] :[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup] SeriesData[x, 0, {1, 0, -Pi^2/8}, 0, 4, 1] ;[o] 2 2 Pi x 4 1 - ------ + O[x] 8 :[font = input; preserveAspect; fontSize = 14] Plot[{1-(Pi x)^2 /8, Cos[Pi x/2]}, {x,-1,1}] :[font = input; Cclosed; preserveAspect; fontSize = 14; startGroup] N[RMSDist[1-(Pi x)^2 /8, Cos[Pi x/2], -1,1]] :[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup] 0.1118357132736736 ;[o] 0.111836 :[font = text; inactive; preserveAspect; fontSize = 14] Perhaps it is surprising to learn that there are other polynomial expressions which do a much better job of approximating a function such as the sine than does the Taylor series about any point. They are called the Legendre polynomials :[font = input; preserveAspect; fontSize = 14] P0[x_] := LegendreP[0,x] P1[x_] := LegendreP[1,x] P2[x_] := LegendreP[2,x] P3[x_] := LegendreP[3,x] :[font = input; preserveAspect; fontSize = 14] f[x_] := Cos[Pi x/2] :[font = input; Cclosed; preserveAspect; fontSize = 14; startGroup] Projf[f[x], P0[x], {x,-1,1}] ;[s] 5:0,0;1,1;6,0;13,1;15,0;30,-1; 2:3,14,11,Courier,1,14,0,0,0;2,14,11,Courier,1,14,65535,0,0; :[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup] 2/Pi ;[o] 2 -- Pi :[font = input; Cclosed; preserveAspect; fontSize = 14; startGroup] Projf[f[x], P1[x], {x,-1,1}] ;[s] 5:0,0;1,1;6,0;13,1;15,0;30,-1; 2:3,14,11,Courier,1,14,0,0,0;2,14,11,Courier,1,14,65535,0,0; :[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup] 0 ;[o] 0 :[font = input; Cclosed; preserveAspect; fontSize = 14; startGroup] Projf[f[x], P2[x], {x,-1,1}] ;[s] 5:0,0;1,1;6,0;13,1;15,0;30,-1; 2:3,14,11,Courier,1,14,0,0,0;2,14,11,Courier,1,14,65535,0,0; :[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup] (5*(-((24 - 2*Pi^2)/Pi^3) + (-24 + 2*Pi^2)/Pi^3)* (-1 + 3*x^2))/4 ;[o] 2 2 24 - 2 Pi -24 + 2 Pi 2 5 (-(----------) + -----------) (-1 + 3 x ) 3 3 Pi Pi ------------------------------------------- 4 :[font = input; Cclosed; preserveAspect; fontSize = 14; startGroup] % + %% + %%% ;[s] 5:0,1;1,0;4,1;6,0;9,1;13,-1; 2:2,14,11,Courier,1,14,0,0,0;3,14,11,Courier,1,14,65535,0,0; :[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup] 2/Pi + (5*(-((24 - 2*Pi^2)/Pi^3) + (-24 + 2*Pi^2)/Pi^3)* (-1 + 3*x^2))/4 ;[o] 2 2 24 - 2 Pi -24 + 2 Pi 2 5 (-(----------) + -----------) (-1 + 3 x ) 3 3 2 Pi Pi -- + ------------------------------------------- Pi 4 :[font = input; Cclosed; preserveAspect; fontSize = 14; startGroup] N[RMSDist[%, f[x],{x,-1,1}]] ;[s] 7:0,0;2,1;9,0;10,1;11,0;13,1;14,0;29,-1; 2:4,14,11,Courier,1,14,0,0,0;3,14,11,Courier,1,14,65535,0,0; :[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup] 0.02441441133097268 ;[o] 0.0244144 :[font = input; preserveAspect; fontSize = 14] Plot[{%%, f[x]}, {x,-1,1}] ;[s] 5:0,0;6,1;8,0;10,1;11,0;27,-1; 2:3,14,11,Courier,1,14,0,0,0;2,14,11,Courier,1,14,65535,0,0; :[font = text; inactive; preserveAspect; fontSize = 14; endGroup] How are these polynomials constructed? By convention, they are defined on -1 ˛ x ˛ 1, as has been done here. Essentially, we begin with the usual power functions and recombine them to find an orthogonal basis for function space. The procedure for doing this is due to Gram and Schmidt, and is based on a construction for ordinary vectors. :[font = section; inactive; Cclosed; preserveAspect; startGroup] Orthogonalization and the Gram-Schmidt method :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] Orthogonalization in ordinary 3-D or n-D space :[font = text; inactive; preserveAspect; fontSize = 14] If we have a list of orthonormal vectors v1, ...,vn, and we are given a new v, we can get Mathematica to orthogonalize the new vector by projecting away the components along v1, ...,vn. The inputs for the orthogonalization are a vector v, which is represented in Mathematica as a list {.....}, and a set of vectors, which is a list of lists. Here is one way to do it: ;[s] 11:0,0;42,1;43,0;50,1;51,0;90,2;101,0;176,1;177,0;184,1;185,0;371,-1; 3:6,16,12,Times,0,14,0,0,0;4,24,16,Times,64,14,0,0,0;1,16,12,Times,2,14,0,0,0; :[font = input; preserveAspect] Orth[v_,VecList_] := v - Sum[VecList[[k]] VecList[[k]].v, {k,1,Dimensions[VecList][[1]]}] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Explanation of the syntax :[font = text; inactive; preserveAspect; endGroup] Ľ If there are n vectors in the list, the command Dimensions[VecList][[1]] gives the number n, while Dimensions[VecList][[2]] would give the number of entries in each vector. A list of lists is a matrix, and the command Dimensions tells the size of the matrix, {number of rows, number of columns}. We need this because we want the command to work for any size list of vectors v1, ...,vn. ;[s] 11:0,0;80,1;105,0;132,1;170,0;266,1;279,0;437,2;438,0;445,2;446,0;448,-1; 3:6,13,9,Times,0,12,0,0,0;3,13,10,Courier,0,12,0,0,0;2,21,13,Times,64,12,0,0,0; :[font = text; inactive; preserveAspect; fontSize = 14] Warning: To use the command Orth properly, it is neccessary for the vectors in VecList to be orthonormal - orthogonal and normalized. ;[s] 2:0,1;7,0;135,-1; 2:1,16,12,Times,0,14,0,0,0;1,16,12,Times,1,14,0,0,0; :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Examples :[font = input; Cclosed; preserveAspect; startGroup] Orth[{1,1,1},{{1,0,0},{0,1,0}}] ;[s] 2:0,1;4,0;32,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {0, 0, 1} ;[o] {0, 0, 1} :[font = input; Cclosed; preserveAspect; startGroup] Orth[{1,1,0},{{1,0,0}}] ;[s] 2:0,1;4,0;24,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {0, 1, 0} ;[o] {0, 1, 0} :[font = input; Cclosed; preserveAspect; startGroup] Orth[{1,1,1},{{1,0,0},{0,1,0}}] ;[s] 2:0,1;4,0;32,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {0, 0, 1} ;[o] {0, 0, 1} :[font = input; preserveAspect] w1 := (1/Sqrt[3]){1,1,1} :[font = input; Cclosed; preserveAspect; startGroup] Orth[{1,1,0},{w1}] ;[s] 4:0,1;4,0;14,1;16,0;19,-1; 2:2,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {1/3, 1/3, -2/3} ;[o] 1 1 2 {-, -, -(-)} 3 3 3 :[font = text; inactive; preserveAspect] If we wish the result of this to be orthonormalized, we need to scale its length: :[font = input; Cclosed; preserveAspect; startGroup] w2 = %/Sqrt[%.%] ;[s] 8:0,1;2,0;5,1;6,0;12,1;13,0;14,1;15,0;17,-1; 2:4,12,10,Courier,1,12,0,0,0;4,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {6^(-1/2), 6^(-1/2), -(2/3)^(1/2)} ;[o] 1 1 2 {-------, -------, -Sqrt[-]} Sqrt[6] Sqrt[6] 3 :[font = input; Cclosed; preserveAspect; startGroup] Orth[{1,0,0},{w1,w2}] ;[s] 6:0,1;4,0;14,1;16,0;17,1;19,0;22,-1; 2:3,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] {1/2, -1/2, 0} ;[o] 1 1 {-, -(-), 0} 2 2 :[font = input; Cclosed; preserveAspect; startGroup] w3 = %/Sqrt[%.%] ;[s] 7:0,0;5,1;6,0;12,1;13,0;14,1;15,0;17,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup; endGroup; endGroup] {2^(-1/2), -2^(-1/2), 0} ;[o] 1 1 {-------, -(-------), 0} Sqrt[2] Sqrt[2] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] Orthogonalization in function space and construction of the Legendre polynomials :[font = text; inactive; preserveAspect; fontSize = 14] We illustrate the procedure to find a set of orthogonal polynomials, which we might use to approximate a complicated function representing realistic data as an idealized polynomial with the smallest loss of accuracy. The set we get is the Legendre polynomials. :[font = text; inactive; preserveAspect; fontSize = 14] The Legendre polynomials are thus just what you get by applying the Gram-Schmidt procedure to the power functions 1,x,x2, x3, etc., except that Legendre chose a convention which makes them orthogonal but not normalized. The normalized Legendre polynomials can be constructed as follows: ;[s] 5:0,0;119,1;120,0;123,1;124,0;289,-1; 2:3,16,12,Times,0,14,0,0,0;2,24,16,Times,32,14,0,0,0; :[font = input; preserveAspect] v0[x_] := 1 v1[x_] := x v2[x_] := x^2 v3[x_] := x^3 v4[x_] := x^4 :[font = input; Cclosed; preserveAspect; startGroup] p0[x_] = v0[x]/RMSNorm[v0[x],{x-1,1}] ;[s] 5:0,0;9,1;11,0;15,1;22,0;38,-1; 2:3,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] 2^(-1/2) ;[o] 1 ------- Sqrt[2] :[font = input; Cclosed; preserveAspect; startGroup] v1[x] - Projf[v1[x],p0[x],{x,-1,1}] ;[s] 8:0,1;2,0;8,1;13,0;14,1;16,0;20,1;22,0;36,-1; 2:4,12,10,Courier,1,12,0,0,0;4,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] x ;[o] x :[font = input; Cclosed; preserveAspect; startGroup] p1[x_] = %/RMSNorm[%,{x-1,1}] ;[s] 7:0,0;9,1;10,0;11,1;18,0;19,1;20,0;30,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (3/2)^(1/2)*x ;[o] 3 Sqrt[-] x 2 :[font = input; Cclosed; preserveAspect; startGroup] v2[x] - Projf[v2,p0,{x-1,1}] \ - Projf[v2,p1,{x-1,1}] ;[s] 14:0,1;2,0;9,1;14,0;15,1;17,0;18,1;20,0;35,1;40,0;41,1;43,0;44,1;46,0;56,-1; 2:7,12,10,Courier,1,12,0,0,0;7,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] -1/3 + x^2 ;[o] 1 2 -(-) + x 3 :[font = input; Cclosed; preserveAspect; startGroup] p2[x_] = %/RMSNorm[%,{x-1,1}] ;[s] 7:0,0;9,1;10,0;11,1;18,0;19,1;20,0;30,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (3*(5/2)^(1/2)*(-1/3 + x^2))/2 ;[o] 5 1 2 3 Sqrt[-] (-(-) + x ) 2 3 --------------------- 2 :[font = input; Cclosed; preserveAspect; startGroup] v3[x] - Projf[v3,p0,{x-1,1}] \ - Projf[v3,p1,{x-1,1}] \ - Projf[v3,p2,{x-1,1}] ;[s] 20:0,1;2,0;9,1;14,0;15,1;17,0;18,1;20,0;39,1;44,0;45,1;47,0;48,1;50,0;68,1;73,0;74,1;76,0;77,1;79,0;89,-1; 2:10,12,10,Courier,1,12,0,0,0;10,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (-3*x)/5 + x^3 ;[o] -3 x 3 ---- + x 5 :[font = input; Cclosed; preserveAspect; startGroup] p3[x_] = %/RMSNorm[%,{x-1,1}] ;[s] 7:0,0;9,1;10,0;11,1;18,0;19,1;20,0;30,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (5*(7/2)^(1/2)*((-3*x)/5 + x^3))/2 ;[o] 7 -3 x 3 5 Sqrt[-] (---- + x ) 2 5 --------------------- 2 :[font = input; Cclosed; preserveAspect; startGroup] v4[x] - Projf[v4,p0,{x-1,1}] \ - Projf[v4,p1,{x-1,1}] \ - Projf[v4,p2,{x-1,1}] \ - Projf[v4,p3,{x-1,1}] ;[s] 26:0,1;2,0;8,1;13,0;14,1;16,0;17,1;19,0;38,1;43,0;44,1;46,0;47,1;49,0;67,1;72,0;73,1;75,0;76,1;78,0;96,1;101,0;102,1;104,0;105,1;107,0;117,-1; 2:13,12,10,Courier,1,12,0,0,0;13,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] -1/5 + x^4 - (3*(5/2)^(1/2)*((3*(5/2)^(1/2))/7 - 10^(-1/2))*(-1/3 + x^2))/2 ;[o] 5 3 Sqrt[-] 5 2 1 1 2 3 Sqrt[-] (--------- - --------) (-(-) + x ) 1 4 2 7 Sqrt[10] 3 -(-) + x - -------------------------------------------- 5 2 :[font = input; Cclosed; preserveAspect; startGroup] Simplify[%] ;[s] 3:0,0;9,1;10,0;12,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] 3/35 - (6*x^2)/7 + x^4 ;[o] 2 3 6 x 4 -- - ---- + x 35 7 :[font = input; Cclosed; preserveAspect; startGroup] p4[x_] = %/RMSNorm[%,{x-1,1}] ;[s] 7:0,0;9,1;10,0;11,1;18,0;19,1;20,0;30,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (105*(3/35 - (6*x^2)/7 + x^4))/(8*2^(1/2)) ;[o] 2 3 6 x 4 105 (-- - ---- + x ) 35 7 -------------------- 8 Sqrt[2] :[font = text; inactive; preserveAspect; fontSize = 14] It may be interesting to Plot the Legendre polynomials and try to visualize their orthogonality: :[font = input; preserveAspect] Plot[{LegendreP[0,x], LegendreP[1,x], LegendreP[2,x], \ LegendreP[3,x], LegendreP[4,x]}, {x,-1,1}] :[font = text; inactive; preserveAspect; fontSize = 14] What happens if the interval is not -1 ˛ x ˛ 1? Let's examine a similar example on another interval. Specifically, let's try to approximate the Sin[x] on the interval 0 ˛ x ˛ š. :[font = text; inactive; preserveAspect; fontSize = 14] One possibility would be to construct a set of orthogonal polynomials on this new interval by beginning with the functions 1, x, x2, x3, etc., and then applying the Gram-Schmidt procedure to get a set of orthogonal polynomials. Because the interval is different, these will definitely be different from the Legendre polynomials we found above - orthogonality depends on the interval as much as it depends on the function. In particular, the new polynomials won't be nicely even or odd like the ones given above, because the interval is not balanced around 0. Here is a shortcut, however. Let us change the variable x in a simple, linear way in order to shift the interval [0, š] to the one we already understand. Our new polynomials will be: ;[s] 5:0,0;130,1;131,0;134,1;135,0;747,-1; 2:3,16,12,Times,0,14,0,0,0;2,24,16,Times,32,14,0,0,0; :[font = input; Cclosed; preserveAspect; startGroup] P0[x_] = LegendreP[0,2 (x-Pi/2)/Pi] P1[x_] = LegendreP[1,2 (x-Pi/2)/Pi] P2[x_] = LegendreP[2,2 (x-Pi/2)/Pi] P3[x_] = LegendreP[3,2 (x-Pi/2)/Pi] P4[x_] = LegendreP[4,2 (x-Pi/2)/Pi] :[font = output; output; inactive; preserveAspect] 1 ;[o] 1 :[font = output; output; inactive; preserveAspect] (-Pi + 2*x)/Pi ;[o] -Pi + 2 x --------- Pi :[font = output; output; inactive; preserveAspect] (Pi^2 - 6*Pi*x + 6*x^2)/Pi^2 ;[o] 2 2 Pi - 6 Pi x + 6 x ------------------- 2 Pi :[font = output; output; inactive; preserveAspect] (-Pi^3 + 12*Pi^2*x - 30*Pi*x^2 + 20*x^3)/Pi^3 ;[o] 3 2 2 3 -Pi + 12 Pi x - 30 Pi x + 20 x ---------------------------------- 3 Pi :[font = output; output; inactive; preserveAspect; endGroup] (Pi^4 - 20*Pi^3*x + 90*Pi^2*x^2 - 140*Pi*x^3 + 70*x^4)/ Pi^4 ;[o] 4 3 2 2 3 4 Pi - 20 Pi x + 90 Pi x - 140 Pi x + 70 x ---------------------------------------------- 4 Pi :[font = text; inactive; preserveAspect; fontSize = 14] It is not hard to check directly that these shifted Legendre polynomials are orthogonal. :[font = text; inactive; preserveAspect; fontSize = 14; endGroup; endGroup; endGroup] Study question. What is the general rule to shift to an interval a ˛ x ˛ b? ;[s] 2:0,1;14,0;77,-1; 2:1,16,12,Times,0,14,0,0,0;1,16,12,Times,1,14,0,0,0; ^*)