(************** Content-type: application/mathematica ************** CreatedBy='Mathematica 5.2' 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[ 26525, 1005]*) (*NotebookOutlinePosition[ 27592, 1040]*) (* CellTagsIndexPosition[ 27548, 1036]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Power Series", "Title", Evaluatable->False, TextAlignment->Center, AspectRatioFixed->True], Cell[TextData[{ "Sean Mauch\nsean@caltech.edu\n", ButtonBox["http://www.its.caltech.edu/~sean", ButtonData:>{ URL[ "http://www.its.caltech.edu/~sean"], None}, ButtonStyle->"Hyperlink"], "\n", "This work is distributed under the GNU FDL. See ", ButtonBox["license.nb ", ButtonData:>{"license.nb", None}, ButtonStyle->"Hyperlink"], "for details." }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell[CellGroupData[{ Cell["The Series[] Function", "Section"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], " generates power series expansions with the ", StyleBox["Series[]", "Input"], " function. " }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(\(?Series\)\)], "Input"], Cell[TextData[{ "Here is the Taylor series expansion of ", Cell[BoxData[ \(TraditionalForm\`\(sin\ x\)\/x\)]], " about the point ", Cell[BoxData[ \(TraditionalForm\`x = 0\)]], " accurate to terms of the order ", Cell[BoxData[ \(TraditionalForm\`x\^3\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[Sin[x]\/x, {x, 0, 3}]\)], "Input"], Cell[TextData[{ "Here is the Taylor series expansion of ", Cell[BoxData[ \(TraditionalForm\`\[ExponentialE]\^x\)]], " about the point ", Cell[BoxData[ \(TraditionalForm\`x = 1\)]], ", accurate to terms of the order ", Cell[BoxData[ \(TraditionalForm\`\((x - 1)\)\^5\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[Exp[x], {x, 1, 5}]\)], "Input"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], " will also handle Laurent series." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[Cos[x]/x\^2, {x, 0, 2}]\)], "Input"], Cell[TextData[{ "If you ask for a Laurent series with an infinite number of negative terms, \ then ", StyleBox["Mathematica", FontSlant->"Italic"], " will inform you that the point is an essential singularity." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[Exp[1\/x], {x, 0, 2}]\)], "Input"], Cell[TextData[{ "For series expansions about infinity, you get a power series in ", Cell[BoxData[ \(TraditionalForm\`1\/x\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[1\/\(x\^2 + 1\), {x, \[Infinity], 2}]\)], "Input"], Cell[BoxData[ \(Series[x\^3\/\(x + 1\), {x, Infinity, 3}]\)], "Input"], Cell[TextData[{ "When you perform mathematical operations on series, ", StyleBox["Mathematica", FontSlant->"Italic"], " keeps track of the order of accuracy for you." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(f[x_] = Series[Cos[x], {x, 0, 10}]\)], "Input"], Cell[BoxData[ \(g[x_] = Series[Sin[x], {x, 0, 5}]\)], "Input"], Cell[BoxData[ \(f[x] + g[x]\)], "Input"], Cell[BoxData[ \(f[x]/g[x]\)], "Input"], Cell[BoxData[ \(g[x]\^2\)], "Input"], Cell[BoxData[ \(Sin[g[x]]\)], "Input"], Cell[TextData[{ "It is an error to evaluate a series for any value of the argument. ", Cell[BoxData[ \(TraditionalForm\`O[x]\)]], " means \"terms of the order ", Cell[BoxData[ \(TraditionalForm\`x\)]], "\"; it is not a function that can be evaluated." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(f[1]\)], "Input"], Cell[TextData[{ "In order to evaluate a series you must use ", StyleBox["Normal[]", "Input"], " to convert the series to a function." }], "Text", TextAlignment->Left, TextJustification->1], Cell[TextData[{ "Here we plot a polynomial approximation of ", Cell[BoxData[ \(TraditionalForm\`cos\ x\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(nf[x_] = Normal[f[x]]\)], "Input"], Cell[BoxData[ \(Plot[{nf[x], Cos[x]}, {x, \(-2\) \[Pi], 2 \[Pi]}]; \n Clear[g, f, nf]\)], "Input"], Cell[TextData[{ "The ", StyleBox["Coefficient[]", FontWeight->"Bold"], " function is useful in manipulating power series." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(\(?Coefficient\)\)], "Input"], Cell[BoxData[ \(Series[Sin[y[x]], {x, 0, 2}]\)], "Input"], Cell[TextData[{ "The coefficient of the ", Cell[BoxData[ \(TraditionalForm\`x\^2\)]], " term in the above expansion is" }], "Text"], Cell[BoxData[ \(Coefficient[%, x, 2]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Finite Difference Schemes", "Section"], Cell[TextData[{ "In solving differential equations with difference schemes, one constructs \ a sequence of values ", Cell[BoxData[ \(TraditionalForm\`{y\_n}\)]], " that approximate the solution ", Cell[BoxData[ \(TraditionalForm\`y(x)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`y\_n \[TildeEqual] y(n\ h)\)]], ", where ", Cell[BoxData[ \(TraditionalForm\`h\)]], " is the discretization or stepsize." }], "Text", TextJustification->1], Cell["\<\ A first order accurate difference scheme for calculating the first \ derivative is\ \>", "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(TraditionalForm\`f' \((x)\) \[TildeEqual] \(f(x + h) - f(x)\)\/h\)], "DisplayFormula", TextAlignment->Center], Cell["\<\ This is known as the forward difference. Note that if we take the \ limit as the stepsize approaches zero then the above formula becomes the \ definition of the derivative.\ \>", "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm \`f' \((x)\) = lim\+\(h \[Rule] 0\)\ \(f(x + h) - f(x)\)\/h\)], "DisplayFormula", TextAlignment->Center], Cell[TextData[{ "If ", Cell[BoxData[ \(TraditionalForm\`h\)]], " is small, we would expect the forward difference to approximate the first \ derivative. This difference scheme is called first order accurate because \ the error is ", Cell[BoxData[ \(TraditionalForm\`O(h)\)]], ". ", Cell[BoxData[ \(TraditionalForm\`O(h)\)]], ", (pronounced \"big O of h\"), means \"terms no bigger than ", Cell[BoxData[ \(TraditionalForm\`h\)]], ". For example, ", Cell[BoxData[ \(TraditionalForm\`50\ h\)]], ", ", Cell[BoxData[ \(TraditionalForm\`\(f(x)\)\ h\)]], " and ", Cell[BoxData[ \(TraditionalForm\`2\ h\^2\)]], " are all ", Cell[BoxData[ \(TraditionalForm\`O(h)\)]], ". Below we do a series expansion of the scheme to verify that it is first \ order accurate." }], "Text", TextJustification->1], Cell[BoxData[ \(Series[\(f[x + h] - f[x]\)\/h, {h, 0, 1}]\)], "Input"], Cell[TextData[{ "The error is ", Cell[BoxData[ \(TraditionalForm\`\(1\/2\) \(\(f'\)'\)[x]\ h + \[CenterEllipsis]\)]], ", which is ", Cell[BoxData[ \(TraditionalForm\`O(h)\)]], ". A second order scheme for the first derivative is the centered \ difference scheme:" }], "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm \`f' \((x)\) \[TildeEqual] \(f(x + h) - f(x - h)\)\/\(2 h\)\)], "DisplayFormula", TextAlignment->Center], Cell[BoxData[ \(Series[\(f[x + h] - f[x - h]\)\/\(2 h\), {h, 0, 2}]\)], "Input"], Cell[TextData[{ "Note that when we are examining the order of the scheme, there is no loss \ of generality in taking ", Cell[BoxData[ \(TraditionalForm\`x = 0\)]], "." }], "Text", TextJustification->1], Cell[BoxData[ \(Series[\(f[h] - f[\(-h\)]\)\/\(2 h\), {h, 0, 2}]\)], "Input"], Cell["\<\ Now we see if we can find a difference scheme to approximate the \ second derivative. We look for a scheme of the form,\ \>", "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm \`\(f'\)' \((x)\) \[TildeEqual] \(a\_1\ \(f(x + h)\) + a\_0\ \(f(x)\) + a\_\(-1\)\ \(f(x - h)\)\)\/h\^2 \)], "DisplayFormula", TextAlignment->Center], Cell["\<\ and determine the coefficients. First we do a series expansion of \ the scheme.\ \>", "Text", TextJustification->1], Cell[BoxData[ \(Series[ \(a[1] f[h] + a[0] f[0] + a[\(-1\)] f[\(-h\)]\)\/h\^2, {h, 0, 2}] // Simplify\)], "Input"], Cell[TextData[{ "Since we have three parameters, we will be able to place constraints on \ the first three terms in the expansion. Requiring it to be a first order \ approximation of ", Cell[BoxData[ \(TraditionalForm\`\(f'\)' \((x)\)\)]], " gives us three equations for the coefficients of powers of ", Cell[BoxData[ \(TraditionalForm\`h\)]], ". (Note that if the coefficient of the ", Cell[BoxData[ \(TraditionalForm\`h\^\(-1\)\)]], " term vanishes then so will the coefficient of the ", Cell[BoxData[ \(TraditionalForm\`h\)]], " term. Thus we will have a second order accurate scheme.)" }], "Text", TextJustification->1], Cell[BoxData[ \(Solve[{\n\t\tCoefficient[%, h, \(-2\)] == 0, \n\t\t Coefficient[%, h, \(-1\)] == 0, \n\t\t Coefficient[%, h, 0] == \(\(f'\)'\)[0]}, {a[1], a[0], a[\(-1\)]}]\)], "Input"], Cell["Thus our difference scheme is:", "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm \`\(f'\)' \((x)\) \[TildeEqual] \(f(x + h) - 2 \( f(x)\) + f(x - h)\)\/h\^2\)], "DisplayFormula", TextAlignment->Center], Cell["Below we verify that it is second order accurate.", "Text", TextJustification->1], Cell[BoxData[ \(Series[\(f[h] - 2 f[0] + f[\(-h\)]\)\/h\^2, {h, 0, 2}]\)], "Input"], Cell[CellGroupData[{ Cell["Exercise 1", "Subsubsection"], Cell[TextData[{ "What is the most accurate finite difference approximation of ", Cell[BoxData[ \(TraditionalForm\`\(f'\)' \((x)\)\)]], " you can construct with the five points: ", Cell[BoxData[ \(TraditionalForm\`f(x + 2 h)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x + h)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x - h)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x - 2 h)\)]], "?" }], "Text", TextJustification->1] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell["Runge-Kutta Methods", "Section"], Cell["\<\ Runge-Kutta methods are popular for numerically solving \ differential equations. Consider the first order equation:\ \>", "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm\`y' = f(t, y), \ \ \ \ y(t\_0) = \(y\_0 . \)\)], "DisplayFormula", TextAlignment->Center], Cell[TextData[{ "We construct a sequence ", Cell[BoxData[ \(TraditionalForm\`{y\_n}\)]], " that approximates the solution at time ", Cell[BoxData[ \(TraditionalForm\`t\_n = t\_0 + n\ h\)]], ", ", Cell[BoxData[ \(TraditionalForm\`y\_n \[TildeEqual] y(t\_0 + n\ h) = y(t\_n)\)]], ". One scheme for calculating this sequence is defined recursively by:" }], "Text", TextJustification->1], Cell[BoxData[ FormBox[GridBox[{ {\(y\_\(n + 1\) = y\_n + \(1\/2\) \((k\_1 + k\_2)\)\)}, {\(k\_1 = h\ \(f(t\_n, \ y\_n)\), \)}, {\(k\_2 = h\ \(f(t\_n + h, y\_n + k\_1) . \)\)} }, ColumnAlignments->{Left}], TraditionalForm]], "DisplayFormula", TextAlignment->Center], Cell[TextData[{ "This is a two-stage Runge-Kutta method because it involves two evaluations \ of the function ", Cell[BoxData[ \(TraditionalForm\`f\)]], ". The method is second order accurate. That is, if ", Cell[BoxData[ \(TraditionalForm\`y\_n\)]], " had the same value as the exact solution, ", Cell[BoxData[ \(TraditionalForm\`y\_n = y(t\_0 + n\ h)\)]], ", then the difference between ", Cell[BoxData[ \(TraditionalForm\`y\_\(n + 1\)\)]], " and the exact solution ", Cell[BoxData[ \(TraditionalForm\`y(t\_0 + \((n + 1)\) h)\)]], " is ", Cell[BoxData[ \(TraditionalForm\`O(h\^3)\)]], ", (order of ", Cell[BoxData[ \(TraditionalForm\`h\^3\)]], "). This implies that sequence ", Cell[BoxData[ \(TraditionalForm\`{y\_n}\)]], " agrees with the exact solution up to ", Cell[BoxData[ \(TraditionalForm\`O(h\^2)\)]], ", i.e. the error is ", Cell[BoxData[ \(TraditionalForm\`O(h\^3)\)]], "." }], "Text", TextJustification->1], Cell["\<\ For an application of of this Runge-Kutta algorithm we numerically \ solve an Abel equation of the first kind:\ \>", "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm\`y' = sin(t) + y - y\^3, \ \ \ \ \ \ \ \ y(0) = 1. \)], "DisplayFormula", TextAlignment->Center], Cell[TextData[{ "First we solve the equation numerically with the built-in functions on the \ range ", Cell[BoxData[ \(TraditionalForm\`t \[Element] \([0, 10]\)\)]], "." }], "Text", TextJustification->1], Cell[BoxData[{ \(Clear[y]\), \(soln = NDSolve[{\(y'\)[t] == Sin[t] + y[t] - y[t]\^3, y[0] == 1}, y[t], {t, 0, 20}]; \ny[t_] = y[t] /. First[soln]; \n numericalPlot = Plot[y[t], {t, 0, 10}]; \nClear[y, soln]; \)}], "Input"], Cell[TextData[{ "Next we calculate the sequence ", Cell[BoxData[ \(TraditionalForm\`{y\_n}\)]], " with the Runge-Kutta algorithm." }], "Text", TextJustification->1], Cell[BoxData[ \(Clear[y]; \nh = 0.1; \nf[t_, y_] := Sin[t] + y - y\^3; \ny[0] = 1.0; \n y[n_] := \(y[n] = y[n - 1] + 0.5\ h\ \(( f[\((n - 1)\)\ h, y[n - 1]] + f[n\ h, y[n - 1] + h\ f[\((n - 1)\)\ h, y[n - 1]]])\)\); \n rungeKuttaPlot = ListPlot[Table[{n\ h, y[n]}, {n, 0, 100}]]; \n Clear[h, f, y]; \)], "Input"], Cell[TextData[{ "Note the close agreement between ", StyleBox["Mathematica", FontSlant->"Italic"], "'s solution and that obtained with the Runge-Kutta method." }], "Text", TextJustification->1], Cell[BoxData[ \(Show[numericalPlot, rungeKuttaPlot]; \n Clear[numericalPlot, rungeKuttaPlot]; \)], "Input"], Cell[CellGroupData[{ Cell["Exercise 1", "Subsubsection"], Cell["\<\ There are infinitely many second order accurate, two-stage \ Runge-Kutta methods. The general method has the form:\ \>", "Text", TextJustification->1], Cell[BoxData[ FormBox[GridBox[{ {\(y\_\(n + 1\) = y\_n + a\ k\_1 + b\ k\_2, \)}, {\(k\_1 = h\ \(f(t\_n, \ y\_n)\), \)}, {\(k\_2 = h\ \(f(t\_n + \[Alpha]\ h, y\_n + \[Beta]\ k\_1) . \)\)} }, ColumnAlignments->{Left}], TraditionalForm]], "DisplayFormula", TextAlignment->Center], Cell[TextData[{ "Find the constraints on the parameters ", Cell[BoxData[ \(TraditionalForm\`a\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\)]], ", \[Alpha] and \[Beta] so that the scheme is second order accurate. Do \ this by calculating Taylor series expansions of the ", Cell[BoxData[ \(TraditionalForm\`y\_\(n + 1\)\)]], " given by the scheme and ", Cell[BoxData[ \(TraditionalForm\`y(t + h)\)]], ". (Assume that ", Cell[BoxData[ \(TraditionalForm\`y\_n = y(t\_n)\)]], ". I have let ", Cell[BoxData[ \(TraditionalForm\`t\_0 = 0\)]], " and dropped the subscript on ", Cell[BoxData[ \(TraditionalForm\`t\_n\)]], " to simplify the notation.) By requiring that these two series agree up \ to terms of order ", Cell[BoxData[ \(TraditionalForm\`h\^2\)]], " you will obtain constraints on the parameters. For the series of ", Cell[BoxData[ \(TraditionalForm\`y(t + h)\)]], " you will need to use the fact that ", Cell[BoxData[ \(TraditionalForm\`y(t)\)]], " satisfies the differential equation." }], "Text", TextJustification->1] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell["Series Solution of a Differential Equation", "Section"], Cell["Consider the differential equation", "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(TraditionalForm \`\(y'\)' + \(1\/\(cos\ x\)\) y' + \(\[ExponentialE]\^x\) y = 0, \ \ \ \ y(0) = \(y' \((0)\) = 1. \)\)], "DisplayFormula", TextAlignment->Center, TextJustification->0], Cell[TextData[{ "We wish to find the first few terms in the series expansion of the \ solution about the point ", Cell[BoxData[ \(TraditionalForm\`x = 0\)]], ". From the initial conditions, we know that ", Cell[BoxData[ \(TraditionalForm\`y(x) = 1 + x + \[ScriptCapitalO](x\^2)\)]], ". We expand the differential equation in a Taylor series about the \ origin." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[ \(\(y'\)'\)[x] + \(1\/Cos[x]\) \(y'\)[x] + Exp[x]\ y[x], {x, 0, 1}]\)], "Input"], Cell[TextData[{ "We know the values, ", Cell[BoxData[ \(TraditionalForm\`y(0) = \(y' \((0)\) = 1\)\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(serExp = Normal[%] /. {y[0] -> 1, \(y'\)[0] -> 1}\)], "Input"], Cell[TextData[{ "Now we equation the coefficients of the powers of ", Cell[BoxData[ \(TraditionalForm\`x\)]], " to zero in order to determine ", Cell[BoxData[ \(TraditionalForm\`\(y'\)' \((0)\)\)]], " and ", Cell[BoxData[ \(TraditionalForm\`\(\(y'\)'\)' \((0)\)\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Solve[{Coefficient[serExp, x, 0] == 0, Coefficient[serExp, x, 1] == 0}, {\(\(y'\)'\)[0], \(\(\(y'\)'\)'\)[0]}]\)], "Input"], Cell[TextData[{ "The first few terms in the Taylor series expansion about ", Cell[BoxData[ \(TraditionalForm\`x = 0\)]], " of the solution are" }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[y[x], {x, 0, 3}] /. Join[{y[0] -> 1, \(y'\)[0] -> 1}, First[%]]\)], "Input"], Cell[BoxData[ \(Clear[serExp]\)], "Input"], Cell[CellGroupData[{ Cell["Exercise 3", "Subsubsection"], Cell["\<\ Use Taylor series to approximate the solution of the nonlinear \ differential equation,\ \>", "Text"], Cell[BoxData[ \(TraditionalForm\`y' \((x)\) = x\^2 - \(y\^2\)(x), \ \ \ \ y(0) = 1, \)], "DisplayFormula", TextAlignment->Center, TextJustification->0], Cell[TextData[{ "up to the ", Cell[BoxData[ \(TraditionalForm\`x\^5\)]], " term. Hint: Expand ", Cell[BoxData[ \(TraditionalForm\`y' - x\^2 + y\^2\)]], " in a Taylor series and equate powers of ", Cell[BoxData[ \(TraditionalForm\`x\)]], " to zero. Solve for ", Cell[BoxData[ \(TraditionalForm\`\(y\^\((n)\)\)(0)\)]], " and substitute these values into the series expansion of ", Cell[BoxData[ \(TraditionalForm\`y(x)\)]], "." }], "Text", TextJustification->1] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell["Solutions", "Section"], Cell[CellGroupData[{ Cell["Solution 1", "Subsubsection"], Cell[TextData[{ "What is the most accurate finite difference approximation of ", Cell[BoxData[ \(TraditionalForm\`\(f'\)' \((x)\)\)]], " you can construct with the five points: ", Cell[BoxData[ \(TraditionalForm\`f(x + 2 h)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x + h)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x - h)\)]], ", ", Cell[BoxData[ \(TraditionalForm\`f(x - 2 h)\)]], "?" }], "Text", TextJustification->1], Cell["\<\ The general five point difference scheme for the second derivative \ is:\ \>", "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm \`\(f'\)' \((x)\) \[TildeEqual] \(\(a\_2\) \(f(x + 2 h)\) + \(a\_1\) \(f(x + h)\) + \(a\_0\) \(f(x)\) + \(a\_\(-1\)\) \(f(x - h)\) + \(a\_\(-2\)\) \(f(x - 2 h)\)\)\/h\^2\)], "DisplayFormula", TextAlignment->Center], Cell["\<\ Since there are five undetermined coefficients, we will be able to \ satisfy constraints on at least the first five terms in the series expansion \ of this difference scheme. We find the series below.\ \>", "Text", TextJustification->1], Cell[BoxData[ \(Series[ \((a[2] f[2 h] + a[1] f[h] + a[0] f[0] + a[\(-1\)] f[\(-h\)] + a[\(-2\)] f[\(-2\) h])\)/h\^2, {h, 0, 2}] // Simplify\)], "Input"], Cell["\<\ The following constraints determine the coefficients for a third \ order accurate approximation.\ \>", "Text"], Cell[BoxData[ \(Solve[{Coefficient[%, h, \(-2\)] == 0, \n\t\t Coefficient[%, h, \(-1\)] == 0, \n\t\t Coefficient[%, h, 0] == \(\(f'\)'\)[0], \n\t\t Coefficient[%, h, 1] == 0, \n\t\tCoefficient[%, h, 2] == 0}, { a[\(-2\)], a[\(-1\)], a[0], a[1], a[2]}]\)], "Input"], Cell["When we do a series expansion of this scheme,", "Text"], Cell[BoxData[ \(Series[ \((\(-f[2 h]\) + 16 f[h] - 30 f[0] + 16 f[\(-h\)] - f[\(-2\) h]) \)/\((12 h\^2)\), {h, 0, 4}] // Simplify\)], "Input"], Cell["\<\ we see that we are lucky and the scheme is fourth order accurate.\ \ \>", "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm \`\(f'\)' \((x)\) = \(\(-\(f(x + 2 h)\)\) + 16 \( f(x + h)\) - 30 \( f(x)\) + 16 \( f(x - h)\) - f(x - 2 h)\)\/\(12\ h\^2\) + O(h\^4)\)], "DisplayFormula", TextAlignment->Center] }, Closed]], Cell[CellGroupData[{ Cell["Solution 2", "Subsubsection"], Cell[TextData[{ "First we define the scheme and calculate the Taylor series expansion of ", Cell[BoxData[ \(TraditionalForm\`y\_\(n + 1\)\)]], " in powers of ", Cell[BoxData[ \(TraditionalForm\`h\)]], "." }], "Text", TextJustification->1], Cell[BoxData[ \(k1 = h\ f[t, y[t]]; \nk2 = h\ f[t + \[Alpha]\ h, y[t] + \[Beta]\ k1]; \)], "Input"], Cell[BoxData[ \(Series[y[t] + a\ k1 + b\ k2, {h, 0, 2}] // Simplify\)], "Input"], Cell[TextData[{ "Next we find the Taylor series expansion of ", Cell[BoxData[ \(TraditionalForm\`y(t + h)\)]], "," }], "Text"], Cell[BoxData[ \(Series[y[t + h], {h, 0, 2}]\)], "Input"], Cell[TextData[{ "and use the differential equation to write this in terms of ", Cell[BoxData[ \(TraditionalForm\`f\)]], "." }], "Text"], Cell[BoxData[ \(% /. {\(y'\)[t] -> f[t, y[t]], \(\(y'\)'\)[t] -> D[f[t, y[t]], t]} // Simplify\)], "Input"], Cell[TextData[{ "Comparing the Taylor series expansions of ", Cell[BoxData[ \(TraditionalForm\`y\_\(n + 1\)\)]], " and ", Cell[BoxData[ \(TraditionalForm\`y(t + h)\)]], " we see that the parameters must satisfy the constraints:" }], "Text", TextJustification->1], Cell[BoxData[ \(TraditionalForm\`a + b = 1, \ \ \ \ \[Alpha]\ b = 1\/2, \ \ \ \ \[Beta]\ b = 1\/\(2 . \)\)], "DisplayFormula", TextAlignment->Center], Cell[BoxData[ \(Clear[k1, k2]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Solution 3", "Subsubsection"], Cell[TextData[{ "First we expand the differential equation ", Cell[BoxData[ \(TraditionalForm\`y' - x\^2 + y\^2\)]], " in a Taylor series." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(eqnSeries = Series[\(y'\)[x] - x\^2 + y[x]\^2, {x, 0, 4}]\)], "Input"], Cell[TextData[{ "From the initial condition we know that ", Cell[BoxData[ \(TraditionalForm\`y(0) = 1\)]], ". We make a list of equations by equating the terms in the above Taylor \ series to zero." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(coeffEqns = Join[{y[0] == 1}, Table[Coefficient[eqnSeries, x, n] == 0, {n, 0, 4}]]; \ncoeffEqns // TableForm\)], "Input"], Cell[TextData[{ "Now we solve these equations for the derivatives of ", Cell[BoxData[ \(TraditionalForm\`y\)]], " at ", Cell[BoxData[ \(TraditionalForm\`x = 0\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(coeffs = Solve[coeffEqns, {y[0], \(y'\)[0], \(\(y'\)'\)[0], \(\(\(y'\)'\)'\)[0], \(\(\(\(y'\)'\)'\)'\)[0], \(\(\(\(\(y'\)'\)'\)'\)'\)[0]}]\)], "Input"], Cell[TextData[{ "Below is the Taylor series expansion of ", Cell[BoxData[ \(TraditionalForm\`y(x)\)]], "." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(Series[y[x], {x, 0, 5}]\)], "Input"], Cell[TextData[{ "We substitute in the calculated values of ", Cell[BoxData[ \(TraditionalForm\`\(y\^\((n)\)\)(0)\)]], " to obtain an approximation of the solution." }], "Text", TextAlignment->Left, TextJustification->1], Cell[BoxData[ \(% /. coeffs[\([1]\)] // Normal\)], "Input"] }, Closed]] }, Closed]] }, Open ]] }, FrontEndVersion->"5.2 for Microsoft Windows", ScreenRectangle->{{0, 1680}, {0, 963}}, WindowToolbars->"EditBar", CellGrouping->Automatic, WindowSize->{772, 476}, WindowMargins->{{2, Automatic}, {Automatic, 2}}, PrintingPageRange->{Automatic, Automatic}, PrintingOptions->{"PaperSize"->{612, 792}, "PaperOrientation"->"Portrait", "Magnification"->1}, PrivateNotebookOptions->{"ColorPalette"->{RGBColor, 128}}, ShowCellLabel->True, ShowCellTags->False, RenderingOptions->{"ObjectDithering"->True, "RasterDithering"->False}, CharacterEncoding->"XAutomaticEncoding", Magnification->1.5 ] (******************************************************************* 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[CellGroupData[{ Cell[1776, 53, 102, 3, 141, "Title", Evaluatable->False], Cell[1881, 58, 439, 14, 125, "Text", Evaluatable->False], Cell[CellGroupData[{ Cell[2345, 76, 40, 0, 109, "Section"], Cell[2388, 78, 225, 8, 48, "Text"], Cell[2616, 88, 44, 1, 42, "Input"], Cell[2663, 91, 349, 13, 78, "Text"], Cell[3015, 106, 61, 1, 62, "Input"], Cell[3079, 109, 365, 13, 73, "Text"], Cell[3447, 124, 58, 1, 42, "Input"], Cell[3508, 127, 164, 6, 47, "Text"], Cell[3675, 135, 63, 1, 42, "Input"], Cell[3741, 138, 279, 8, 73, "Text"], Cell[4023, 148, 61, 1, 63, "Input"], Cell[4087, 151, 201, 7, 52, "Text"], Cell[4291, 160, 77, 1, 64, "Input"], Cell[4371, 163, 74, 1, 67, "Input"], Cell[4448, 166, 235, 7, 73, "Text"], Cell[4686, 175, 67, 1, 42, "Input"], Cell[4756, 178, 66, 1, 42, "Input"], Cell[4825, 181, 44, 1, 42, "Input"], Cell[4872, 184, 42, 1, 42, "Input"], Cell[4917, 187, 40, 1, 42, "Input"], Cell[4960, 190, 42, 1, 42, "Input"], Cell[5005, 193, 333, 10, 73, "Text"], Cell[5341, 205, 37, 1, 42, "Input"], Cell[5381, 208, 198, 6, 48, "Text"], Cell[5582, 216, 182, 7, 47, "Text"], Cell[5767, 225, 54, 1, 42, "Input"], Cell[5824, 228, 107, 2, 68, "Input"], Cell[5934, 232, 191, 7, 47, "Text"], Cell[6128, 241, 49, 1, 42, "Input"], Cell[6180, 244, 61, 1, 42, "Input"], Cell[6244, 247, 143, 5, 47, "Text"], Cell[6390, 254, 53, 1, 42, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[6480, 260, 44, 0, 62, "Section"], Cell[6527, 262, 478, 16, 73, "Text"], Cell[7008, 280, 153, 5, 47, "Text"], Cell[7164, 287, 134, 3, 54, "DisplayFormula"], Cell[7301, 292, 221, 5, 73, "Text"], Cell[7525, 299, 149, 4, 55, "DisplayFormula"], Cell[7677, 305, 874, 30, 125, "Text"], Cell[8554, 337, 74, 1, 62, "Input"], Cell[8631, 340, 322, 10, 78, "Text"], Cell[8956, 352, 150, 4, 54, "DisplayFormula"], Cell[9109, 358, 85, 1, 62, "Input"], Cell[9197, 361, 215, 7, 73, "Text"], Cell[9415, 370, 82, 1, 62, "Input"], Cell[9500, 373, 168, 4, 73, "Text"], Cell[9671, 379, 208, 5, 54, "DisplayFormula"], Cell[9882, 386, 128, 4, 47, "Text"], Cell[10013, 392, 133, 3, 93, "Input"], Cell[10149, 397, 670, 17, 125, "Text"], Cell[10822, 416, 208, 4, 120, "Input"], Cell[11033, 422, 70, 1, 47, "Text"], Cell[11106, 425, 169, 4, 54, "DisplayFormula"], Cell[11278, 431, 89, 1, 47, "Text"], Cell[11370, 434, 88, 1, 63, "Input"], Cell[CellGroupData[{ Cell[11483, 439, 35, 0, 61, "Subsubsection"], Cell[11521, 441, 542, 21, 73, "Text"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[12112, 468, 38, 0, 62, "Section"], Cell[12153, 470, 165, 4, 73, "Text"], Cell[12321, 476, 130, 3, 35, "DisplayFormula"], Cell[12454, 481, 419, 12, 73, "Text"], Cell[12876, 495, 318, 7, 94, "DisplayFormula"], Cell[13197, 504, 1027, 34, 151, "Text"], Cell[14227, 540, 158, 4, 73, "Text"], Cell[14388, 546, 137, 3, 37, "DisplayFormula"], Cell[14528, 551, 217, 7, 47, "Text"], Cell[14748, 560, 251, 5, 172, "Input"], Cell[15002, 567, 178, 6, 47, "Text"], Cell[15183, 575, 379, 8, 224, "Input"], Cell[15565, 585, 205, 6, 73, "Text"], Cell[15773, 593, 115, 2, 68, "Input"], Cell[CellGroupData[{ Cell[15913, 599, 35, 0, 61, "Subsubsection"], Cell[15951, 601, 163, 4, 73, "Text"], Cell[16117, 607, 330, 7, 82, "DisplayFormula"], Cell[16450, 616, 1135, 35, 177, "Text"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[17634, 657, 61, 0, 62, "Section"], Cell[17698, 659, 97, 2, 47, "Text"], Cell[17798, 663, 216, 5, 53, "DisplayFormula"], Cell[18017, 670, 444, 12, 99, "Text"], Cell[18464, 684, 117, 3, 66, "Input"], Cell[18584, 689, 178, 7, 47, "Text"], Cell[18765, 698, 82, 1, 42, "Input"], Cell[18850, 701, 361, 13, 73, "Text"], Cell[19214, 716, 161, 3, 94, "Input"], Cell[19378, 721, 214, 7, 47, "Text"], Cell[19595, 730, 110, 2, 68, "Input"], Cell[19708, 734, 46, 1, 42, "Input"], Cell[CellGroupData[{ Cell[19779, 739, 35, 0, 61, "Subsubsection"], Cell[19817, 741, 111, 3, 47, "Text"], Cell[19931, 746, 164, 4, 37, "DisplayFormula"], Cell[20098, 752, 516, 18, 73, "Text"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[20663, 776, 28, 0, 62, "Section"], Cell[CellGroupData[{ Cell[20716, 780, 35, 0, 61, "Subsubsection"], Cell[20754, 782, 542, 21, 73, "Text"], Cell[21299, 805, 120, 4, 47, "Text"], Cell[21422, 811, 287, 6, 54, "DisplayFormula"], Cell[21712, 819, 249, 5, 73, "Text"], Cell[21964, 826, 184, 4, 120, "Input"], Cell[22151, 832, 120, 3, 47, "Text"], Cell[22274, 837, 300, 5, 146, "Input"], Cell[22577, 844, 61, 0, 47, "Text"], Cell[22641, 846, 171, 3, 68, "Input"], Cell[22815, 851, 115, 4, 47, "Text"], Cell[22933, 857, 245, 6, 54, "DisplayFormula"] }, Closed]], Cell[CellGroupData[{ Cell[23215, 868, 35, 0, 38, "Subsubsection"], Cell[23253, 870, 262, 9, 47, "Text"], Cell[23518, 881, 108, 2, 68, "Input"], Cell[23629, 885, 84, 1, 42, "Input"], Cell[23716, 888, 138, 5, 47, "Text"], Cell[23857, 895, 60, 1, 42, "Input"], Cell[23920, 898, 147, 5, 47, "Text"], Cell[24070, 905, 118, 2, 68, "Input"], Cell[24191, 909, 286, 9, 73, "Text"], Cell[24480, 920, 160, 3, 53, "DisplayFormula"], Cell[24643, 925, 46, 1, 42, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[24726, 931, 35, 0, 38, "Subsubsection"], Cell[24764, 933, 210, 7, 47, "Text"], Cell[24977, 942, 90, 1, 42, "Input"], Cell[25070, 945, 270, 8, 73, "Text"], Cell[25343, 955, 153, 3, 94, "Input"], Cell[25499, 960, 248, 10, 47, "Text"], Cell[25750, 972, 192, 4, 68, "Input"], Cell[25945, 978, 177, 7, 47, "Text"], Cell[26125, 987, 56, 1, 42, "Input"], Cell[26184, 990, 235, 7, 47, "Text"], Cell[26422, 999, 63, 1, 42, "Input"] }, Closed]] }, Closed]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)