(************** Content-type: application/mathematica ************** Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 33792, 959]*) (*NotebookOutlinePosition[ 34562, 985]*) (* CellTagsIndexPosition[ 34518, 981]*) (*WindowFrame->Normal*) Notebook[{ Cell["Lab 10: Integrate This!", "Title", Background->RGBColor[0, 0, 1]], Cell["Math 342\tSpring 2003", "Subtitle"], Cell[BoxData[{ \(\(<< Graphics`FilledPlot`;\)\), "\n", \(\(Off[General::spell];\)\), "\[IndentingNewLine]", \(\(Off[General::spell1];\)\), "\n", \(\(Off[NIntegrate::ploss];\)\), "\n", \(\(Off[NIntegrate::ncvb];\)\), "\n", \(\(\(Off[$MaxExtraPrecision::meprec];\)\(\n\) \)\), "\n", \(\(Clear[riemannsum];\)\), "\n", \(\(riemannsum[fnc_, tind_, to_, tf_, nint_, rtorlt_, fncicin_: Automatic] := Block[{}, If[fncicin === Automatic, fncic = 0, fncic = fncicin]; \[IndentingNewLine]\n f[t_] = fnc /. tind -> t; \n\t\t\n F[t_] = \(NDSolve[{\(y'\)[t] == f[t], y[to] == fncic}, y[t], {t, to - 0.1*\((tf - to)\), tf + 0.1*\((tf - to)\)}]\)[\([1, 1, 2]\)]; \n\n h = \((tf - to)\)/nint; \n If[rtorlt == "\", discvalues = Table[f[t] // N, {t, to, tf - h, h}], discvalues = Table[f[t] // N, {t, to + h, tf, h}]]; \n\n Fdiscvalues = Table[{to + k*h, fncic + Sum[discvalues[\([i]\)]*h, {i, 1, k}]}, {k, 1, nint}]; \n\t\tFdiscvalues = Prepend[Fdiscvalues, {to, fncic}]; \n\t\n fvalues = Table[f[t], {t, to, tf, \((tf - to)\)/500}]; \n Fvalues = Table[F[t], {t, to, tf, \((tf - to)\)/500}]; \n Fmax = Max[ Flatten[{Fvalues, Table[Fdiscvalues[\([i, 2]\)], {i, 1, nint}]}]]; \n Fmin = Min[ Flatten[{Fvalues, Table[Fdiscvalues[\([i, 2]\)], {i, 1, nint}]}]]; \n\n fmaxt = Max[fvalues]; \nfmint = Min[fvalues]; \nfmin = fmint; \n fmax = fmaxt; \n\nchangesign = {}; \n Do[If[Sign[discvalues[\([k]\)]] != Sign[discvalues[\([k + 1]\)]], changesign = Append[changesign, to + k*h]], {k, 1, nint - 1}]; \nchangesign = Prepend[changesign, to]; \n changesign = Append[changesign, tf]; \n\t\nsize = 7; \t\t\n steprat = 0.020; \n\t\thrange = If[Fmax <= 0, {1.1*Fmin, 0}, If[Fmin >= 0, {0, 1.1*Fmax}, {1.1*Fmin, 1.1*Fmax}]]; \n\t\tfrange = If[fmax <= 0, {1.1*fmin, 0}, If[fmin >= 0, {0, 1.1*fmax}, {1.1*fmin, 1.1*fmax}]]; \n\t\t\n p1 = Plot[f[t], {t, to, tf}, PlotRange -> {{to, tf}, frange}, AspectRatio -> 1, PlotStyle -> {RGBColor[0, 0, 0], Thickness[0.008]}, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, DisplayFunction -> Identity]; \n\t\t\n\t p2t = {}; \n\t\t\n\t\t\t\tDo[ If[discvalues[\([k]\)] > 0, fillcolor = RGBColor[1, 0, 0], fillcolor = RGBColor[0, 0, 1]]; \n\t\t\t\t\tpinter = FilledPlot[ discvalues[\([k]\)], {t, to + \((k - 1)\)*h, to + k*h}, PlotRange -> {{to, tf}, frange}, PlotStyle -> {RGBColor[0, 0, 0], Thickness[0.008]}, Fills -> {fillcolor}, AxesFront -> True, Curves -> None, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, DisplayFunction -> Identity]; \n\t\t\t\t\tp2t = Append[p2t, pinter], {k, 1, nint}]; \n\t\t\n\t\tp2t = Append[p2t, p1]; \n\t\t\n p2 = Show[p2t, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, AxesLabel \[Rule] {x, fnc}, \ DisplayFunction -> Identity]; \n\t\t\n p3 = Plot[F[t], {t, to, tf}, PlotRange -> hrange, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, DisplayFunction -> Identity]; \n p4 = ListPlot[Fdiscvalues, PlotRange -> hrange, PlotStyle -> {PointSize[0.02], RGBColor[0.398444, \ 0.847669, \ 0.347662]}, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, DisplayFunction -> Identity]; \n p5 = Show[{p4, p3}, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, \ DisplayFunction -> Identity]; \n\np6t = {}; \n\t\t\n Do[If[discvalues[\([k]\)] - f[to + \((k - 0.5)\)*h] > 0, fillcolor = RGBColor[1, 0, 0], fillcolor = RGBColor[0, 0, 1]]; \n\t\t\t\t\tpinter = FilledPlot[{discvalues[\([k]\)], f[t]}, {t, to + \((k - 1)\)*h, to + k*h}, PlotRange -> {{to, tf}, frange}, PlotStyle -> {RGBColor[0, 0, 0], Thickness[0.008]}, Fills -> {fillcolor}, AxesFront -> True, Curves -> None, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, DisplayFunction -> Identity]; \n\t\t\t\t\tp6t = Append[p6t, pinter], {k, 1, nint}]; \n\t\t\n (*\(p6t = Append[p6t, p1];\)*) \n\t\t\n p6 = Show[p6t, ImageSize -> {72*size, 72*size/3}, PlotRegion -> {{0, 1}, {0, 1}}, AxesLabel \[Rule] {x, "\"}, \ DisplayFunction -> Identity]; \n\n accumerror = Table[{to + \((k - 1)\)*h, Fdiscvalues[\([k, 2]\)] - F[to + \((k - 1)\)*h]}, {k, 1, nint}]; \t\t\n p7 = ListPlot[accumerror, PlotStyle -> {PointSize[0.020]}, AxesLabel \[Rule] {x, f}, DisplayFunction -> Identity]; \t\n\t\t\t\t\t\t\n Show[GraphicsArray[{p2, p6}], ImageSize -> {72*size, 72*size/2}, DisplayFunction -> $DisplayFunction];\n\n\t\t\n\t\t\t\t\t\t\t\t\t\ \t];\)\)}], "Input", PageWidth->PaperWidth, CellOpen->False, InitializationCell->True], Cell[CellGroupData[{ Cell["Introduction", "Section"], Cell[TextData[{ "Today, we'll explore the Newton-Cotes numerical integration formulas, \ their derivation, and error analysis. Then, we'll look at how ", StyleBox["Mathematica", FontSlant->"Italic"], " integrates..." }], "Text", TextAlignment->Left, TextJustification->0] }, Closed]], Cell[CellGroupData[{ Cell[TextData[{ "A Picture Paints a Thousand Words \n", StyleBox["Bonus: find the hidden homework problems!", FontSize->12] }], "Section"], Cell[TextData[{ "Back in the early days of calculus... you first learned how to numerically \ integrate using Riemann Sums. We can define a command that shows the areas \ of the rectangles that are accumulated in the Riemann sum and the errors in \ the estimate. The command is:\n\t\t ", StyleBox["riemannsum[f_, x_, a_, b_, n_, rightleft_]", FontWeight->"Bold"], ". \nThe arguments are the function, ", StyleBox["f", FontWeight->"Bold"], ", the independent variable, ", StyleBox["x", FontWeight->"Bold"], ", the lower bound on the integral, ", StyleBox["a", FontWeight->"Bold"], ", the upper bound, ", StyleBox["b", FontWeight->"Bold"], ", the number of rectangles, ", StyleBox["n", FontWeight->"Bold"], ", and a \"right\" or \"left\" sum indicator. Red areas are positive, and \ blue areas are negative. The error is taken as the Riemann estimate minus the \ exact value of the integral.\nHere is a visualization of what happens when we \ use Riemann Sums to integrate. " }], "Text", TextAlignment->Left, TextJustification->0], Cell["\<\ Thanks to Dr. John Scharf for the code for this and several other commands we \ are using today!\ \>", "Commentary", FontFamily->"Arial", FontColor->RGBColor[0, 0, 1]], Cell[BoxData[{ \(\(Clear[n];\)\), "\[IndentingNewLine]", \(\(n = 5;\)\), "\n", \(\(riemannsum[Sin[x], x, 0, Pi/2, n, "\"];\)\), "\n", \(\(riemannsum[Sin[x], x, 0, Pi/2, n, "\"];\)\)}], "Input", PageWidth->PaperWidth], Cell[TextData[{ "This gives us a good visualization of what is happening... \n", StyleBox["Factoid: The \"average\" of the left and right Riemann sums \ gives us the Trapezoid Rule!\n", FontColor->RGBColor[1, 0, 0]], "Which gives me an idea... let's define:" }], "Text"], Cell[TextData[{ " ", StyleBox["riemSumLeft[f_,x_,a_,b_,n_] ", FontWeight->"Bold"], "and", StyleBox[" riemSumRight[f_,x_,a_,b_,n_] ", FontWeight->"Bold"], "commands to calculate left-hand and right-hand Riemann sums. The arguments \ are ", StyleBox["f", FontWeight->"Bold"], ", the integrand function, ", StyleBox["x", FontWeight->"Bold"], ", the independent variable, ", StyleBox["a", FontWeight->"Bold"], ", the lower bound, ", StyleBox["b", FontWeight->"Bold"], ", the upper bound, and ", StyleBox["n", FontWeight->"Bold"], ", the number of subintervals." }], "Text", PageWidth->PaperWidth], Cell[BoxData[{ \(\(\(riemSumLeft[f_, x_, a_, b_, n_] := Sum[\((f /. x \[Rule] \((a + \((i - 1)\)*\((b - a)\)/n)\))\)*\((b - a)\)/ n, {i, 1, n}];\)\(\n\) \)\), "\n", \(\(riemSumRight[f_, x_, a_, b_, n_] := Sum[\((f /. x \[Rule] \((a + i*\((b - a)\)/n)\))\)*\((b - a)\)/n, {i, 1, n}];\)\)}], "Input", PageWidth->PaperWidth], Cell["Here is how we use it:", "Text"], Cell[BoxData[{ \(\[IndentingNewLine]\(f[x_] = Sin[x];\)\), "\[IndentingNewLine]", \(\(a = 0;\)\), "\[IndentingNewLine]", \(\(b = \[Pi]\/2;\)\), "\[IndentingNewLine]", \(\(exactvalue = \[Integral]\_a\%b f[ x] \[DifferentialD]x;\)\), "\[IndentingNewLine]", \(\(\(Print["\", exactvalue];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(t1 = Table[{2^i, \(\[Pi]/2\)/2^i, RleftI = riemSumLeft[f[x], x, a, b, 2^i] // N, RrightI = riemSumRight[f[x], x, a, b, 2^i] // N, RleftI - exactvalue, RrightI - exactvalue}, {i, 1, 11}];\)\(\n\) \)\), "\n", \(TableForm[t1, TableHeadings \[Rule] {None, {"\", "\", "\", \ "\", "\", "\"}}]\)}], "Input", PageWidth->PaperWidth], Cell["Now, we can define the Trapezoid rule:", "Text"], Cell[BoxData[ \(\(trapezoid[f_, x_, a_, b_, n_] := \((riemSumLeft[f, x, a, b, n] + riemSumRight[f, x, a, b, n])\)/2;\)\)], "Input", PageWidth->PaperWidth], Cell["Now we calculate the estimates of the integral.", "Text", PageWidth->PaperWidth], Cell[BoxData[{ \(\(n = 10;\)\), "\[IndentingNewLine]", \(trapezoid[f[x], x, a, b, n] // N\), "\n", \(riemSumLeft[f[x], x, a, b, n] // N\), "\n", \(riemSumRight[f[x], x, a, b, n] // N\)}], "Input", PageWidth->PaperWidth], Cell["And we compute the errors.", "Text", PageWidth->PaperWidth], Cell[BoxData[{ \(\((trapezoid[f[x], x, a, b, n] // N)\) - exactvalue\), "\n", \(\((riemSumLeft[f[x], x, a, b, n] // N)\) - exactvalue\), "\n", \(\((riemSumRight[f[x], x, a, b, n] // N)\) - exactvalue\)}], "Input", PageWidth->PaperWidth], Cell[TextData[{ "An intermediate improvement that we can make is to define the ", StyleBox["Riemann Midpoint Formula", FontWeight->"Bold"], ". The idea here is that we use the midpoint of each subinterval in a \ Riemann sum. " }], "Text"], Cell[TextData[StyleBox["Question: Why would you want to do this? Think in \ terms of how the error will be affected by using the midpoint... say if a \ function is increasing or decreasing over the interval of integration...", FontWeight->"Bold"]], "Text", Background->RGBColor[0, 1, 1]], Cell[BoxData[ \(\(riemSumMid[f_, x_, a_, b_, n_] := Sum[\((f /. x \[Rule] \((a + \((i - 1/2)\)*\((b - a)\)/n)\))\)*\((b - a)\)/n, {i, 1, n}];\)\)], "Input", PageWidth->PaperWidth], Cell["Check it out!", "Text"], Cell[BoxData[{ \(\(n = 10;\)\), "\[IndentingNewLine]", \(riemSumMid[f[x], x, a, b, n] // N\), "\[IndentingNewLine]", \(trapezoid[f[x], x, a, b, n] // N\), "\n", \(riemSumLeft[f[x], x, a, b, n] // N\), "\n", \(riemSumRight[f[x], x, a, b, n] // N\)}], "Input", PageWidth->PaperWidth], Cell["Again, we compute the errors.", "Text", PageWidth->PaperWidth], Cell[BoxData[{ \(\((riemSumMid[f[x], x, a, b, n] // N)\) - exactvalue\), "\[IndentingNewLine]", \(\((trapezoid[f[x], x, a, b, n] // N)\) - exactvalue\), "\n", \(\((riemSumLeft[f[x], x, a, b, n] // N)\) - exactvalue\), "\n", \(\((riemSumRight[f[x], x, a, b, n] // N)\) - exactvalue\)}], "Input", PageWidth->PaperWidth], Cell[TextData[StyleBox["Question:\tWhich technique gives you the best \ solution for a given number of intervals?", FontWeight->"Bold"]], "Text", PageWidth->PaperWidth, CellMargins->{{Inherited, -199}, {Inherited, Inherited}}, Background->RGBColor[0, 1, 1]], Cell[TextData[{ "Now, let's define Simpson's 1/3 rule, with a new perspective: Since the \ magnitude of the error for the trapezoid method is approximately two times \ the magnitude of the error for a midpoint Riemann sum with the same number of \ subintervals, we can add ", StyleBox["two", FontWeight->"Bold"], " midpoint estimates and ", StyleBox["one", FontWeight->"Bold"], " trapezoid estimate and divide the result by ", StyleBox["three", FontWeight->"Bold"], ". The errors for the two methods should nearly cancel. What we end up \ with is Simpson's 1/3 method." }], "Text"], Cell[BoxData[ \(\(simpson[f_, x_, a_, b_, n_] := \((2*riemSumMid[f, x, a, b, n/2] + trapezoid[f, x, a, b, n/2])\)/3;\)\)], "Input", PageWidth->PaperWidth], Cell["Here is what it can do for us!", "Text"], Cell[BoxData[{ \(\(Print["\", exactvalue];\)\), "\[IndentingNewLine]", \(\(Print["\", simpson[f[x], x, a, b, 6] // N];\)\), "\n", \(\(Print["\", trapezoid[f[x], x, a, b, 6] // N];\)\), "\n", \(\(Print["\", riemSumMid[f[x], x, a, b, 6] // N];\)\)}], "Input", PageWidth->PaperWidth], Cell["Now we compare the errors.", "Text", PageWidth->PaperWidth], Cell[BoxData[{ \(\(Print["\", \((simpson[f[x], x, a, b, 6] // N)\) - exactvalue];\)\), "\n", \(\(Print["\", \((trapezoid[ f[x], x, a, b, 6] // N)\) - exactvalue];\)\), "\n", \(\(Print["\", \((riemSumMid[ f[x], x, a, b, 6] // N)\) - exactvalue];\)\)}], "Input", PageWidth->PaperWidth], Cell[CellGroupData[{ Cell["For next week\tHand in next Thursday 1 May 2003:", "Subsubsection"], Cell[TextData[{ "Don't be working on this now, or you'll ", StyleBox["never", FontWeight->"Bold"], " finish the lab by Lunch Time.\nUse", " the functions f1, f2, f3,f4,and f5 defined below, the interval [0, 4], \ and create a table comparing the results and errors of each of the \ integration techniques we've looked at. Then summarize what you observe..." }], "Text"], Cell[BoxData[{ \(\ \(Clear[f1, f2, f3, f4, f5];\)\), "\[IndentingNewLine]", \(\(f1[x_] = x^2;\)\), "\[IndentingNewLine]", \(\(f2[x_] = x^4;\)\), "\[IndentingNewLine]", \(\(f3[x_] = 1/\((x + 1)\);\)\), "\[IndentingNewLine]", \(\(f4[x_] = Sqrt[1 + x^2];\)\), "\[IndentingNewLine]", \(\(f5[x_] = Exp[x];\)\)}], "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Some Theory from Lagrange...", "Section"], Cell["\<\ Get your ZZ-Top beards out and sing along \"...that small town, La \ Grange...\"\ \>", "Commentary", FontColor->RGBColor[0, 0, 1]], Cell[TextData[{ "The Newton-Cotes formulas that are derived in the book come from replacing \ the integrand, ", StyleBox["f,", FontSlant->"Italic"], " ", "by the ", StyleBox["n", FontSlant->"Italic"], "th degree Lagrange polynomial that interpolates the values of ", StyleBox["f", FontSlant->"Italic"], " at evenly spaced points. So... let's see where it comes from:" }], "Text", TextAlignment->Left, TextJustification->0], Cell["Let's find a generic Lagrange polynomial for two points:", "Text", TextAlignment->Left, TextJustification->0], Cell[BoxData[{ \(Clear[data, \ P1, x0, \ x1, fx0, fx1]\), "\[IndentingNewLine]", \(\(data = {{x0, fx0}, {x1, fx1}};\)\), "\[IndentingNewLine]", \(P1[x_] = InterpolatingPolynomial[data, x]\)}], "Input"], Cell["\<\ Now, what happens if we integrate this on the interval [x0, x1]?\ \>", "Text"], Cell[BoxData[ \(Integrate[P1[x], {x, x0, x1}]\)], "Input"], Cell[TextData[{ "Hmmmm.... according to the book, this should be the Trapezoid Rule: ", Cell[BoxData[ \(TraditionalForm\`\(\(x1 - x0\)\/2\) \((fx0 + fx1)\)\)]], ". " }], "Text"], Cell[TextData[StyleBox["Question:\tCan you simplify it to look that way?", FontWeight->"Bold"]], "Text", Background->RGBColor[0, 1, 1]], Cell[BoxData[""], "Input"], Cell[CellGroupData[{ Cell["Simpson's:", "Subsubsection"], Cell["What happens if we have more data points?", "Text"], Cell[BoxData[{ \(Clear[data, \ P2, x0, \ x1, x2, fx0, fx1, fx2]\), "\[IndentingNewLine]", \(\(data = {{x0, fx0}, {x1, fx1}, {x2, fx2}};\)\), "\[IndentingNewLine]", \(P2[x_] = InterpolatingPolynomial[data, x]\)}], "Input"], Cell["\<\ Now, what happens if we integrate this on the interval [a,b]?\ \>", "Text"], Cell[BoxData[ \(Integrate[P2[x], {x, x0, x2}] // Simplify\)], "Input"], Cell[TextData[{ "Hmmmm.... according to the book, this should be Simpson's 1/3 Rule: ", Cell[BoxData[ \(TraditionalForm\`\(h\/3\) \((fx0 + 4 fx1 + fx2)\)\)]], ". " }], "Text"], Cell[TextData[StyleBox["Question:\tCan you simplify it to look that way?", FontWeight->"Bold"]], "Text", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "Hint: what do you know about the width of the intervals [x0, x1] and [x1, \ x2]?\nHint 2: Look at the optional forms of the ", StyleBox["Simplify", FontWeight->"Bold"], " command... maybe try something like: Simplify[whatever, {x0-x1==x1-x2}] \ this tells ", StyleBox["Mathematica", FontSlant->"Italic"], " that when simplifying, it can treat the quantities x0-x1 and x1-x2 as the \ same value..." }], "Text", FontColor->RGBColor[0, 0, 1]], Cell[BoxData[ \(\(Yup ... \)\ You' ve\ got\ to\ do\ some\ typing\ in\ here ... \)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["You Try it!", "Subsubsection"], Cell[TextData[{ "So... what happens if we do this with 4 data points? Can you end up with \ the Simpson's 3/8 rule:\n\t\t\t\t\t", Cell[BoxData[ \(TraditionalForm\`\(\(3 h\)\/8\) \((fx0 + 3 fx1 + 3 fx2 + fx3)\)\)]], "?" }], "Text"], Cell[BoxData[""], "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Error... Not Computing...", "Section"], Cell[TextData[{ "In reality... we don't use this stuff much on functions we can ", StyleBox["actually", FontWeight->"Bold"], " integrate. So... what if we are controlling a satellite's orbit... if we \ know the degree of accuracy we desire, what technique (or how many intervals) \ do we need? Let's look at the ", StyleBox["theoretical", FontWeight->"Bold"], " error bounds:", " " }], "Text"], Cell[CellGroupData[{ Cell["Trapezoids...", "Subsubsection"], Cell[TextData[{ "In the last section, we saw that the Trapezoid rule, ", Cell[BoxData[ \(TraditionalForm\`\(\(x1 - x0\)\/2\) \((fx0 + fx1)\)\)]], ", was really just the integral of the degree 1 Interpolating Polynomial:" }], "Text", TextAlignment->Left, TextJustification->0], Cell[BoxData[{ \(Clear[data, \ P1, x0, \ x1, fx0, fx1]\), "\[IndentingNewLine]", \(\(data = {{x0, fx0}, {x1, fx1}};\)\), "\[IndentingNewLine]", \(\(P1[x_] = InterpolatingPolynomial[data, x];\)\), "\[IndentingNewLine]", \(Simplify[Integrate[P1[x], {x, x0, x1}]]\)}], "Input"], Cell[TextData[{ "So, it might make sense to say that the ", StyleBox["error", FontWeight->"Bold"], " from this integration is just the error of the interpolating polynomial, \ which is related back to the error term from a first order Taylor Series \ (remember the first couple of fun labs!?!?!). This error is: ", Cell[BoxData[ \(TraditionalForm\`\(\(f\^\((2)\)\)( z)\) \((\(\((x - x0)\) \((x - x1)\)\)\/2)\)\)]], ", where ", StyleBox["z", FontSlant->"Italic"], " is some value in the interval (x0, x1) that maximizes the absolute value \ of ", StyleBox["f''", FontSlant->"Italic"], ". Let's let ", StyleBox["Mathematica", FontSlant->"Italic"], " do that for us:" }], "Text"], Cell[BoxData[{ \(Clear[data, \ E1, x0, \ x1, fx0, fx1, \ fdpz]\), "\[IndentingNewLine]", \(E1[x_] = Simplify[Integrate[ fdpz*\(\((x - x0)\) \((x - x1)\)\)\/2, {x, x0, x1}]]\)}], "Input"], Cell[TextData[{ "In my warped notation, ", StyleBox["fdpz", FontWeight->"Bold"], " is ", Cell[BoxData[ \(TraditionalForm\`\(f\^\((2)\)\)(z)\)]], "." }], "Text", FontColor->RGBColor[0, 0, 1]], Cell[TextData[{ "Thus, the error for using a ", StyleBox["single", FontWeight->"Bold"], " trapezoid to integrate over the interval [x0, x1] is exactly as we see in \ the book (box 21.2, page 589):\n\t ", Cell[BoxData[ \(TraditionalForm\`\(\(-1\)\/12\) \(\((x1 - x0)\)\^3\) \(\(f\^\((2)\)\)( z)\)\)]], ", \nand if we let ", Cell[BoxData[ \(TraditionalForm\`h = \((x1 - x0)\)\)]], ", then we can say that the error is ", StyleBox["O", FontWeight->"Bold", FontSlant->"Italic"], Cell[BoxData[ \(TraditionalForm\`\((h)\)\^3\)], FontWeight->"Bold"], "." }], "Text"], Cell[TextData[{ "If we use the ", StyleBox["composite ", FontWeight->"Bold"], "(or", StyleBox[" Multiple-Application", FontWeight->"Bold"], ")", StyleBox[" ", FontWeight->"Bold"], "trapezoid rule - where we integrate over the interval [a, b] and break up \ that interval into several subintervals, each of which we apply a single \ trapezoid rule to integrate - we end up with an error term of: ", StyleBox["O", FontSlant->"Italic"], Cell[BoxData[ \(TraditionalForm\`\((h)\)\^2\)]], ". What is going on?\nWell, the ", StyleBox["global", FontWeight->"Bold"], " error that we have when applying the trapezoid rule multiple times is: \ ", Cell[BoxData[ \(TraditionalForm\`\(\(-1\)\/12\) \(\((x1 - x0)\)\^3\) \(\(f\^\((2)\)\)( z1)\) + \(\(-1\)\/12\) \(\((x1 - x0)\)\^3\) \(\(f\^\((2)\)\)( z2)\) + \[CenterEllipsis] + \(\(-1\)\/12\) \(\((x1 - x0)\)\^3\) \(\(f\^\((2)\)\)(zn)\)\)]], ", which we can rewrite as ", Cell[BoxData[ \(TraditionalForm\`\(\(-1\)\/12\) \(\((x1 - x0)\)\^3\)[\(f\^\((2)\)\)( z1) + \(f\^\((2)\)\)(z2) + \[CenterEllipsis] + \(f\^\((2)\)\)( zn)]\)]], ". If we further assume that ", StyleBox["f''", FontSlant->"Italic"], " is continuous (which it should be if we are to be using this \ technique...), then the intermediate value theorem says there is some point, \ ", StyleBox["z", FontSlant->"Italic"], ", for which ", Cell[BoxData[ \(TraditionalForm\`n*\(\(f\^\((2)\)\)(z)\) = \(f\^\((2)\)\)( z1) + \(f\^\((2)\)\)(z2) + \[CenterEllipsis] + \(f\^\((2)\)\)( zn)\)]], ". But if we have ", StyleBox["n", FontSlant->"Italic"], " of these intervals that we are subdividing [a, b] into, then ", Cell[BoxData[ \(TraditionalForm\`n*\((x1 - x0)\) = b - a\)]], ". So we can rewrite the error as:\n", Cell[BoxData[ \(TraditionalForm\`\(\(-1\)\/12\) \((b - a)\) \(\((x1 - x0)\)\^2\) \(\(f\^\((2)\)\)(z)\)\^\((2)\)\)]], ", which is an ", StyleBox["O", FontSlant->"Italic"], Cell[BoxData[ \(TraditionalForm\`\((h)\)\^2\)]], " error! \n", StyleBox["note: ", FontWeight->"Bold"], "here, ", StyleBox["h", FontSlant->"Italic"], " should be ", StyleBox["much smaller", FontWeight->"Bold"], " than it is above... we control the interval width ", StyleBox["h", FontSlant->"Italic"], " by choosing the number of intervals ", StyleBox["n", FontSlant->"Italic"], ", which is related to the number of data points we sample, ", StyleBox["n+1", FontSlant->"Italic"], "." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Simpson's 1/3 rule:", "Subsubsection"], Cell[TextData[{ "Recall that Simpson's 1/3 rule is ", Cell[BoxData[ \(TraditionalForm\`\(\(x1 - x0\)\/3\) \((fx0 + 4 fx1 + fx2)\)\)]], ":" }], "Text"], Cell[BoxData[{ \(Clear[data, \ P2, x0, \ x1, x2, fx0, fx1, fx2]\), "\[IndentingNewLine]", \(\(data = {{x0, fx0}, {x1, fx1}, {x2, fx2}};\)\), "\[IndentingNewLine]", \(\(P2[x_] = InterpolatingPolynomial[data, x];\)\), "\[IndentingNewLine]", \(Simplify[ Integrate[ P2[x], {x, x0, x2}], {\((x2 - x1)\) \[Equal] \((x1 - x0)\)}]\)}], "Input"], Cell[TextData[{ "Since the error term for a degree 2 interpolating polynomial is ", Cell[BoxData[ \(TraditionalForm\`\(\(f\^\((3)\)\)( z)\) \((\(\((x - x0)\) \((x - x1)\) \((x - x2)\)\)\/6)\)\)]], ", what happens if we try to duplicate the above argurment to obtain the \ error term for Simpson's 1/3 rule (see box 21.3 on page 597):" }], "Text"], Cell[BoxData[{ \(Clear[\ E2, x0, \ x1, x2, fx0, fx1, fx2, \ ftpz]\), "\[IndentingNewLine]", \(E2[x_] = Simplify[Integrate[ ftpz*\(\((x - x0)\) \((x - x1)\) \((x - x2)\)\)\/6, {x, x0, x2}], {x3 - x2 \[Equal] x1 - x0, x2 - x1 \[Equal] x1 - x0}]\)}], "Input"], Cell[TextData[{ "Say what? No Error... ", StyleBox["NO WAY!\n", FontWeight->"Bold"], "Looks like the third order term drops out.... hmmm... cool. The generic \ formula for a degree ", StyleBox["n", FontSlant->"Italic"], " interpolating polynomial is: ", " ", Cell[BoxData[ \(TraditionalForm\`\(\(f\^\((n + 1)\)\)( z)\) \((\(\((x - x0)\) \((x - x1)\) \((x - x2)\) \ \(\[CenterEllipsis](x - xn)\)\)\/\(\((n + 1)\)!\))\)\)]], ". What if we check out the fourth order error term, but only integrate \ from x0 to x3..." }], "Text"], Cell[BoxData[ \(Once\ again, \ \ type\ something\ here, \(\(anything ... \)\ and\ then\ \ hit\ NUM\ ENTER ... \)\ and\ watch\ the\ errors\ fly ... \)], "Input"], Cell[TextData[{ "In my warped notation, ", StyleBox["ffpz", FontWeight->"Bold"], " is ", Cell[BoxData[ \(TraditionalForm\`\(f\^\((4)\)\)(z)\)]], "." }], "Text", FontColor->RGBColor[0, 0, 1]] }, Closed]], Cell[CellGroupData[{ Cell["Simpson's 3/8 rule:", "Subsubsection"], Cell["\<\ So... what exactly does including that extra data point buy for us in terms \ of accuracy? Here is Simpson's 3/8 rule:\ \>", "Text"], Cell[BoxData[{ \(Clear[data, \ P3, x0, \ x1, x2, x3, fx0, fx1, fx2, fx3]\), "\[IndentingNewLine]", \(\(data = {{x0, fx0}, {x1, fx1}, {x2, fx2}, {x3, fx3}};\)\), "\[IndentingNewLine]", \(\(P3[x_] = InterpolatingPolynomial[data, x];\)\), "\[IndentingNewLine]", \(Simplify[ Integrate[P3[x], {x, x0, x3}], {x2 - x1 \[Equal] x1 - x0, x3 - x2 \[Equal] x1 - x0}]\)}], "Input"], Cell["... and the error term for it:", "Text"], Cell[BoxData[{ \(Clear[\ E3, x0, \ x1, x2, x3, \ ftpz]\), "\[IndentingNewLine]", \(E3[x_] = Simplify[Integrate[ ftpz*\(\((x - x0)\) \((x - x1)\) \((x - x2)\) \((x - x3)\)\)\/24, \ {x, x0, x3}], {x3 - x2 \[Equal] x1 - x0, x2 - x1 \[Equal] x1 - x0}]\)}], "Input"], Cell[TextData[{ "The reason that you've heard of Simpson's 1/3 rule and probably ", StyleBox["not", FontWeight->"Bold"], " heard of Simpson's 3/8 rule is that for the extra work involved of having \ to work with a degree 3 interpolating polynomial, your error term is on the \ same order of magnitude as Simpson's 1/3 rule. \nWhenever you have an ", StyleBox["even", FontWeight->"Bold"], " degree interpolating polynomial, the error terms will cancel when \ integrating... meaning that you have to treat them as \"the next degree \ higher\"... so an integration scheme based upon a degree 4 interpolating \ polynomial will have error bounds similar to a degree 5 interpolating \ polynomial... ", StyleBox["cool", FontWeight->"Bold"], ". ", " " }], "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["How does Mathematica Integrate?", "Section"], Cell[TextData[{ "How does ", StyleBox["Mathematica", FontSlant->"Italic"], " perform integration. Let's look at the following example from \ Statistics... For any positive \[Sigma], it can be shown that \[Integral]", Cell[BoxData[ \(TraditionalForm\`\(1\/\(\[Sigma] \@\( 2 \[Pi]\)\)\) \(e\^\(\(-\((1/ 2)\)\) \((x/\[Sigma])\)\^2\)\) dx = 1\)]], "(the integral is from negative to positive infinity). The function ", Cell[BoxData[ \(TraditionalForm\`f(x) = \(1\/\(\[Sigma] \@\( 2 \[Pi]\)\)\) e\^\(\(-\((1/2)\)\) \((x/\[Sigma])\)\^2\)\)]], " is the normal density function wtih mean 0 and standard deviation \ \[Sigma]. The probability that a randomly chosen value described by this \ distribution lies in [a,b] is simply the integral of ", StyleBox["f", FontSlant->"Italic"], " over the interval. What is the probability that a randomly chosen value \ will lie in [-2\[Sigma], 2\[Sigma]] or [-3\[Sigma], 3\[Sigma]]?\t \nFirst, \ try ", StyleBox["Mathematica", FontSlant->"Italic"], "'s integrate command, then try the nintegrate command. Also try using the \ integration formulas from ", StyleBox["A Picture Paints a Thousand Words", FontWeight->"Bold"], " - here they are again:\t" }], "Text", TextAlignment->Left, TextJustification->0], Cell[BoxData[{ \(\(\(riemSumLeft[f_, x_, a_, b_, n_] := Sum[\((f /. x \[Rule] \((a + \((i - 1)\)*\((b - a)\)/n)\))\)*\((b - a)\)/ n, {i, 1, n}];\)\(\n\) \)\), "\n", \(\(\(riemSumRight[f_, x_, a_, b_, n_] := Sum[\((f /. x \[Rule] \((a + i*\((b - a)\)/n)\))\)*\((b - a)\)/n, {i, 1, n}];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(trapezoid[f_, x_, a_, b_, n_] := \((riemSumLeft[f, x, a, b, n] + riemSumRight[f, x, a, b, n])\)/2;\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(\(riemSumMid[f_, x_, a_, b_, n_] := Sum[\((f /. x \[Rule] \((a + \((i - 1/2)\)*\((b - a)\)/n)\))\)*\((b - a)\)/n, {i, 1, n}];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(simpson[f_, x_, a_, b_, n_] := \((2*riemSumMid[f, x, a, b, n/2] + trapezoid[f, x, a, b, n/2])\)/3;\)\)}], "Input", PageWidth->PaperWidth], Cell["Remember, here is how they are used:", "Text"], Cell[BoxData[{ \(\(Print["\", exactvalue];\)\), "\[IndentingNewLine]", \(\(Print["\", simpson[f[x], x, a, b, 6] // N];\)\), "\n", \(\(Print["\", trapezoid[f[x], x, a, b, 6] // N];\)\), "\n", \(\(Print["\", riemSumMid[f[x], x, a, b, 6] // N];\)\)}], "Input", PageWidth->PaperWidth], Cell[TextData[StyleBox["Question:\tWhich technique gives you the best answer \ - how does Mathematica compare?", FontWeight->"Bold"]], "Text", Background->RGBColor[0, 1, 1]] }, Closed]] }, FrontEndVersion->"4.1 for Microsoft Windows", ScreenRectangle->{{0, 1024}, {0, 677}}, AutoGeneratedPackage->None, WindowToolbars->{"RulerBar", "EditBar"}, CellGrouping->Manual, WindowSize->{1016, 644}, WindowMargins->{{0, Automatic}, {Automatic, 0}}, StyleDefinitions -> "DemoText.nb" ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[1705, 50, 74, 1, 70, "Title"], Cell[1782, 53, 41, 0, 41, "Subtitle"], Cell[1826, 55, 5804, 108, 19, "Input", CellOpen->False, InitializationCell->True], Cell[CellGroupData[{ Cell[7655, 167, 31, 0, 54, "Section"], Cell[7689, 169, 286, 8, 48, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[8012, 182, 146, 4, 51, "Section"], Cell[8161, 188, 1084, 28, 160, "Text"], Cell[9248, 218, 181, 5, 23, "Commentary"], Cell[9432, 225, 252, 5, 100, "Input"], Cell[9687, 232, 281, 6, 91, "Text"], Cell[9971, 240, 652, 25, 86, "Text"], Cell[10626, 267, 398, 9, 120, "Input"], Cell[11027, 278, 38, 0, 29, "Text"], Cell[11068, 280, 896, 18, 386, "Input"], Cell[11967, 300, 54, 0, 29, "Text"], Cell[12024, 302, 186, 4, 60, "Input"], Cell[12213, 308, 88, 1, 29, "Text"], Cell[12304, 311, 240, 5, 100, "Input"], Cell[12547, 318, 67, 1, 29, "Text"], Cell[12617, 321, 252, 4, 80, "Input"], Cell[12872, 327, 249, 6, 29, "Text"], Cell[13124, 335, 293, 4, 64, "Text"], Cell[13420, 341, 230, 5, 80, "Input"], Cell[13653, 348, 29, 0, 29, "Text"], Cell[13685, 350, 307, 6, 120, "Input"], Cell[13995, 358, 70, 1, 29, "Text"], Cell[14068, 361, 345, 6, 100, "Input"], Cell[14416, 369, 265, 5, 45, "Text"], Cell[14684, 376, 610, 15, 67, "Text"], Cell[15297, 393, 186, 4, 80, "Input"], Cell[15486, 399, 46, 0, 29, "Text"], Cell[15535, 401, 507, 9, 160, "Input"], Cell[16045, 412, 67, 1, 29, "Text"], Cell[16115, 415, 471, 7, 140, "Input"], Cell[CellGroupData[{ Cell[16611, 426, 73, 0, 44, "Subsubsection"], Cell[16687, 428, 382, 8, 79, "Text"], Cell[17072, 438, 345, 6, 140, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[17466, 450, 47, 0, 34, "Section"], Cell[17516, 452, 142, 4, 23, "Commentary"], Cell[17661, 458, 455, 15, 48, "Text"], Cell[18119, 475, 119, 2, 29, "Text"], Cell[18241, 479, 215, 3, 80, "Input"], Cell[18459, 484, 88, 2, 29, "Text"], Cell[18550, 488, 62, 1, 40, "Input"], Cell[18615, 491, 191, 5, 32, "Text"], Cell[18809, 498, 139, 2, 45, "Text"], Cell[18951, 502, 26, 0, 40, "Input"], Cell[CellGroupData[{ Cell[19002, 506, 35, 0, 44, "Subsubsection"], Cell[19040, 508, 57, 0, 29, "Text"], Cell[19100, 510, 255, 5, 80, "Input"], Cell[19358, 517, 85, 2, 29, "Text"], Cell[19446, 521, 74, 1, 40, "Input"], Cell[19523, 524, 190, 5, 32, "Text"], Cell[19716, 531, 139, 2, 45, "Text"], Cell[19858, 535, 481, 12, 79, "Text"], Cell[20342, 549, 105, 2, 40, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[20484, 556, 36, 0, 28, "Subsubsection"], Cell[20523, 558, 259, 7, 63, "Text"], Cell[20785, 567, 26, 0, 40, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[20860, 573, 44, 0, 34, "Section"], Cell[20907, 575, 415, 11, 48, "Text"], Cell[CellGroupData[{ Cell[21347, 590, 38, 0, 44, "Subsubsection"], Cell[21388, 592, 291, 7, 32, "Text"], Cell[21682, 601, 302, 5, 100, "Input"], Cell[21987, 608, 735, 21, 70, "Text"], Cell[22725, 631, 219, 5, 73, "Input"], Cell[22947, 638, 214, 9, 29, "Text"], Cell[23164, 649, 623, 20, 94, "Text"], Cell[23790, 671, 2675, 79, 225, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[26502, 755, 44, 0, 28, "Subsubsection"], Cell[26549, 757, 163, 5, 32, "Text"], Cell[26715, 764, 409, 10, 100, "Input"], Cell[27127, 776, 368, 7, 51, "Text"], Cell[27498, 785, 313, 7, 73, "Input"], Cell[27814, 794, 570, 16, 84, "Text"], Cell[28387, 812, 164, 2, 40, "Input"], Cell[28554, 816, 214, 9, 29, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[28805, 830, 44, 0, 28, "Subsubsection"], Cell[28852, 832, 143, 3, 29, "Text"], Cell[28998, 837, 432, 9, 100, "Input"], Cell[29433, 848, 46, 0, 29, "Text"], Cell[29482, 850, 298, 6, 73, "Input"], Cell[29783, 858, 784, 18, 98, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[30616, 882, 50, 0, 34, "Section"], Cell[30669, 884, 1335, 30, 158, "Text"], Cell[32007, 916, 1024, 22, 360, "Input"], Cell[33034, 940, 52, 0, 29, "Text"], Cell[33089, 942, 507, 9, 160, "Input"], Cell[33599, 953, 177, 3, 45, "Text"] }, Closed]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)