(************** 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[ 28823, 921]*) (*NotebookOutlinePosition[ 30459, 976]*) (* CellTagsIndexPosition[ 30300, 967]*) (*WindowFrame->Normal*) Notebook[{ Cell[TextData[{ "Introduction to ", StyleBox["Mathematica \[MathematicaIcon]", FontSlant->"Italic"] }], "Title"], Cell["3. Differential Equations", "Subtitle"], Cell["\<\ P. S. Cally, School of Mathematical Sciences, Monash \ University\ \>", "Author"], Cell["\<\ Our main tools for solving differential equations are DSolve and \ NDSolve. They are both very extensive and powerful programs, but neither can \ do everything. Nevertheless, they can do a lot \[Ellipsis]\ \>", "Text"], Cell[CellGroupData[{ Cell["Exact Solutions of ODEs", "Section"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], " can find general solutions to many ODEs -- most of those for which exact \ solutions are known" }], "Text"], Cell[BoxData[ \(DSolve[\(y''\)[x] + x\ \(y'\)[x] + y[x] \[Equal] Sin[x], y[x], x]\)], "Input"], Cell["We can also specify boundary conditions", "Text"], Cell[BoxData[ \(sol = DSolve[{\(y''\)[x] + x\ \(y'\)[x] + y[x] \[Equal] Sin[x], y[0] \[Equal] 0, y[1] \[Equal] 0}, y, x]\)], "Input"], Cell[TextData[{ StyleBox["Note: we have requested y and not y[x] as the solution.", FontWeight->"Bold"], " It is usually much more convenient to get back the pure function. Why?\n\ Check the solution:" }], "Text"], Cell[BoxData[ \(\(y''\)[x] + x\ \(y'\)[x] + y[x] /. sol // FullSimplify\)], "Input"], Cell[BoxData[ \({y[0], y[1]} /. sol\)], "Input"], Cell["Plot", "Text"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. sol], {x, 0, 1}];\)\)], "Input"], Cell[BoxData[ \(Clear[sol]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Numerical Solution of ODEs: IVPs", "Section"], Cell[TextData[{ "Not surprisingly, ", StyleBox["Mathematica", FontSlant->"Italic"], " cannot find an exact solution to this horrible nonlinear equation." }], "Text"], Cell[BoxData[ \(DSolve[{\(y''\)[x] + y[x] \[Equal] x\ y[x]\ Exp[\(-y[x]\^2\)], y[0] \[Equal] 1, \(y'\)[0] \[Equal] 0}, y, x]\)], "Input"], Cell["\<\ But numerical solution is easy! It returns an answer in the form of \ an interpolating function object, which may be evaluated wherever we wish (in \ its domain). Timing[\[Ellipsis]] is wrapped around the call to NDSolve in \ this case to show just how quickly the numerical integration can be \ done.\ \>", "Text"], Cell[BoxData[ \(Timing[\(sol = Flatten[NDSolve[{\(y''\)[x] + y[x] \[Equal] x\ y[x]\ Exp[\(-y[x]\^2\)], y[0] \[Equal] 1, \(y'\)[0] \[Equal] 0}, y, {x, 0, 20}]];\)]\)], "Input"], Cell[TextData[{ "The Flatten is not crucial; try it without to see what happens. Here we \ see how to evaluate ", StyleBox["y", FontSlant->"Italic"], " (or its derivatives) at a particular value of ", StyleBox["x", FontSlant->"Italic"], " in the domain." }], "Text"], Cell[BoxData[ \({y[5], \(y'\)[5]} /. sol\)], "Input"], Cell["And this is how we can plot the solution. ", "Text"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. sol], {x, 0, 20}];\)\)], "Input"], Cell[TextData[{ "You can plot ", StyleBox["y'", FontSlant->"Italic"], " just as easily." }], "Text"], Cell[BoxData[ \(\(Plot[Evaluate[\(y'\)[x] /. sol], {x, 0, 20}];\)\)], "Input"], Cell["\<\ Various (unnecessary) stylistic options are included here, just to \ show you what can be done.\ \>", "Text"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. sol], {x, 0, 20}, AxesLabel \[Rule] TraditionalForm /@ {x, y}, TextStyle \[Rule] {FontFamily \[Rule] "\", FontSlant -> "\", FontSize \[Rule] 14, FontWeight \[Rule] "\"}, ImageSize \[Rule] 400, AxesStyle \[Rule] RGBColor[1, 0, 0], Background \[Rule] RGBColor[ .9, .9, 0], PlotStyle \[Rule] {{RGBColor[0, .5, 0], Dashing[{ .06, .02, .005, .02}], Thickness[ .007]}}, PlotLabel -> "\<\[MathematicaIcon] Silly Graph in Mathematica \ \[MathematicaIcon]\>"];\)\)], "Input"], Cell["\<\ NDSolve has a vast number of options, allowing us to choose many \ different numerical methods and accuracies. See the help utility for details.\ \ \>", "Text"], Cell[BoxData[ \(Clear[sol]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Numerical Solution of PDEs: Linear and Nonlinear BVPs", "Section"], Cell[TextData[{ "For linear Boundary Value Problems, ", StyleBox["Mathematica", FontSlant->"Italic"], " uses the Gel'fand\[Hyphen]Lokutsiyevskii chasing method." }], "Text"], Cell[BoxData[ \(bvpsol = NDSolve[{\(y''\)[x] + 25 Cos[x] y[x] \[Equal] 1, y[0] \[Equal] 0, \(y'\)[1] \[Equal] 0}, y, {x, 0, 1}]\)], "Input"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. bvpsol], {x, 0, 1}, PlotRange \[Rule] All];\)\)], "Input"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], " does not have a built-in method to automatically solve nonlinear BVPs." }], "Text"], Cell[BoxData[ \(NDSolve[{\(y''\)[x] \[Equal] y[x]\^2, y[0] \[Equal] 0, y[1] \[Equal] 2}, y, {x, 0, 1}]\)], "Input"], Cell["\<\ However, we can write our own shooting method (or relaxation method \ if you prefer). Here is a typical setup for shooting.\ \>", "Text"], Cell["\<\ We know that one boundary value at x=0 is y=0, so of course we \ adopt that. But we don't know what y'(0) should be. This utility score takes \ a guess at this, y'(0)=dy, and returns a measure of how well we have done in \ hitting our other BC at x=1.\ \>", "Text"], Cell[BoxData[ \(score[dy_?NumberQ] := \((y[1] - 2)\) /. Flatten[NDSolve[{\(y''\)[x] \[Equal] y[x]\^2, y[0] \[Equal] 0, \(y'\)[0] \[Equal] dy}, y, {x, 0, 1}]]\)], "Input"], Cell["Plot score against dy to see where is has a zero (or zeros).", "Text"], Cell[BoxData[ \(\(Plot[score[dy], {dy, \(-70\), 20}];\)\)], "Input"], Cell["\<\ Clearly, dy\[TildeTilde]2 and -39 are what we want. Since the \ equation is nonlinear, it is not surprising that there is more than one \ solution. Home in on them with FindRoot.\ \>", "Text"], Cell[BoxData[ \(dy1 = dy /. FindRoot[score[dy] \[Equal] 0, {dy, \(-39\)}]\)], "Input"], Cell[BoxData[ \(dy2 = dy /. FindRoot[score[dy] \[Equal] 0, {dy, 2}]\)], "Input"], Cell["\<\ Check to make sure that score vanishes (to high accuracy) at these \ values.\ \>", "Text"], Cell[BoxData[ \(score /@ {dy1, dy2}\)], "Input"], Cell["Now shoot one more time, with the \"correct\" y'(0).", "Text"], Cell[BoxData[ \(\(sol1 = NDSolve[{\(y''\)[x] \[Equal] y[x]\^2, y[0] \[Equal] 0, \(y'\)[0] \[Equal] dy1}, y, {x, 0, 1}];\)\)], "Input"], Cell[BoxData[ \(\(sol2 = NDSolve[{\(y''\)[x] \[Equal] y[x]\^2, y[0] \[Equal] 0, \(y'\)[0] \[Equal] dy2}, y, {x, 0, 1}];\)\)], "Input"], Cell["Plot the resulting solutions.", "Text"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. {sol1, sol2}], {x, 0, 1}, PlotStyle \[Rule] {{}, {Hue[0]}}];\)\)], "Input"], Cell["\<\ There are other methods, but shooting is usually quick and \ accurate.\ \>", "Text"], Cell[BoxData[ \(Clear[bvpsol, score, dy1, dy2, sol1, sol2]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Exact Solution of PDEs", "Section"], Cell["\<\ Exact solutions of PDEs are not so common! However, first order \ PDEs can often be solved, primarily because they can be reduced to ODEs. Here \ are a few:\ \>", "Text"], Cell[BoxData[ \(DSolve[\[PartialD]\_x u[x, y] + x \[PartialD]\_y u[x, y] \[Equal] u[x, y], u, {x, y}]\)], "Input"], Cell[BoxData[ \(DSolve[ y \[PartialD]\_x u[x, y] + x \[PartialD]\_y u[x, y] \[Equal] u[x, y], u, {x, y}]\)], "Input"], Cell[BoxData[ \(DSolve[ y \[PartialD]\_x u[x, y] + x \[PartialD]\_y u[x, y] \[Equal] u[x, y]\^2, u, {x, y}]\)], "Input"], Cell[TextData[{ "Here, C[1] represents an ", StyleBox["arbitrary function", FontWeight->"Bold"], ". \nNow to higher order equations. Even the 2D Laplace equation has an \ exact solution (derivable by the method of characteristics, or a cunningly \ chosen change of variables)." }], "Text"], Cell[BoxData[ \(DSolve[\[PartialD]\_\(x, x\)u[x, y] + \[PartialD]\_\(y, y\)u[x, y] \[Equal] 0, u, {x, y}]\)], "Input"], Cell["\<\ Unfortunately, not very useful when boundary conditions are \ imposed. The wave equation\ \>", "Text"], Cell[BoxData[ \(DSolve[\[PartialD]\_\(t, t\)u[x, t] == \[PartialD]\_\(x, x\)u[x, t], u, {x, t}]\)], "Input"], Cell["\<\ This is the well-known D'Alembert solution. The diffusion equation\ \>", "Text"], Cell[BoxData[ \(DSolve[{\[PartialD]\_t u[x, t] == \[PartialD]\_\(x, x\)u[x, t]}, u, {x, t}]\)], "Input"], Cell["Sorry, can't do it!", "Text"], Cell["\<\ Don't expect to find exact solutions of too many second or higher \ order PDEs! And even when you can, they may not be what you want.\ \>", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Numerical Solution of PDEs", "Section"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], "'s numerical PDE scheme is the method of lines. It is therefore unsuitable \ for elliptic equations. It does a good job for wave and diffusion equations \ though." }], "Text"], Cell["\<\ Here is a straightforward application with Mathematica's default \ settings (4th order finite differences in x). It complains about \ accuracy/precision, but returns a solution.\ \>", "Text"], Cell[BoxData[ \(sol4 = NDSolve[{\[PartialD]\_\(t, t\)u[x, t] == \[PartialD]\_\(x, x\)u[x, t], u[0, t] \[Equal] 0, u[1, t] \[Equal] 0, u[x, 0] \[Equal] \(x\^2\) \((1 - x)\), \(\(Derivative[0, 1]\)[u]\)[ x, 0] \[Equal] 0}, u, {x, 0, 1}, {t, 0, 5}]\)], "Input"], Cell["Now let's set up a plotting procedure.", "Text"], Cell[BoxData[ \(plt[t_, soln_, opts___] := Plot[Evaluate[u[x, t] /. soln], {x, 0, 1}, opts, PlotRange \[Rule] {\(- .16\), .16}, PlotLabel \[Rule] ToString[t], Frame \[Rule] True]\)], "Input"], Cell["Here is the solution at t=5.", "Text"], Cell[BoxData[ \(\(plt[5, sol4];\)\)], "Input"], Cell[BoxData[ \(\(Table[plt[t, sol4], {t, 0, 5, .1}];\)\)], "Input"], Cell["\<\ Collapse these by double-clicking on the frame bar at the right \ which encloses them all (i.e., 2nd from left). Now double-click on the \ picture. Voila!\ \>", "Text"], Cell[TextData[{ "Now let's do a 2D example, vibration of a driven square drum. In this \ example I set ", StyleBox["periodic", FontWeight->"Bold"], " boundary conditions on all four sides (can you see that?). Again, we use \ the default ", StyleBox["Mathematica", FontSlant->"Italic"], " settings, and check how long the integration takes." }], "Text"], Cell[BoxData[ RowBox[{"Timing", "[", RowBox[{"solwv", "=", RowBox[{"First", "[", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{"{", RowBox[{\(\[PartialD]\_\(t, t\)u[t, x, y] == \[PartialD]\_\(x, x\)u[t, x, y] + \[PartialD]\_\(y, y\)u[t, x, y] - Sin[u[t, x, y]]\), ",", \(u[0, x, y] == \[ExponentialE]\^\(-\((x\^2 + y\^2)\)\)\), ",", RowBox[{ RowBox[{ SuperscriptBox["u", TagBox[\((1, 0, 0)\), Derivative], MultilineFunction->None], "[", \(0, x, y\), "]"}], "==", "0"}], ",", \(u[t, \(-5\), y] == u[t, 5, y]\), ",", \(u[t, x, \(-5\)] == u[t, x, 5]\)}], "}"}], ",", "u", ",", \({t, 0, 15}\), ",", \({x, \(-5\), 5}\), ",", \({y, \(-5\), 5}\)}], "]"}], "]"}]}], "]"}]], "Input", CellTags->"NDSolve"], Cell["\<\ Pretty slow. Try again with pseudospectral derivatives (which are \ ideally suited to periodic BCs),\ \>", "Text"], Cell[BoxData[ RowBox[{"Timing", "[", RowBox[{"solwv", "=", RowBox[{"First", "[", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{"{", RowBox[{\(\[PartialD]\_\(t, t\)u[t, x, y] \[Equal] \[PartialD]\_\(x, x\)u[t, x, y] + \[PartialD]\_\(y, y\)u[t, x, y] - Sin[u[t, x, y]]\), ",", \(u[0, x, y] \[Equal] \[ExponentialE]\^\(-\((x\^2 + y\^2)\)\)\), ",", RowBox[{ RowBox[{ SuperscriptBox["u", TagBox[\((1, 0, 0)\), Derivative], MultilineFunction->None], "[", \(0, x, y\), "]"}], "\[Equal]", "0"}], ",", \(u[t, \(-5\), y] \[Equal] u[t, 5, y]\), ",", \(u[t, x, \(-5\)] \[Equal] u[t, x, 5]\)}], "}"}], ",", "u", ",", \({t, 0, 15}\), ",", \({x, \(-5\), 5}\), ",", \({y, \(-5\), 5}\), ",", \(Method \[Rule] {"\", \ \ "\" \[Rule] {"\", (*\ Specify\ use\ of\ a\ tensor\ product\ \(\(grid\)\(.\)\ \)\ *) \ DifferenceOrder \[Rule] "\"\ (*\ Specify\ use\ of\ pseudospectral\ approximations\ for\ \ x\ \ and\ \(\(y\)\(.\)\)*) }}\)}], "]"}], "]"}]}], "]"}]], "Input", CellTags->"NDSolve"], Cell["\<\ Wow, what an improvement! The superior spatial smoothness must have \ allowed the temporal integration to use far fewer points. Get set to plot the \ solution.\ \>", "Text"], Cell[BoxData[ \(Clear[pwv]\)], "Input"], Cell[BoxData[ \(pwv[t_, opts___] := Plot3D[Evaluate[u[t, x, y] /. solwv], {x, \(-5\), 5}, {y, \(-5\), 5}, opts, PlotLabel \[Rule] "\" <> ToString[t], PlotRange \[Rule] {\(- .5\), 1}]\)], "Input", CellTags->"NDSolve"], Cell["\<\ Make a movie, and collapse the frames again. Then play by \ double-clicking.\ \>", "Text"], Cell[BoxData[ \(\(Table[ pwv[t, ColorFunctionScaling \[Rule] False], {t, 0, 15, .2}];\)\)], "Input"], Cell["\<\ Isn't that pretty! Here is the last frame drawn without a mesh (if \ you prefer):\ \>", "Text"], Cell[BoxData[ \(\(pwv[15, Mesh \[Rule] False, ColorFunctionScaling \[Rule] False];\)\)], "Input"], Cell["\<\ Let's redo this with rigid boundary conditions (and appropriately \ modified ICs). We will again use pseudospectral derivatives. However, note \ that pseudospectral derivatives should ONLY be used on either (1) a regular \ grid if there are periodic BCs; or (2) a Gauss-Lobatto grid (zeros of \ Chebyshev polynomials) otherwise. If you try to use pseudospectral \ derivatives with say Dirichlet BCs, you are likely to get smashed by the \ dreaded (shudder) Runge Phenomenon (aaagh!). But this is getting too far into \ mathematics for our current purposes. Ask me if you want to know more, or \ read Fornberg or Boyd. Here's a simple utility for calculating a Gauss-Lobatto grid of length L \ starting at x0, with n points.\ \>", "Text"], Cell[BoxData[ \(CGLGrid[x0_, \ L_, \ n_Integer\ /; \ n\ > \ 1]\ := \ x0\ + 1\/2\ L \((1\ - \ \ Cos[\[Pi]\ Range[0, n - 1]/\((n - 1)\)])\)\)], "Input", CellTags->{"b:7.1.2", "ndsg:2.0.1.2"}], Cell["We use the same grid for x and y in this case.", "Text"], Cell[BoxData[ \(xygrid = CGLGrid[\(-5\), 10, 31] // N\)], "Input"], Cell[TextData[{ "Tell NDSolve to use this grid using the ", StyleBox["Coordinates", FontFamily->"Courier"], " option." }], "Text"], Cell[BoxData[ RowBox[{"Timing", "[", RowBox[{"sol0", "=", RowBox[{"First", "[", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{"{", RowBox[{\(\[PartialD]\_\(t, t\)u[t, x, y] \[Equal] \[PartialD]\_\(x, x\)u[t, x, y] + \[PartialD]\_\(y, y\)u[t, x, y] - Sin[u[t, x, y]]\), ",", \(u[0, x, y] \[Equal] Cos[\(\[Pi]\/10\) x] Cos[\(\[Pi]\/10\) y] \[ExponentialE]\^\(-\((x\^2 + y\^2)\)\)\), ",", RowBox[{ RowBox[{ SuperscriptBox["u", TagBox[\((1, 0, 0)\), Derivative], MultilineFunction->None], "[", \(0, x, y\), "]"}], "\[Equal]", "0"}], ",", \(u[t, \(-5\), y] \[Equal] u[t, 5, y] \[Equal] 0\), ",", \(u[t, x, \(-5\)] \[Equal] u[t, x, 5] \[Equal] 0\)}], "}"}], ",", "u", ",", \({t, 0, 15}\), ",", \({x, \(-5\), 5}\), ",", \({y, \(-5\), 5}\), ",", \(Method \[Rule] {"\", \ \ "\" \[Rule] {"\", (*\ Specify\ use\ of\ a\ tensor\ product\ \(\(grid\)\(.\)\ \)\ *) \ DifferenceOrder \[Rule] "\"\ (*\ Specify\ use\ of\ pseudospectral\ approximations\ for\ \ x\ \ and\ \(\(y\)\(.\)\)*) , Coordinates \[Rule] {xygrid, xygrid}}}\)}], "]"}], "]"}]}], "]"}]], "Input", CellTags->"NDSolve"], Cell["Now plot.", "Text"], Cell[BoxData[ \(pwv0[t_] := Plot3D[Evaluate[u[t, x, y] /. sol0], {x, \(-5\), 5}, {y, \(-5\), 5}, PlotLabel \[Rule] "\" <> ToString[t], PlotRange \[Rule] {\(- .5\), 1}]\)], "Input", CellTags->"NDSolve"], Cell[BoxData[ \(\(Table[pwv0[t], {t, 0, 15, .2}];\)\)], "Input"], Cell[BoxData[ \(Clear[plt, pwv, pwv0, xygrid, sol0, solwv, pdesol]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Exercises", "ExerciseMain"], Cell["\<\ Do not look at the solution until you have completed the \ exercise!\ \>", "Text", TextAlignment->Center, FontFamily->"Helvetica", FontSize->18, FontWeight->"Bold", FontColor->RGBColor[1, 0, 0]], Cell[CellGroupData[{ Cell["Airy Equation", "Exercise"], Cell["\<\ The Airy equation y''=x y pops up in many contexts, from asymptotic \ matching to the propagation of wavefronts. Solve it and plot both solutions \ over (-15,5). Do you understand how Ai may arise in modelling \ wavefronts?\ \>", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell["DSolve[y''[x] - x y[x] == 0, y, x]", "Input", CellTags->"S3.5.11"], Cell[BoxData[ \(\(Plot[AiryAi[x], {x, \(-15\), 5}];\)\)], "Input"], Cell[BoxData[ \(\(Plot[AiryBi[x], {x, \(-15\), 5}];\)\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Bessel Equation", "Exercise"], Cell[TextData[{ "Find the exact solutions of Bessel's equation of order", Cell[BoxData[ RowBox[{ StyleBox[" ", FontSlant->"Plain"], RowBox[{ StyleBox["n", FontSlant->"Italic"], ",", StyleBox[" ", FontSlant->"Italic"], \(\(x\^2\) y'' + x\ y' + \((x\^2 - n\^2)\) y = 0\)}]}]]], ". Plot each of the solutions for ", StyleBox["n", FontSlant->"Italic"], "=0,1,2,3,4,5 on ", " 0False]}], SeriesData[ x, 0, {1, 0, Rational[ -1, 4], 0, Rational[ 1, 64], 0, Rational[ -1, 2304], 0, Rational[ 1, 147456]}, 0, 10, 1], Editable->False]]], " to get started. Why?" }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell["\<\ We need to step away from the singularity at x=0 in order to start \ our numerical integration.\ \>", "Text"], Cell[BoxData[ \(ys[x_] := 1 - x\^2\/4 + x\^4\/64 - x\^6\/2304 + x\^8\/147456\)], "Input"], Cell[BoxData[ \(With[{\[Epsilon] = 0.0001}, bess = NDSolve[{x\ \(y''\)[x] + \ \(y'\)[x] + x\ y[x] \[Equal] 0, y[\[Epsilon]] \[Equal] ys[\[Epsilon]], \(y'\)[\[Epsilon]] \[Equal] \ \(ys'\)[\[Epsilon]]}, y, {x, \[Epsilon], 15}]]\)], "Input"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. bess], {x, 0.0001, 15}];\)\)], "Input"], Cell["Try it without the \[Epsilon].", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Linear BVP", "Exercise"], Cell[TextData[{ "Numerically solve ", Cell[BoxData[ \(TraditionalForm\`y'' + \(\[ExponentialE]\^x\) y = 0\)]], ", with BCs y(0)=0, y(1)+y'(1)=1. Plot The solution." }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell[BoxData[ \(lbvp = NDSolve[{\(y''\)[x] + \ Exp[x] y[x] \[Equal] 0, y[0] \[Equal] 0, y[1] + \(y'\)[1] \[Equal] 1}, y, {x, 0, 5}]\)], "Input"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. lbvp], {x, 0, 5}];\)\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Nonlinear BVP", "Exercise"], Cell[TextData[{ "Numerically solve ", Cell[BoxData[ \(TraditionalForm\`y'' + \(\[ExponentialE]\^x\) y = sin\ y\)]], ", with BCs y(0)=0, y(5)+y'(5)=1. Plot The solution." }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell[BoxData[ \(score[dy_?NumberQ] := \((y[5] + \(y'\)[5] - 1)\) /. Flatten[NDSolve[{\(y''\)[x] + \ Exp[x] y[x] \[Equal] Sin[y[x]], y[0] \[Equal] 0, \(y'\)[0] \[Equal] dy}, y, {x, 0, 5}]]\)], "Input"], Cell[BoxData[ \(\(Plot[score[dy], {dy, \(-1\), 1}];\)\)], "Input"], Cell[BoxData[ \(dy1 = dy /. FindRoot[score[dy] \[Equal] 0, {dy, \(- .3\)}]\)], "Input"], Cell[BoxData[ \(\(sol1 = NDSolve[{\(y''\)[x] + \ Exp[x] y[x] \[Equal] Sin[y[x]], y[0] \[Equal] 0, \(y'\)[0] \[Equal] dy1}, y, {x, 0, 5}];\)\)], "Input"], Cell[BoxData[ \(\(Plot[Evaluate[y[x] /. sol1], {x, 0, 5}];\)\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Annular Drum", "Exercise"], Cell[TextData[{ "An annular drumhead, 0.1"Italic"], " and \[Theta] with 15 points in each direction. (Ignore a warning message \ about \"Scaled local spatial error\".) Make an animation of your solution." }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell["\<\ Since the solution is obviously periodic in \[Theta], we use a \ regular grid in that direction.\ \>", "Text"], Cell[BoxData[ \(\(\[Theta]grid = Range[0, 2 \[Pi], \(2 \[Pi]\)\/15];\)\)], "Input"], Cell[TextData[{ "However, the ", StyleBox["r", FontSlant->"Italic"], "-grid must be Gauss-Lobatto:" }], "Text"], Cell[BoxData[ \(rgrid = CGLGrid[ .1, .9, 15] // N\)], "Input"], Cell["Solve the PDE.", "Text"], Cell[BoxData[ RowBox[{"Timing", "[", RowBox[{"solp", "=", RowBox[{"First", "[", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{"{", RowBox[{\(\[PartialD]\_\(t, t\)u[t, r, \[Theta]] \[Equal] \[PartialD]\_\(r, r\)u[t, r, \[Theta]] + \(1\/r\) \[PartialD]\_r u[t, r, \[Theta]] + \ \(1\/r\^2\) \[PartialD]\_\(\ \[Theta], \[Theta]\)u[t, r, \[Theta]]\), ",", \(u[0, r, \[Theta]] \[Equal] 20 \((r - .1)\) \((1 - r)\) Exp[Sin[\[Theta]]]\), ",", RowBox[{ RowBox[{ SuperscriptBox["u", TagBox[\((1, 0, 0)\), Derivative], MultilineFunction->None], "[", \(0, r, \[Theta]\), "]"}], "\[Equal]", "0"}], ",", \(u[t, .1, \[Theta]] \[Equal] u[t, 1, \[Theta]] \[Equal] 0\), ",", \(u[t, r, 0] \[Equal] u[t, r, 2 \[Pi]]\)}], "}"}], ",", "u", ",", \({t, 0, 5}\), ",", \({r, .1, 1}\), ",", \({\[Theta], 0, 2 \[Pi]}\), ",", \(Method \[Rule] {"\", \ "\" \[Rule] {"\", DifferenceOrder \[Rule] "\"\ , \ \[IndentingNewLine]Coordinates \[Rule] {rgrid, \[Theta]grid}\ \[IndentingNewLine]}}\)}], "]"}], "]"}]}], "]"}]], "Input", CellTags->"NDSolve"], Cell["Plotting utility:", "Text"], Cell[BoxData[ \(p[t_] := ParametricPlot3D[ Evaluate[{r\ Cos[\[Theta]], r\ Sin[\[Theta]], u[t, r, \[Theta]]} /. solp], {r, .1, 1}, {\[Theta], 0, 2 \[Pi]}, PlotLabel \[Rule] "\" <> ToString[t], PlotRange \[Rule] {\(-12\), 12}, Boxed \[Rule] False, Axes \[Rule] False, ImageSize \[Rule] 600, BoxRatios \[Rule] {1, 1, .4}]\)], "Input", CellTags->"NDSolve"], Cell["Make athe animation:", "Text"], Cell[BoxData[ \(\(Table[p[t], {t, 0, 5, .2}];\)\)], "Input"] }, Closed]] }, Open ]] }, FrontEndVersion->"5.2 for X", ScreenRectangle->{{0, 1280}, {0, 1024}}, WindowSize->{817, 935}, WindowMargins->{{1, Automatic}, {Automatic, 2}}, ShowSelection->True, Magnification->1, StyleDefinitions -> "Classroom.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->{ "NDSolve"->{ Cell[12590, 426, 1107, 24, 128, "Input", CellTags->"NDSolve"], Cell[13827, 457, 1560, 32, 193, "Input", CellTags->"NDSolve"], Cell[15622, 500, 254, 5, 63, "Input", CellTags->"NDSolve"], Cell[17594, 558, 1705, 34, 204, "Input", CellTags->"NDSolve"], Cell[19330, 596, 239, 5, 63, "Input", CellTags->"NDSolve"], Cell[26624, 867, 1593, 32, 207, "Input", CellTags->"NDSolve"], Cell[28256, 903, 432, 9, 79, "Input", CellTags->"NDSolve"]}, "b:7.1.2"->{ Cell[17080, 540, 230, 4, 64, "Input", CellTags->{"b:7.1.2", "ndsg:2.0.1.2"}]}, "ndsg:2.0.1.2"->{ Cell[17080, 540, 230, 4, 64, "Input", CellTags->{"b:7.1.2", "ndsg:2.0.1.2"}]}, "S3.5.11"->{ Cell[20395, 640, 74, 1, 47, "Input", CellTags->"S3.5.11"]} } *) (*CellTagsIndex CellTagsIndex->{ {"NDSolve", 29474, 939}, {"b:7.1.2", 29992, 954}, {"ndsg:2.0.1.2", 30101, 957}, {"S3.5.11", 30205, 960} } *) (*NotebookFileOutline Notebook[{ Cell[1754, 51, 122, 4, 70, "Title"], Cell[1879, 57, 45, 0, 41, "Subtitle"], Cell[1927, 59, 91, 3, 39, "Author"], Cell[2021, 64, 228, 4, 46, "Text"], Cell[CellGroupData[{ Cell[2274, 72, 42, 0, 64, "Section"], Cell[2319, 74, 179, 5, 28, "Text"], Cell[2501, 81, 105, 2, 47, "Input"], Cell[2609, 85, 55, 0, 28, "Text"], Cell[2667, 87, 155, 3, 47, "Input"], Cell[2825, 92, 221, 5, 58, "Text"], Cell[3049, 99, 88, 1, 47, "Input"], Cell[3140, 102, 52, 1, 47, "Input"], Cell[3195, 105, 20, 0, 28, "Text"], Cell[3218, 107, 76, 1, 47, "Input"], Cell[3297, 110, 43, 1, 47, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[3377, 116, 51, 0, 44, "Section"], Cell[3431, 118, 175, 5, 28, "Text"], Cell[3609, 125, 150, 2, 51, "Input"], Cell[3762, 129, 325, 6, 46, "Text"], Cell[4090, 137, 241, 5, 69, "Input"], Cell[4334, 144, 284, 9, 46, "Text"], Cell[4621, 155, 57, 1, 47, "Input"], Cell[4681, 158, 58, 0, 28, "Text"], Cell[4742, 160, 77, 1, 47, "Input"], Cell[4822, 163, 110, 5, 28, "Text"], Cell[4935, 170, 82, 1, 47, "Input"], Cell[5020, 173, 119, 3, 28, "Text"], Cell[5142, 178, 637, 11, 111, "Input"], Cell[5782, 191, 170, 4, 28, "Text"], Cell[5955, 197, 43, 1, 47, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[6035, 203, 72, 0, 44, "Section"], Cell[6110, 205, 183, 5, 28, "Text"], Cell[6296, 212, 163, 3, 47, "Input"], Cell[6462, 217, 111, 2, 47, "Input"], Cell[6576, 221, 155, 4, 28, "Text"], Cell[6734, 227, 128, 2, 51, "Input"], Cell[6865, 231, 147, 3, 28, "Text"], Cell[7015, 236, 275, 5, 46, "Text"], Cell[7293, 243, 211, 4, 67, "Input"], Cell[7507, 249, 76, 0, 28, "Text"], Cell[7586, 251, 72, 1, 47, "Input"], Cell[7661, 254, 202, 4, 46, "Text"], Cell[7866, 260, 90, 1, 47, "Input"], Cell[7959, 263, 84, 1, 47, "Input"], Cell[8046, 266, 100, 3, 28, "Text"], Cell[8149, 271, 52, 1, 47, "Input"], Cell[8204, 274, 68, 0, 28, "Text"], Cell[8275, 276, 171, 4, 51, "Input"], Cell[8449, 282, 171, 4, 51, "Input"], Cell[8623, 288, 45, 0, 28, "Text"], Cell[8671, 290, 128, 2, 47, "Input"], Cell[8802, 294, 94, 3, 28, "Text"], Cell[8899, 299, 75, 1, 47, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[9011, 305, 41, 0, 44, "Section"], Cell[9055, 307, 180, 4, 46, "Text"], Cell[9238, 313, 127, 2, 47, "Input"], Cell[9368, 317, 134, 3, 47, "Input"], Cell[9505, 322, 139, 3, 51, "Input"], Cell[9647, 327, 299, 7, 76, "Text"], Cell[9949, 336, 137, 2, 47, "Input"], Cell[10089, 340, 113, 4, 58, "Text"], Cell[10205, 346, 119, 2, 47, "Input"], Cell[10327, 350, 90, 3, 58, "Text"], Cell[10420, 355, 115, 2, 47, "Input"], Cell[10538, 359, 35, 0, 28, "Text"], Cell[10576, 361, 157, 3, 28, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[10770, 369, 45, 0, 44, "Section"], Cell[10818, 371, 246, 6, 46, "Text"], Cell[11067, 379, 201, 4, 46, "Text"], Cell[11271, 385, 305, 5, 70, "Input"], Cell[11579, 392, 54, 0, 28, "Text"], Cell[11636, 394, 221, 4, 63, "Input"], Cell[11860, 400, 44, 0, 28, "Text"], Cell[11907, 402, 50, 1, 47, "Input"], Cell[11960, 405, 73, 1, 47, "Input"], Cell[12036, 408, 178, 4, 28, "Text"], Cell[12217, 414, 370, 10, 46, "Text"], Cell[12590, 426, 1107, 24, 128, "Input", CellTags->"NDSolve"], Cell[13700, 452, 124, 3, 28, "Text"], Cell[13827, 457, 1560, 32, 193, "Input", CellTags->"NDSolve"], Cell[15390, 491, 183, 4, 46, "Text"], Cell[15576, 497, 43, 1, 47, "Input"], Cell[15622, 500, 254, 5, 63, "Input", CellTags->"NDSolve"], Cell[15879, 507, 100, 3, 28, "Text"], Cell[15982, 512, 124, 3, 47, "Input"], Cell[16109, 517, 105, 3, 28, "Text"], Cell[16217, 522, 110, 2, 47, "Input"], Cell[16330, 526, 747, 12, 130, "Text"], Cell[17080, 540, 230, 4, 64, "Input", CellTags->{"b:7.1.2", "ndsg:2.0.1.2"}], Cell[17313, 546, 62, 0, 28, "Text"], Cell[17378, 548, 70, 1, 47, "Input"], Cell[17451, 551, 140, 5, 28, "Text"], Cell[17594, 558, 1705, 34, 204, "Input", CellTags->"NDSolve"], Cell[19302, 594, 25, 0, 28, "Text"], Cell[19330, 596, 239, 5, 63, "Input", CellTags->"NDSolve"], Cell[19572, 603, 69, 1, 47, "Input"], Cell[19644, 606, 83, 1, 47, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[19764, 612, 33, 0, 46, "ExerciseMain"], Cell[19800, 614, 214, 8, 35, "Text"], Cell[CellGroupData[{ Cell[20039, 626, 33, 0, 48, "Exercise"], Cell[20075, 628, 247, 5, 46, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[20359, 638, 33, 0, 48, "Subsubsection"], Cell[20395, 640, 74, 1, 47, "Input", CellTags->"S3.5.11"], Cell[20472, 643, 70, 1, 47, "Input"], Cell[20545, 646, 70, 1, 47, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[20652, 652, 35, 0, 34, "Exercise"], Cell[20690, 654, 513, 18, 46, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[21240, 677, 33, 0, 48, "Subsubsection"], Cell[21276, 679, 45, 1, 47, "Input"], Cell[21324, 682, 158, 4, 51, "Input"], Cell[21485, 688, 290, 6, 63, "Input"], Cell[21778, 696, 324, 6, 75, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[22139, 707, 57, 0, 34, "Exercise"], Cell[22199, 709, 742, 19, 60, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[22978, 733, 33, 0, 48, "Subsubsection"], Cell[23014, 735, 119, 3, 28, "Text"], Cell[23136, 740, 93, 1, 66, "Input"], Cell[23232, 743, 277, 5, 63, "Input"], Cell[23512, 750, 83, 1, 47, "Input"], Cell[23598, 753, 46, 0, 28, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[23681, 758, 30, 0, 34, "Exercise"], Cell[23714, 760, 188, 5, 28, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[23939, 770, 33, 0, 48, "Subsubsection"], Cell[23975, 772, 166, 3, 47, "Input"], Cell[24144, 777, 77, 1, 47, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[24258, 783, 33, 0, 34, "Exercise"], Cell[24294, 785, 193, 5, 28, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[24524, 795, 33, 0, 48, "Subsubsection"], Cell[24560, 797, 241, 4, 79, "Input"], Cell[24804, 803, 70, 1, 47, "Input"], Cell[24877, 806, 91, 1, 47, "Input"], Cell[24971, 809, 189, 4, 47, "Input"], Cell[25163, 815, 77, 1, 47, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[25277, 821, 32, 0, 34, "Exercise"], Cell[25312, 823, 797, 17, 84, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[26146, 845, 33, 0, 48, "Subsubsection"], Cell[26182, 847, 120, 3, 28, "Text"], Cell[26305, 852, 89, 1, 64, "Input"], Cell[26397, 855, 121, 5, 28, "Text"], Cell[26521, 862, 67, 1, 47, "Input"], Cell[26591, 865, 30, 0, 28, "Text"], Cell[26624, 867, 1593, 32, 207, "Input", CellTags->"NDSolve"], Cell[28220, 901, 33, 0, 28, "Text"], Cell[28256, 903, 432, 9, 79, "Input", CellTags->"NDSolve"], Cell[28691, 914, 36, 0, 28, "Text"], Cell[28730, 916, 65, 1, 47, "Input"] }, Closed]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)