(************** Content-type: application/mathematica **************
CreatedBy='Mathematica 5.0'
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[ 67305, 2039]*)
(*NotebookOutlinePosition[ 67997, 2063]*)
(* CellTagsIndexPosition[ 67953, 2059]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[BoxData[
\(<< Graphics`\)], "Input",
InitializationCell->True],
Cell[CellGroupData[{
Cell[TextData[{
"Using ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" to Solve Differential Equations"
}], "Subtitle"],
Cell[CellGroupData[{
Cell["Single differential equations", "Section"],
Cell[CellGroupData[{
Cell["Symbolic solutions", "Subsection"],
Cell[TextData[{
"The basic command in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" for solving a differential equation symbolically is DSolve. This works \
for single differential equations to find a general solution, initial value \
problems, and systems of differential equations. The one thing you ",
StyleBox["must",
FontSlant->"Italic"],
" be careful to do when working with differential equations in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" is to explicitly specify each function and its independent variable. So, \
in order to enter the differential equation: ",
Cell[BoxData[
\(TraditionalForm\`y'' - 2\ y' + y = sin\ x\)]],
", you would have to enter it as: ",
Cell[BoxData[
\(TraditionalForm\`\(y''\)[x] - 2\ \(y'\)[x] + y[x] \[Equal]
Sin[x]\)]],
" (notice that I also converted it into \"",
StyleBox["Mathematica",
FontSlant->"Italic"],
"-speak\" at the same time; you must you == instead of =, square brackets, \
and capitalize all built-in functions)."
}], "Text"],
Cell[CellGroupData[{
Cell["General solutions", "Subsubsection"],
Cell[BoxData[
\(DSolve[\(y''\)[x] + 2\ \(y'\)[x] + y[x] \[Equal] Sin[x], y[x],
x]\)], "Input"],
Cell[TextData[{
"Notice that, in addition to specifying the differential equation to solve, \
you must also specify the function to solve for (",
Cell[BoxData[
\(TraditionalForm\`y[x]\)]],
") and the independent variable (x). In the solution, C[1] and C[2] are \
the two arbitrary constants for the general solution to this differential \
equation.\n\nIf all you need to do is see your answer, the above works fine. \
However, normally, you want to take this solution and do other things with it \
(plug numbers into it, graph it, make it into small table decorations, etc.). \
The form that ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" gives you the answer in here (",
Cell[BoxData[
\(TraditionalForm\`y[x] \[Rule] \ ... \)]],
") is what is called a \"replacement rule\". On the one hand, this is a \
very powerful form; you can use it in some very fancy ways. However, it is a \
little tricky to use if you aren't pretty familiar with ",
StyleBox["Mathematica",
FontSlant->"Italic"],
". So, most of the time we will want to simply assign this solution to a \
function we can work with later:"
}], "Text"],
Cell[BoxData[
\(Clear[ySol, diffEqn]\)], "Input"],
Cell["\<\
(You should get in the habit of Clear-ing all new names that you are about to \
use. Otherwise, old definitions you have forgotten about could really screw \
things up...)\
\>", "Text"],
Cell[BoxData[
\(diffEqn = \(y''\)[x] + 2\ \(y'\)[x] + y[x] \[Equal] Sin[x]\)], "Input"],
Cell[BoxData[
\(ySol[x_] =
y[x] /. First[
DSolve[\(y''\)[x] + 2\ \(y'\)[x] + y[x] \[Equal] Sin[x], y[x],
x]]\)], "Input"],
Cell[TextData[{
"This defines a function ",
Cell[BoxData[
\(TraditionalForm\`ySol[x]\)]],
" to be our solution function. (The /. is an operator that says \"replace \
all the y[x]'s to the left of this using the rule that follows\", which in \
this case is the solution to the differential equation, which said ",
Cell[BoxData[
\(y[
x] \[Rule] \[ExponentialE]\^\(-x\)\ C[
1] + \[ExponentialE]\^\(-x\)\ x\ C[2] - Cos[x]\/2\)]],
". The First command strips out the extra {}'s in the solution.) \n\nNow, \
we check our answer:"
}], "Text"],
Cell[BoxData[
\(diffEqn /. y \[Rule] ySol // FullSimplify\)], "Input"],
Cell["\<\
The above line says \"Take the original diffEqn, replace all the y's in it \
with ySol (our solution) and then simplify the result.\" If the answer is \
\"True\" then the solution satisfies the differential equation.\
\>", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["\<\
*Useful hint: If you are getting complex solutions (which you presumably \
don't want), read this section.\
\>", "Subsubsection"],
Cell["\<\
Try adding one of the following to the end of the DSolve line:
//ComplexExpand//Simplify
OR
//ComplexExpand//FullSimplify
Unfortunately, even this doesn't always work.At that point,you have 4 options \
(that I can think of, anyway):
1) Solve the problem numerically using NDSolve (probably the best idea, if \
you can get enough accuracy).
2) Plot the graph anyway (it will just have some blank spots in it)
3) Assume the imaginary part is ignorable and take the Real part of your \
answer (using the Re function). To be honest,since I'm not entirely sure \
howMathematica's algorithms work for diff equations, I would be a little \
hesitant to do this without knowing a bit more. (Notice that this is really a \
fancier version of option 2.)
4) Use some other method to solve the problem symbolically (i.e., solve the \
characteristic polynomial numerically, then use what you know about the form \
of the homogeneous solution and the particular solution). This is probably \
the best way, but it can be a lot of extra work (and probably won't get you \
much extra payoff). I would save this for when you are going to do some \
serious work with that particular diff. equation.\
\>", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Initial value problems", "Subsubsection"],
Cell[TextData[{
"Most of the time when you work with a differential equation on a computer \
system, you will be interested in \"initial value problems\", rather than \
just finding the general solution. This is handled very simply: you just \
include the initial conditions in a list with your differential equation to \
solve. So, to solve ",
Cell[BoxData[
\(TraditionalForm\`y'' - 2\ y' + y = sin\ x\)]],
", ",
Cell[BoxData[
\(TraditionalForm\`y(0) = 1, \ y' \((0)\) = \(-2\)\)]],
", you would:"
}], "Text"],
Cell[BoxData[
\(Clear[x, y, ySol, diffEqn]\)], "Input"],
Cell[BoxData[
\(diffEqn = \(y''\)[x] + 2\ \(y'\)[x] + y[x] \[Equal] Sin[x]\)], "Input"],
Cell[BoxData[
\(ySol[x_] =
y[x] /. First[
DSolve[{\(y''\)[x] + 2\ \(y'\)[x] + y[x] \[Equal] Sin[x],
y[0] \[Equal] 1, \(y'\)[0] \[Equal] \(-2\)}, y[x],
x]]\)], "Input"],
Cell["\<\
(Again, be sure to use == instead of = everywhere inside DSolve.)
We could graph this solution curve:\
\>", "Text"],
Cell[BoxData[
\(Plot[ySol[x], {x, 0, 4 \[Pi]}]\)], "Input"],
Cell["\<\
What if we wanted to graph several solution curves, for different initial \
conditions? Let's examine the impact of curves with different initial y \
values. To allow ourselves the most flexibility, let's use some dummy \
variables:\
\>", "Text"],
Cell[BoxData[
\(Clear[y0, v0]\)], "Input"],
Cell[BoxData[
\(ySol[x_] =
y[x] /. First[
DSolve[{\(y''\)[x] + 2\ \(y'\)[x] + y[x] \[Equal] Sin[x],
y[0] \[Equal] y0, \(y'\)[0] \[Equal] v0}, y[x], x]]\)], "Input"],
Cell[BoxData[
\(Plot[
ySol[x] /. {y0 \[Rule] \(-1\), v0 \[Rule] 2}, {x, 0,
4 \[Pi]}]\)], "Input"],
Cell[TextData[{
"This says to plot ySol again, but replace all the y0's with -1 and v0's \
with 2. A more interesting question would be to hold the v0's constant (say \
at ",
Cell[BoxData[
\(TraditionalForm\`v0 = 1\)]],
") and let the y0's vary from -5 to 5. There is a really cool way to do \
this:"
}], "Text"],
Cell[BoxData[
\(Clear[v0, y0]\)], "Input"],
Cell[BoxData[
\(v0 = 1; y0 = Range[\(-5\), 5, 1]\)], "Input"],
Cell[BoxData[
\(Plot[Evaluate[ySol[x]], {x, 0, 4 \[Pi]},
PlotRange \[Rule] All]\)], "Input"],
Cell[TextData[{
"This plots ySol for all the different values we specified for y0. (You ",
StyleBox["must",
FontSlant->"Italic"],
" use the Evaluate command around ySol[x] or ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" won't plug in the numbers for y0 before trying to graph and it won't \
work.) So, what happens if we let y0 be constant (say ",
Cell[BoxData[
\(TraditionalForm\`y0 = 2\)]],
") and let v0 vary from -5 to 5?"
}], "Text"],
Cell[BoxData[
\(Clear[v0, y0]\)], "Input"],
Cell[BoxData[
\(y0 = 2; v0 = Range[\(-5\), 5, 1]\)], "Input"],
Cell[BoxData[
\(Plot[Evaluate[ySol[x]], {x, 0, 4 \[Pi]},
PlotRange \[Rule] All]\)], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["Unknown parameters", "Subsubsection"],
Cell[TextData[{
"This trick is also very handy if you want to examine the effect of \
different parameters on a differential equation. So, let's consider a new \
equation (oh boy...): ",
Cell[BoxData[
\(TraditionalForm\`x'' + k\ x = sin\ t\)]]
}], "Text"],
Cell[BoxData[
\(Clear[diffEqn, x0, v0, k, xSol, x, t]\)], "Input"],
Cell[BoxData[
\(diffEqn = \(x''\)[t] + k\ x[t] \[Equal] Sin[t]\)], "Input"],
Cell[BoxData[
\(xSol[t_] =
FullSimplify[
x[t] /. First[
DSolve[{diffEqn, x[0] == x0, \(x'\)[0] \[Equal] v0}, x[t], t]],
k \[Element] Reals]\)], "Input"],
Cell[TextData[{
"(The ",
Cell[BoxData[
\(TraditionalForm\`FullSimplify[\(...\)\(,\)\(k \[Element]
Reals\)]\)]],
" is a useful trick to let ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" know that your k is a real number and simplify accordingly. It doesn't \
actually matter in this particular problem, but sometimes it does...)"
}], "Text"],
Cell[BoxData[
\(x0 = 1; v0 = 1; k = Range[0.1, 3.1, 0.5]\)], "Input"],
Cell[BoxData[
\(Plot[Evaluate[xSol[t]], {t, 0, 4 \[Pi]}]\)], "Input"],
Cell[TextData[{
"Notice that I avoid ",
Cell[BoxData[
\(TraditionalForm\`k = 0\)]],
". So, what if k is negative?"
}], "Text"],
Cell[BoxData[
\(x0 = 1; v0 = 1; k = Range[\(-3.1\), \(-0.1\), 0.5]\)], "Input"],
Cell[BoxData[
\(Plot[Evaluate[xSol[t]], {t, 0, 4 \[Pi]}]\)], "Input"]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell["Numeric approxiamtions", "Subsection"],
Cell[TextData[{
"Sometimes, ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" can't solve a differential equation (if it is evil enough). Or, it might \
take a very long time for it to solve and you might not really have any need \
for a complete symbolic solution. Often, a good numerical approximation is \
all you really need. In ",
StyleBox["Mathematica",
FontSlant->"Italic"],
", you you NDSolve (Numeric DSolve) to compute these. One of the nice \
things about the way ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" does this is that it doesn't just give you a numeric resulet (i.e., ",
Cell[BoxData[
\(TraditionalForm\`x[2] = 3.2222152\)]],
"), it gives you a numerical \"interpolating function\" which you can then \
use to compute the value of the function at any point in your computed \
domain. Thus, you can graph this solution, find its intercepts, etc. Not \
surprisingly, you ",
StyleBox["must",
FontSlant->"Italic"],
" provide enough initial conditions to make sure that the solution is \
unique (i.e., no unknown constants).\n\nLet's solve the following \
differential equation: ",
Cell[BoxData[
\(TraditionalForm\`y\^\((3)\) + y'' + y' = \(-y\^3\)\)]],
" (remember that ",
Cell[BoxData[
\(TraditionalForm\`y\^\((3)\)\)]],
" means the 3rd derivative), where ",
Cell[BoxData[
\(TraditionalForm\`y(0) = 1, \ y' \((0)\) = 0, \ y'' \((0)\) = 0\)]],
". We must also specify a domain over which we want the solution to be \
valid, say ",
Cell[BoxData[
\(TraditionalForm\`0 \[LessEqual] x \[LessEqual] 20\)]],
"."
}], "Text"],
Cell[BoxData[
\(Clear[x, y, ySol]\)], "Input"],
Cell[BoxData[
RowBox[{\(ySol[x_]\), "=",
RowBox[{\(y[x]\), "/.",
RowBox[{"First", "[",
RowBox[{"NDSolve", "[",
RowBox[{
RowBox[{"{",
RowBox[{
RowBox[{
RowBox[{
RowBox[{
SuperscriptBox["y",
TagBox[\((3)\),
Derivative],
MultilineFunction->None], "[", "x", "]"}], "+",
RowBox[{
SuperscriptBox["y", "\[DoublePrime]",
MultilineFunction->None], "[", "x", "]"}], "+",
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "[", "x", "]"}]}],
"==", \(-y[x]\^3\)}], ",", \(y[0] \[Equal] 1\), ",",
RowBox[{
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "[", "0", "]"}], "\[Equal]",
RowBox[{
SuperscriptBox["y", "\[DoublePrime]",
MultilineFunction->None], "[", "0", "]"}], "\[Equal]",
"0"}]}], "}"}], ",", \(y[x]\), ",", \({x, 0, 20}\)}],
"]"}], "]"}]}]}]], "Input"],
Cell[BoxData[
\(ySol[20]\)], "Input"],
Cell[BoxData[
\(Plot[ySol[x], {x, 0, 20}]\)], "Input"],
Cell["\<\
If you need greater accuracy in your approximation, you can set the \
WorkingPrecision to be higher (actually, if you really want to investigate \
this in more detail, you should read up on WorkingPrecision, PrecisionGoal, \
and AccuracyGoal in the online help; it can get a bit complicated - setting \
WorkingPrecision to a large number will often help, however). The following \
does all internal computations to an internal accuracy of 30 digits (so your \
answer is probably accurate to about 20 digits):\
\>", "Text"],
Cell[BoxData[
RowBox[{\(ySol[x_]\), "=",
RowBox[{\(y[x]\), "/.",
RowBox[{"First", "[",
RowBox[{"NDSolve", "[",
RowBox[{
RowBox[{"{",
RowBox[{
RowBox[{
RowBox[{
RowBox[{
SuperscriptBox["y",
TagBox[\((3)\),
Derivative],
MultilineFunction->None], "[", "x", "]"}], "+",
RowBox[{
SuperscriptBox["y", "\[DoublePrime]",
MultilineFunction->None], "[", "x", "]"}], "+",
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "[", "x", "]"}]}],
"==", \(-y[x]\^3\)}], ",", \(y[0] \[Equal] 1\), ",",
RowBox[{
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "[", "0", "]"}], "\[Equal]",
RowBox[{
SuperscriptBox["y", "\[DoublePrime]",
MultilineFunction->None], "[", "0", "]"}], "\[Equal]",
"0"}]}], "}"}], ",", \(y[x]\), ",", \({x, 0, 20}\),
",", \(WorkingPrecision \[Rule] 30\)}], "]"}],
"]"}]}]}]], "Input"],
Cell["\<\
It's a bit trickier here if you want to graph one of these for a variety of \
initial conditions:\
\>", "Text"],
Cell[BoxData[
RowBox[{"ySolutionList", "=",
RowBox[{"Table", "[",
RowBox[{
RowBox[{\(y[x]\), "/.",
RowBox[{"First", "[",
RowBox[{"NDSolve", "[",
RowBox[{
RowBox[{"{",
RowBox[{
RowBox[{
RowBox[{
RowBox[{
SuperscriptBox["y",
TagBox[\((3)\),
Derivative],
MultilineFunction->None], "[", "x", "]"}], "+",
RowBox[{
SuperscriptBox["y", "\[DoublePrime]",
MultilineFunction->None], "[", "x", "]"}], "+",
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "[", "x", "]"}]}],
"==", \(-y[x]\^3\)}], ",", \(y[0] \[Equal] y0\), ",",
RowBox[{
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "[", "0", "]"}],
"\[Equal]",
RowBox[{
SuperscriptBox["y", "\[DoublePrime]",
MultilineFunction->None], "[", "0", "]"}],
"\[Equal]", "0"}]}], "}"}], ",", \(y[x]\),
",", \({x, 0, 20}\)}], "]"}], "]"}]}],
",", \({y0, 0, 1.25, 0.25}\)}], "]"}]}]], "Input"],
Cell[BoxData[
\(Plot[Evaluate[ySolutionList], {x, 0, 20}]\)], "Input"],
Cell["\<\
If you want to solve for multiple initial conditions (like above) or for \
different values of a parameter, the important thing to remember is that \
NDSolve will not work with unknown constants in it (unlike DSolve), so you \
have to plug those constants in multiple times (using the Table command like \
I did above is a good way to do this). (You could actually approach DSolve \
the same way if you wish, but if the differential equation takes a long time \
to solve, it could really add up to a long solution time in that case).\
\>", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Plotting solutions and slope fields", "Subsection"],
Cell[CellGroupData[{
Cell["Slope fields", "Subsubsection"],
Cell[TextData[{
"To graph a slope-field in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
", you use the command PlotVectorField. The \"vector field\" you are \
actually plotting has the coordinates ",
Cell[BoxData[
\(TraditionalForm\`{1, y'}\)]],
" at each point (since each vector goes \"over\" 1 and \"up\" by the \
slope). So, to graph the slope field for: ",
Cell[BoxData[
\(TraditionalForm\`y' = x\/\(2 + sin\ y\)\)]]
}], "Text"],
Cell[BoxData[
\(<< Graphics`\)], "Input"],
Cell[BoxData[
\(PlotVectorField[{1, x\/\(2 + Sin[y]\)}, {x, \(-10\), 10}, {y, \(-10\),
10}]\)], "Input"],
Cell["\<\
To clean it up a bit and make it look nicer, try the following options:\
\>", "Text"],
Cell[BoxData[
\(PlotVectorField[{1, x\/\(2 + Sin[y]\)}, {x, \(-10\), 10}, {y, \(-10\),
10}, Axes \[Rule] True]\)], "Input"],
Cell["\<\
If you prefer to make all of the vectors have the same length, try the \
following:\
\>", "Text"],
Cell[BoxData[
\(slopeField =
PlotVectorField[{1, x\/\(2 + Sin[y]\)}, {x, \(-10\), 10}, {y, \(-10\),
10}, Axes \[Rule] True, ScaleFunction \[Rule] \((0.4 &)\),
ScaleFactor \[Rule] 0.9]\)], "Input"],
Cell[TextData[{
"It's a bit complicated how the ",
Cell[BoxData[
\(TraditionalForm\`ScaleFunction \[Rule] \((0.4 &)\),
ScaleFactor \[Rule] 0.9\)]],
" options work, but if you play around with those numbers a bit you can get \
something that looks quite nice (be ",
StyleBox["sure",
FontSlant->"Italic"],
" that you copy ScaleFunction exactly and then play with the number only; \
the & is very important)."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Putting the solutions and the slope fields together", "Subsubsection"],
Cell[TextData[{
"Often, you will want to plot the solution curves for different initial \
conditions along with the slope field together on one graph. You can do that \
by combining the solution graphing we did above with the slope fields.\n\n",
"Because I picked a pretty nasty equation, DSolve doesn't work like I want, \
so we will do this one using NDSolve. Let's solve for the different initial \
values for the curves we want to plot (say from ",
Cell[BoxData[
\(TraditionalForm\`y(0) = \(-8\)\)]],
" to ",
Cell[BoxData[
\(TraditionalForm\`y(0) = 8\)]],
", every 2 units):"
}], "Text"],
Cell[BoxData[
\(Clear[x, y, y0, ySolutions, plotSolutions]\)], "Input"],
Cell[BoxData[
\(ySolutions =
Table[y[x] /.
NDSolve[{\(y'\)[x] \[Equal] x\/\(2 + Sin[y[x]]\),
y[0] \[Equal] y0}, y[x], {x, \(-10\), 10}], {y0, \(-8\), 8,
2}]\)], "Input"],
Cell[BoxData[
\(plotSolutions =
Plot[Evaluate[ySolutions], {x, \(-10\), 10},
PlotStyle \[Rule] Blue]\)], "Input"],
Cell["Now, let's combine the slope field with the solution curves:", "Text"],
Cell[BoxData[
\(Show[slopeField, plotSolutions]\)], "Input"],
Cell["\<\
This is, of course, hideous. To fix that, we can explicitly tell it what \
region to plot: \
\>", "Text"],
Cell[BoxData[
\(Show[slopeField, plotSolutions,
PlotRange \[Rule] {{\(-10\), 10}, {\(-10\), 10}}]\)], "Input"]
}, Closed]]
}, Closed]]
}, Open ]],
Cell[CellGroupData[{
Cell["Systems of differential equations", "Section"],
Cell[TextData[{
"You will notice that I use several different methods to solve systems of \
equations here. This is partly because some of them (the sections on \
eigensystems and exponentiation of operators) are techniques covered in our \
text; I show you how to work those methods in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
". Honestly, however, most of the time you will probably use DSolve and \
NDSolve. Be ",
StyleBox["sure",
FontSlant->"Italic"],
" to read the first section on setting up systems in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
", however."
}], "Text"],
Cell[CellGroupData[{
Cell["Setting up the systems - working with matrices", "Subsection"],
Cell[TextData[{
"There are a lot of different ways to handle systems of differential \
equations in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
". The most simplistic method is to just enter them as lists of \
differential equations:"
}], "Text"],
Cell[BoxData[
\(DSolve[{\(x'\)[t] \[Equal] x[t] - 2 y[t], \(y'\)[t] \[Equal]
2 x[t] + y[t]}, {x[t], y[t]}, t]\)], "Input"],
Cell[TextData[{
"This method works (and is easy to understand), but starts to get kind of \
cumbersome for larger systems (and certainly doesn't win any cool points...). \
Personally, I prefer to set systems up from the outset as matrix equations.\n\
\nTo do this, you need to know a few things about how ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" handles matrices and vectors. As far as ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" is concerned, a vector is just a ",
StyleBox["list",
FontSlant->"Italic"],
" of items: ",
Cell[BoxData[
\(TraditionalForm\`{2, x, y}\)]],
" A matrix is just a list of lists: ",
Cell[BoxData[
\(TraditionalForm\`{{1, 2, 3}, {3, 5, 2}, {7, 2, 0}}\)]],
", where each inside list is a row of the matrix. So this last matrix would \
look like:"
}], "Text"],
Cell[BoxData[
\({{1, 2, 3}, {3, 5, 2}, {7, 2, 0}} // MatrixForm\)], "Input"],
Cell[TextData[{
"In fact, you can enter matrices directly (so they actually look like \
matrices) by using the \"Create Table/Matrix/Palette\" command under the \
\"Input\" menu (then just choose \"Matrix\" and the number of rows/columns) \
or, if you need a 2 by 2 matrix, there is a button on the standard tool \
palette for this.\n\nTo multiply two matrices, you ",
StyleBox["must",
FontSlant->"Italic"],
" use . between them; if you use implied multiplation (or the standard \
\[Times] symbol), ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" won't use matrix multiplication. So:"
}], "Text"],
Cell[BoxData[
RowBox[{
RowBox[{
RowBox[{"(", GridBox[{
{"2", "4"},
{"6", "8"}
}], ")"}], ".",
RowBox[{"(", GridBox[{
{\(-3\), "4"},
{\(-5\), "2"}
}], ")"}]}], "//", "MatrixForm"}]], "Input"],
Cell["Or:", "Text"],
Cell[BoxData[
RowBox[{
RowBox[{"A", "=",
RowBox[{"(", GridBox[{
{"2", "4"},
{"6", "8"}
}], ")"}]}], ";",
RowBox[{"B", "=",
RowBox[{"(", GridBox[{
{\(-3\), "4"},
{\(-5\), "2"}
}], ")"}]}]}]], "Input"],
Cell[BoxData[
\(A . B // MatrixForm\)], "Input"],
Cell[TextData[{
"(Sticking //MatrixForm at the end of the line just makes the answer look \
nice. It isn't necessary.)\n\nOne other thing ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" does that is kind of nice is that it handles multiplication of a matrix \
times a vector properly (if you use .), even if the vector isn't really the \
right kind of vector (i.e., you can be sloppy about row vectors vs. column \
vectors). We will take advantage of this to simplify some work that ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" would otherwise force you to do. (I will treat all vectors as row \
vectors, even when it would be more correct to use column vectors. Trust me, \
this simplifies things a lot...)\n\nIn order for vector equations to work \
correctly in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
", you must first execute the following command (this is ",
StyleBox["important",
FontSlant->"Italic"],
", so pay attention):"
}], "Text"],
Cell[BoxData[
\(SetAttributes[Equal, Listable]\)], "Input",
InitializationCell->True],
Cell[TextData[{
"(If you are interested, this makes it so the == sign is \"listable\", \
i.e., it works with matrices and other lists.)\n\nNow, let's pull all this \
together and start using matrix notation. First, let's just set up a system \
of algebraic equations. The system: \n",
Cell[BoxData[
FormBox[GridBox[{
{\(2 x - 3 y + 4 z = 7\)},
{\(2 x - 5 z = 4\)},
{\(7 y + 4 x = 9\)}
}], TraditionalForm]]],
" can be written in matrix form with the following definitions:"
}], "Text"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"2", \(-3\), "4"},
{"2", "0", \(-5\)},
{"0", "7", "4"}
}], ")"}]}], ";", \(vC = {7, 4, 9}\),
";", \(vX = {x, y, z}\)}]], "Input"],
Cell["\<\
One note about notation: For clarity sake (and to avoid using variable names \
that begin with capital letters), I will name all matrices beginning with an \
m (so mA, mC, mWhatever) and all vectors beginning with a v (vX, vSolution, \
vSomething).
So the system becomes:\
\>", "Text"],
Cell[BoxData[
\(system = mA . vX \[Equal] vC\)], "Input"],
Cell[BoxData[
\(Solve[system]\)], "Input"],
Cell[TextData[{
"So we can solve the system very easily this way. \n\nNow, let's set up a \
system of differential equations. Let's set up:\n",
Cell[BoxData[GridBox[{
{\(x' = x - 2 y\)},
{\(y' = 2 x + y\)}
}]]],
" with initial conditions ",
Cell[BoxData[
\(TraditionalForm\`x(0) = 6, \ y(0) = 4\)]],
".\t"
}], "Text"],
Cell[BoxData[
\(Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"1", \(-2\)},
{"2", "1"}
}], ")"}]}], ";", \(vX[t_] = {x[t], y[t]}\)}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t]\)], "Input"],
Cell[BoxData[
\(vSolution[t_] = vX[t] /. DSolve[system, vX[t], t]\)], "Input"],
Cell["\<\
Or, if you prefer to include the initial conditions in the solution:\
\>", "Text"],
Cell[BoxData[
\(mX0 = {x0, y0}\)], "Input"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. DSolve[{system, vX[0] \[Equal] mX0}, vX[t], t]\)], "Input"],
Cell["A more complicated system might look something like this:", "Text"],
Cell[BoxData[
\(Clear[mA, vX, vC, x, y, t, x0, y0, vX0, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"1", \(-2\)},
{"2", "1"}
}], ")"}]}], ";", \(vC = {2, t}\),
";", \(vX[t_] = {x[t], y[t]}\)}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t] + vC\)], "Input"],
Cell[BoxData[
\(vSolution[t_] = vX[t] /. DSolve[system, vX[t], t]\)], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["Working with matrices and eigensystems", "Subsection"],
Cell[TextData[{
"Sometimes, the solution that DSolve gives to a system isn't really what \
you were looking for (it may look much more complicated than you expected). \
This is because ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" has to be able to solve some very complicated systems. Sometimes, if \
your system is relatively simple (an autonomous homogeneous system, for \
example), you might be able to get a simpler solution by finding the \
eigenvalues and eigenvectors. For example, if you are trying to solve the \
following system:"
}], "Text"],
Cell[BoxData[
\(Clear[mA, vX, x1, x2, x3, x4, system, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"2", "3", \(-2\), "1"},
{"5", "0", "6", \(-1\)},
{"7", "5", "9", \(-7\)},
{"8", "4", "2", "1"}
}], ")"}]}],
";", \(vX[t_] = {x1[t], x2[t], x3[t], x4[t]}\)}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t]\)], "Input"],
Cell["\<\
Let's find the eigenvalues and eigenvectors of this system (I stuck the \
N[...] in there to get a numeric result because the exact values are pretty \
insane...):\
\>", "Text"],
Cell[BoxData[
\(eigValues = Eigenvalues[N[mA]]\)], "Input"],
Cell[BoxData[
\(eigVectors = Eigenvectors[N[mA]]\)], "Input"],
Cell["\<\
Notice that, since all the coefficients of the original system are real, the \
complex ones occur in complex conjugate pairs. Now, you can just cut and \
paste the appropriate eigenvalues/eigenvectors to construct the solution, or \
with a little more finesse (I'll take them from last to first; \
eigVectors[[4]] means the 4th vector in the list):\
\>", "Text"],
Cell[BoxData[
\(vX4 =
eigVectors[\([4]\)]\ \[ExponentialE]\^\(eigValues[\([4]\)]\ t\)\)], \
"Input"],
Cell[BoxData[
\(vX3 =
eigVectors[\([3]\)]\ \[ExponentialE]\^\(eigValues[\([3]\)]\ t\)\)], \
"Input"],
Cell["\<\
Now, the other two eigenvalues are complex, so recall that we handle those a \
bit differently:\
\>", "Text"],
Cell[BoxData[
\(vX2 = \(\[ExponentialE]\^\(Re[eigValues[\([1]\)]]
t\)\) \((Re[eigVectors[\([1]\)]]
Cos[Im[eigValues[\([1]\)]] t] -
Im[eigVectors[\([1]\)]]
Sin[Im[eigValues[\([1]\)]] t])\)\)], "Input"],
Cell[BoxData[
\(vX1 = \(\[ExponentialE]\^\(Re[eigValues[\([1]\)]]
t\)\) \((Re[eigVectors[\([1]\)]]
Cos[Im[eigValues[\([1]\)]] t] +
Im[eigVectors[\([1]\)]]
Sin[Im[eigValues[\([1]\)]] t])\)\)], "Input"],
Cell["\<\
(Wasn't that fun? Make sure you understand where all of these came from \
before you continue.)
Now, we can find the general solution by combining these:\
\>", "Text"],
Cell[BoxData[
\(vSolution[t_] =
C[1] vX1 + C[2] vX2 + C[3] vX3 + C[4] vX4 // FullSimplify\)], "Input"],
Cell[BoxData[
\(vSolution[t] // MatrixForm\)], "Input"],
Cell["If the initial conditions are given by:", "Text"],
Cell[BoxData[
\(vX0 = {1, 1, 1, 1}\)], "Input"],
Cell["We can solve for the constants:", "Text"],
Cell[BoxData[
\(initConditions = Solve[vSolution[0] \[Equal] vX0]\)], "Input"],
Cell["\<\
The error here is a warning that we probably need to increase our numeric \
precision for all of our calculations (especially finding the eigenvalues and \
eigenvectors); round-off errors could easily screw things up for this \
complicated a situation. However, we finally get:\
\>", "Text"],
Cell[BoxData[
\(\(vSolution[t] /. initConditions // Transpose\) //
MatrixForm\)], "Input"],
Cell[TextData[{
"Of course, if you have eigenvalues with multiplicity greater than 1, then \
you have a lot more pain and suffering to deal with. Unfortunately, ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" doesn't automatically give you all the eigenvalues (\"generalized \
eigenvalues\") you need for this, though you can certainly use ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" to find them. However, there is a much better way to deal with systems \
of equations (with constant coefficients, anyway), which I cover in the next \
section, so let's just skip that, shall we?"
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["\<\
Solving using exponentiation of matrices (or \"exponentiation of operators\")\
\
\>", "Subsection"],
Cell[TextData[{
"When you solve the single differential equation ",
Cell[BoxData[
\(TraditionalForm\`x' = a\ x, \ x(0) = x\_0\)]],
", you get the solution ",
Cell[BoxData[
\(TraditionalForm\`x(t) = \(x\_0\) \[ExponentialE]\^\(a\ t\)\)]],
". It turns out that if you are solving the system of differential \
equations ",
Cell[BoxData[
FormBox[
RowBox[{
RowBox[{
StyleBox[\(X'\),
FontWeight->"Bold"],
StyleBox["=",
FontWeight->"Plain"],
StyleBox[\(A . \ X\),
FontWeight->"Bold"]}],
StyleBox[",",
FontWeight->"Bold"],
StyleBox[" ",
FontWeight->"Bold"],
RowBox[{
RowBox[{
StyleBox["X",
FontWeight->"Bold"],
StyleBox["(",
FontWeight->"Plain"],
StyleBox["0",
FontWeight->"Plain"],
StyleBox[")",
FontWeight->"Plain"]}],
StyleBox["=",
FontWeight->"Plain"],
SubscriptBox[
StyleBox["X",
FontWeight->"Bold"], "0"]}]}], TraditionalForm]]],
" (where ",
StyleBox["A",
FontWeight->"Bold"],
" is a matrix of constant coefficients)",
", the solution winds up being ",
Cell[BoxData[
FormBox[
RowBox[{
RowBox[{
StyleBox["X",
FontWeight->"Bold"],
StyleBox["(",
FontWeight->"Plain"],
StyleBox["t",
FontWeight->"Plain"],
StyleBox[")",
FontWeight->"Plain"]}],
StyleBox["=",
FontWeight->"Plain"],
RowBox[{
SuperscriptBox[
StyleBox["\[ExponentialE]",
FontWeight->"Plain"],
RowBox[{
StyleBox["A",
FontWeight->"Bold"],
StyleBox[" ",
FontWeight->"Plain"],
StyleBox["t",
FontWeight->"Plain"]}]], ".",
SubscriptBox[
StyleBox["X",
FontWeight->"Bold"], "0"]}]}], TraditionalForm]]],
" (needless to say, there is a ",
StyleBox["lot",
FontSlant->"Italic"],
" of math that goes on behind the scenes to make this work out). Of \
course, for this to make any sense, you have to have some idea what it means \
to raise \[ExponentialE] to a matrix power. I won't go into all the details \
here, but the short form is that if you consider one definition of ",
Cell[BoxData[
\(TraditionalForm\`\[ExponentialE]\^a\)]],
" to be ",
Cell[BoxData[
\(TraditionalForm\`\[Sum]\+\(n = 0\)\%\[Infinity] a\^n\/\(n!\)\)]],
", then you can define ",
Cell[BoxData[
FormBox[
SuperscriptBox["\[ExponentialE]",
StyleBox["A",
FontWeight->"Bold"]], TraditionalForm]]],
" as ",
Cell[BoxData[
FormBox[
RowBox[{\(\[Sum]\+\(n = 0\)\%\[Infinity]\),
FractionBox[
SuperscriptBox[
StyleBox["A",
FontWeight->"Bold"], "n"], \(n!\)]}], TraditionalForm]]],
"(where ",
Cell[BoxData[
FormBox[
RowBox[{
SuperscriptBox[
StyleBox["A",
FontWeight->"Bold"], "0"], "=",
StyleBox["I",
FontWeight->"Bold"]}], TraditionalForm]]],
", the identity matrix), assuming the infinite sum actually converges. So, \
for example, ",
Cell[BoxData[
FormBox[
RowBox[{
SuperscriptBox["\[ExponentialE]",
RowBox[{
RowBox[{"(", GridBox[{
{"2", "0"},
{"0", "3"}
}], ")"}], "t"}]], "=",
RowBox[{
RowBox[{\(\[Sum]\+\(n = 0\)\%\[Infinity]\),
FractionBox[
RowBox[{
SuperscriptBox[
RowBox[{"(", GridBox[{
{"2", "0"},
{"0", "3"}
}], ")"}], "n"], "t"}], \(n!\)]}], "=",
RowBox[{
RowBox[{"(", GridBox[{
{\(\[Sum]\+\(n = 0\)\%\[Infinity]\((2 t)\)\^n\/\(n!\)\),
"0"},
{
"0", \(\[Sum]\+\(n = 0\)\%\[Infinity]\((3 \
t)\)\^n\/\(n!\)\)}
}], ")"}], "=",
RowBox[{"(", GridBox[{
{\(\[ExponentialE]\^\(2 t\)\), "0"},
{"0", \(\[ExponentialE]\^\(3 t\)\)}
}], ")"}]}]}]}], TraditionalForm]]],
" (this is, of course, a very simple example). The good news is that ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" has a function built in for dealing with this: MatrixExp. So, we could \
use:"
}], "Text"],
Cell[BoxData[
RowBox[{
RowBox[{"MatrixExp", "[",
RowBox[{
RowBox[{"(", GridBox[{
{"2", "0"},
{"0", "3"}
}], ")"}], "t"}], "]"}], "//", "MatrixForm"}]], "Input"],
Cell["\<\
So, let's work the problem from the previous section again, shall we?\
\>", "Text"],
Cell[BoxData[
\(Clear[mA, vX0, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"2", "3", \(-2\), "1"},
{"5", "0", "6", \(-1\)},
{"7", "5", "9", \(-7\)},
{"8", "4", "2", "1"}
}], ")"}]}], ";", \(vX0 = {1, 1, 1, 1}\), ";"}]], "Input"],
Cell[BoxData[
\(vSolution[t_] = MatrixExp[N[mA]\ t] . vX0 // ComplexExpand\)], "Input"],
Cell["\<\
(Notice that I use N[mA] in order to get a nicely readable solution, sort \
of...)\
\>", "Text"],
Cell[BoxData[
\(vSolution[t] // MatrixForm\)], "Input"],
Cell[TextData[{
"Notice, however, the imaginary components in the answer. This is due to \
some round-off errors. Now, for positive ",
StyleBox["t",
FontSlant->"Italic"],
" values, this turns out to be relatively insignificant. Compare:"
}], "Text"],
Cell[BoxData[
\(vSolution[100. ] // MatrixForm\)], "Input"],
Cell[TextData[{
"However, for negative ",
StyleBox["t",
FontSlant->"Italic"],
" values, it certainly isn't:"
}], "Text"],
Cell[BoxData[
\(vSolution[\(-100. \)] // MatrixForm\)], "Input"],
Cell["\<\
You can fix this in one of two ways: The best way is to keep things symbolic \
as long as possible:\
\>", "Text"],
Cell[BoxData[
\(vSolution[t_] = MatrixExp[mA\ t] . vX0 // ComplexExpand\)], "Input"],
Cell["This isn't a very readable solution, but notice that:", "Text"],
Cell[BoxData[
\(vSolution[100. ] // MatrixForm\)], "Input"],
Cell[BoxData[
\(vSolution[\(-100. \)] // MatrixForm\)], "Input"],
Cell["\<\
The other method is to just take the real components of your solution (and \
hope that nothing important is going on behind the scenes with that imaginary \
component):\
\>", "Text"],
Cell[BoxData[
\(vSolution[t_] = \(MatrixExp[N[mA]\ t] . vX0 // Re\) //
ComplexExpand\)], "Input"],
Cell[BoxData[
\(vSolution[100. ] // MatrixForm\)], "Input"],
Cell[BoxData[
\(vSolution[\(-100. \)] // MatrixForm\)], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["General symbolic solutions", "Subsection"],
Cell["\<\
Of course, if you really don't want to fool with any of the above, you could \
just call up DSolve:\
\>", "Text"],
Cell[BoxData[
\(Clear[mA, vX0, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"2", "3", \(-2\), "1"},
{"5", "0", "6", \(-1\)},
{"7", "5", "9", \(-7\)},
{"8", "4", "2", "1"}
}], ")"}]}], ";", \(vX0 = {1, 1, 1, 1}\),
";", \(vX[t_] = {x1[t], x2[t], x3[t], x4[t]}\), ";"}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t]\)], "Input"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. First[
DSolve[{system, vX[0] \[Equal] vX0}, vX[t], t]]\)], "Input"],
Cell["Just for comparison:", "Text"],
Cell[BoxData[
\(vSolution[100. ] // MatrixForm\)], "Input"],
Cell[BoxData[
\(vSolution[\(-100. \)] // MatrixForm\)], "Input"],
Cell["If you have a nonhomogeneous system:", "Text"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"2", "3"},
{"5", "0"}
}], ")"}]}], ";", \(vX0 = {1, 1}\),
";", \(vX[t_] = {x1[t], x2[t]}\), ";", \(vC = {Sin[t], Cos[t]}\),
";"}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t] + vC\)], "Input"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. First[
DSolve[{system, vX[0] \[Equal] vX0}, vX[t], t]]\)], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["Numeric approxiamtions", "Subsection"],
Cell[TextData[{
"Or, we could have ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" do all the nasty approximations behind our backs and use NDSolve:"
}], "Text"],
Cell[BoxData[
\(Clear[mA, vX0, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"2", "3", \(-2\), "1"},
{"5", "0", "6", \(-1\)},
{"7", "5", "9", \(-7\)},
{"8", "4", "2", "1"}
}], ")"}]}], ";", \(vX0 = {1, 1, 1, 1}\),
";", \(vX[t_] = {x1[t], x2[t], x3[t], x4[t]}\), ";"}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t]\)], "Input"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. First[
NDSolve[{system, vX[0] \[Equal] vX0},
vX[t], {t, 0, 200}]]\)], "Input"],
Cell[TextData[{
"This error message is a warning that you can't trust this approximation \
beyond about ",
Cell[BoxData[
\(TraditionalForm\`t = 21\)]],
". To fix that, you must increase the maximum number of steps (which will \
take longer to compute):"
}], "Text"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. First[
NDSolve[{system, vX[0] \[Equal] vX0}, vX[t], {t, 0, 200},
MaxSteps \[Rule] 10000]]\)], "Input"],
Cell["Now, as before:", "Text"],
Cell[BoxData[
\(vSolution[100. ] // MatrixForm\)], "Input"],
Cell[BoxData[
\(vSolution[\(-100. \)] // MatrixForm\)], "Input"],
Cell[TextData[{
"You had better pay attention to this warning. What is the problem? We \
told NDSolve to solve on the domain ",
Cell[BoxData[
\(TraditionalForm\`0 \[LessEqual] t \[LessEqual] 200\)]],
", ",
StyleBox["not",
FontSlant->"Italic"],
" for negative ",
StyleBox["t",
FontSlant->"Italic"],
". This is something you need to watch out for."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Vector fields/direction fields/phase portraits", "Subsection"],
Cell["\<\
When you are working with systems of differential equations, it is often \
useful to look at either the vector field (direction field) of the system or \
the parametric plot of the solutions (or to just plot the individual \
solutions vs. time). Let's take a look at how to do this with 2-dimensional \
systems (some of this is applicable to 3 dimensional systems as well).
(Sections with an * are optional.)\
\>", "Text"],
Cell[CellGroupData[{
Cell["Vector fields", "Subsubsection"],
Cell["\<\
The vector field for a system of differential equations (first order, \
autonomous) is pretty easy to plot using PlotVectorField:\
\>", "Text"],
Cell[BoxData[
\(<< Graphics`\)], "Input"],
Cell[BoxData[
\(Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"1", \(-2\)},
{"2", "1"}
}], ")"}]}], ";", \(vX[t_] = {x[t], y[t]}\)}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t]\)], "Input"],
Cell[BoxData[
\(PlotVectorField[mA . {x, y}, {x, \(-10\), 10}, {y, \(-10\), 10},
Axes \[Rule] True]\)], "Input"],
Cell["\<\
(Notice, you don't have to actually solve the system to plot this.) If you \
want to plot the \"direction field\", you just need to arrange that all the \
vectors have the same length:\
\>", "Text"],
Cell[BoxData[
\(directionField =
PlotVectorField[mA . {x, y}, {x, \(-10\), 10}, {y, \(-10\), 10},
Axes \[Rule] True, ScaleFunction \[Rule] \((0.4 &)\),
ScaleFactor \[Rule] 0.9]\)], "Input"],
Cell[TextData[{
"Notice that these are plots in the XY plane (i.e., time doesn't enter into \
this). So, what would it mean to plot solution curves on this plane? It \
would mean to graph the parametric equations: ",
Cell[BoxData[
\(TraditionalForm\`{x \((t)\), y \((t)\)}\)]],
", graphing x vs. y in the plane (the \"trajectory\" in the XY plane)."
}], "Text"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. DSolve[{system, vX[0] \[Equal] {x0, y0}}, vX[t], t] //
Flatten\)], "Input"],
Cell[TextData[{
"If ",
Cell[BoxData[
\(TraditionalForm\`x(0) = \(2\ and\ \(y(0)\) = \(-1\)\)\)]],
", then we have:"
}], "Text"],
Cell[BoxData[
\(x0 = 2; y0 = \(-1\)\)], "Input"],
Cell[BoxData[
\(ParametricPlot[Evaluate[vSolution[t]], {t, 0, 2 \[Pi]}]\)], "Input"],
Cell["To show a bunch of curves:", "Text"],
Cell[BoxData[
\(y0 = 1; Clear[x0]\)], "Input"],
Cell[BoxData[
\(solutionCurves = Table[vSolution[t], {x0, \(-5\), 5, 1}]\)], "Input"],
Cell[BoxData[
\(solutionPlot =
ParametricPlot[Evaluate[solutionCurves], {t, \(-2\) \[Pi], 2 \[Pi]},
PlotStyle \[Rule] Blue]\)], "Input"],
Cell[BoxData[
\(Show[directionField, solutionPlot,
PlotRange \[Rule] {{\(-10\), 10}, {\(-10\), 10}}]\)], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["Parametric plot vs. time plots", "Subsubsection"],
Cell["\<\
In interpreting the plots above, remember that these are parametric plots of \
X vs. Y. What if you want to see how the two components (X and Y) behave \
over time?\
\>", "Text"],
Cell[BoxData[
\(Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"1", \(-2\)},
{"2", "1"}
}], ")"}]}], ";", \(vX[t_] = {x[t], y[t]}\)}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t]\)], "Input"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. DSolve[{system, vX[0] \[Equal] {x0, y0}}, vX[t], t] //
Flatten\)], "Input"],
Cell[BoxData[
\(x0 = 2; y0 = \(-1\)\)], "Input"],
Cell[BoxData[
\(ParametricPlot[Evaluate[vSolution[t]], {t, 0, 5}]\)], "Input"],
Cell[BoxData[
\(xTimePlot =
Plot[Evaluate[\(vSolution[t]\)[\([1]\)]], {t, 0, 5},
PlotStyle \[Rule] Red]\)], "Input"],
Cell[BoxData[
\(yTimePlot =
Plot[Evaluate[\(vSolution[t]\)[\([2]\)]], {t, 0, 5},
PlotStyle \[Rule] Green]\)], "Input"],
Cell[BoxData[
\(Show[xTimePlot, yTimePlot]\)], "Input"],
Cell["\<\
(To get the first coordinate of vSolution[t], you use vSolution[t][[1]], and \
similarly for other coordinates.)How does this compare with the parametric \
plot shown earlier?\
\>", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Nonhomogeneous systems and higher dimensional systems*", "Subsubsection"],
Cell["If your system is nonhomogeneous:", "Text"],
Cell[BoxData[
\(Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution, vC]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"mA", "=",
RowBox[{"(", GridBox[{
{"1", \(-2\)},
{"2", "1"}
}], ")"}]}], ";", \(vX[t_] = {x[t], y[t]}\),
";", \(vC = {t, t\^2}\), ";"}]], "Input"],
Cell[BoxData[
\(system = \(vX'\)[t] \[Equal] mA . vX[t] + vC\)], "Input"],
Cell[TextData[{
"The problem is that this system doesn't just vary with X and Y, it also \
varies with ",
StyleBox["t",
FontSlant->"Italic"],
". If you think of ",
StyleBox["t",
FontSlant->"Italic"],
" as time (which it often is), however, you can use things like animation \
to show this. One way to animate a graph is to make a table of \"frames\", \
which you can then double-click on and animate. So, if we use one frame for \
every \"tick\" of the clock we can make some interesting vector fields (we \
let ",
StyleBox["t ",
FontSlant->"Italic"],
"go from 0 to 10, with a stepsize of 1 here):"
}], "Text"],
Cell[BoxData[
\(Table[
PlotVectorField[
Evaluate[mA . {x, y} + vC], {x, \(-10\), 10}, {y, \(-10\), 10},
Axes \[Rule] True, ScaleFunction \[Rule] \((0.4 &)\),
ScaleFactor \[Rule] 0.9], {t, 0, 10, 1}]\)], "Input"],
Cell[TextData[{
"While this is pretty cool, I'm not sure it is really all that useful. \
(You can't, for example, use one of these graphs to trace out a trajectory. \
You would have to trace it out in real-time while it is being animated...)\n\n\
Notice, that we don't have to worry about this when graphing the parametric \
equations: ",
Cell[BoxData[
\(TraditionalForm\`{x(t), y(t)}\)]],
", graphing x vs. y in the plane (the \"trajectory\" in the XY plane). It \
is just another parametric equation (though a bit ugly...)"
}], "Text"],
Cell[BoxData[
\(vSolution[t_] =
vX[t] /. DSolve[{system, vX[0] \[Equal] {x0, y0}}, vX[t], t] //
Flatten\)], "Input"],
Cell[TextData[{
"If ",
Cell[BoxData[
\(TraditionalForm\`x(0) = \(2\ and\ \(y(0)\) = \(-1\)\)\)]],
", then we have:"
}], "Text"],
Cell[BoxData[
\(x0 = 2; y0 = \(-1\);\)], "Input"],
Cell[BoxData[
\(ParametricPlot[Evaluate[vSolution[t]], {t, 0, 2 \[Pi]}]\)], "Input"],
Cell[TextData[{
"Notice how this has been \"distorted\" from the homogeneous case above.\n\n\
So, what if you have a system of more than 2 dimensions? Well, if it is \
three dimensional, you could use the 3-dimensional graphing capabilities of \
",
StyleBox["Mathematica",
FontSlant->"Italic"],
" to draw a vector field or trajectories. (This may or may not be useful \
to you.) Another possibility is to take the components 2 at a time and treat \
them like a 2-dimensional case (to tell you something at least, of how those \
2 components interact). I may put up some examples of this later if I have \
time."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["\<\
Phase portraits of numeric approximations (think nonlinear systems)\
\>", "Subsubsection"],
Cell[TextData[{
"Unfortunately, many nonlinear systems of differential equations can't be \
solved (by ",
StyleBox["Mathematica",
FontSlant->"Italic"],
", at least) in any reasonable sort of manner. You can use the critical \
points of the system (we are talking mainly about 2-dimensional systems here) \
along with the eigenvalues of the linear approximaiton to the system and its \
phase portrait to analyze these systems. While the material given above is \
theoretically sufficient to draw these things, it turns out there are some \
annoying little issues that can really cuase problems if you aren't careful \
(and sometimes even if you are).\n\nSo, let's start with the following \
nonlinear autonomous system:\n\n",
Cell[BoxData[
FormBox[GridBox[{
{\(x' = x\^2 - 2 y\^2 - 2\)},
{\(y' = \(-4\) x\^2 + 3 y\^2\)}
}], TraditionalForm]]],
" or ",
Cell[BoxData[
FormBox[GridBox[{
{\(x' = f(x, y)\)},
{\(y' = g(x, y)\)}
}], TraditionalForm]]],
" where ",
Cell[BoxData[
\(TraditionalForm\`f(x, y) = x\^2 - 2 y\^2 - 2\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`g(x, y) = \(-4\) x\^2 + 3 y\^2\)]]
}], "Text"],
Cell[BoxData[
\(f[x_, y_] := x\^2 + 2 y\^2 - 2;
g[x_, y_] := \(-4\) x\^2 + 3 y\^2\)], "Input"],
Cell["Let's find our critical points:", "Text"],
Cell[BoxData[
\(critPts =
Solve[{f[x, y] \[Equal] 0, g[x, y] \[Equal] 0}] // N\)], "Input"],
Cell["\<\
Now, to classify these critical points, we need to evaluate the Jacobian at \
each of them and find its Eigenvalues (at the origin). So, we need to define \
a Jacobian first:\
\>", "Text"],
Cell[BoxData[
RowBox[{\(jacob[f_, g_]\), ":=",
RowBox[{"(", GridBox[{
{\(\[PartialD]\_x\ f\), \(\[PartialD]\_y\ f\)},
{\(\[PartialD]\_x\ g\), \(\[PartialD]\_y\ g\)}
}], ")"}]}]], "Input"],
Cell[BoxData[
\(jacob[f[x, y], g[x, y]]\)], "Input"],
Cell["\<\
At each critical point, you need to find the Jacobian and then find the \
eigenvalues of that matrix. So, for example, at the first critical point:\
\>", "Text"],
Cell[BoxData[
\(jacob[f[x, y], g[x, y]] /. critPts[\([1]\)] // MatrixForm\)], "Input"],
Cell[BoxData[
\(Eigenvalues[jacob[f[x, y], g[x, y]]] /. critPts[\([1]\)]\)], "Input"],
Cell["\<\
So, telll me, what does it mean that these critical points are imaginary with \
negative real part? (Chcck your answer down below when we graph the phase \
portrait.) Now, it is actually pretty easy to do this for all of the \
critical points at once (each row is a pair of eigenvalues, listed in the \
same order as the critical points in the list above):\
\>", "Text"],
Cell[BoxData[
\(Eigenvalues[jacob[f[x, y], g[x, y]]] /. critPts // TableForm\)], "Input"],
Cell["\<\
Take a moment to go through and classify each critical point. Now, let's \
graph the vector field (well, the direction field, to be precise) for this \
system:\
\>", "Text"],
Cell[BoxData[
\(<< Graphics`\)], "Input"],
Cell[BoxData[
\(vField =
PlotVectorField[{f[x, y], g[x, y]}, {x, \(-2\), 2}, {y, \(-2\), 2},
Axes \[Rule] True, ScaleFunction \[Rule] \((0.4 &)\),
ScaleFactor \[Rule] 0.2,
ColorFunction \[Rule] \((Green\ &)\)]\)], "Input"],
Cell["\<\
(Notice that we choose our domain to plot over so we can see the critical \
points and the region near them. If your vectors are too long or too short, \
you can change the value of the ScaleFactor option.)
Let's graph our critical points:\
\>", "Text"],
Cell[BoxData[
\(critPtsPlot =
ListPlot[{x, y} /. critPts,
PlotStyle -> {Red, PointSize[0.02]}]\)], "Input"],
Cell["\<\
(Notice two things here: We use ListPlot since we just want to plot some \
points, and we choose a PointSize of about 2% to make the points big enough \
to see easily.)
Sometimes, it is also handy to graph your \"null-clines\" as well:\
\>", "Text"],
Cell[BoxData[
\(nullClines = {f[x, y] \[Equal] 0, g[x, y] \[Equal] 0}\)], "Input"],
Cell[BoxData[
\(nullClinePlot =
ImplicitPlot[nullClines, {x, \(-2\), 2}, {y, \(-2\), 2},
PlotStyle \[Rule] {Blue, Yellow}, PlotPoints \[Rule] 100]\)], "Input"],
Cell["Now, it's time to put them all together:", "Text"],
Cell[BoxData[
\(Show[vField, nullClinePlot, critPtsPlot]\)], "Input"],
Cell[TextData[{
"This is pretty cool, but it would also be really helpful to add in some \
trajectories as well. Unfortunately, there is really not a nice automatic \
way to pick out the most useful trajectories to plot. Basically, you have to \
choose these individually (unless you just get lucky). Basically, you pick \
points in interesting regions of the graph to be initial conditions. Of \
course, since that means those would be the starting points of the graph, you \
might want to actually plot them for some negative ",
StyleBox["t",
FontSlant->"Italic"],
" values as well as postive (i.e., backtrack along the curves to see where \
they came from). \n\nNow, if the system can be solved symbolically in a \
reasonable amount of time, this isn't all that difficult (see the examples \
earlier in this notebook). However, if we try to solve this one, we get:"
}], "Text"],
Cell[BoxData[
\(DSolve[{\(x'\)[t] \[Equal] f[x[t], y[t]], \(y'\)[t] \[Equal]
g[x[t], y[t]], x[0] \[Equal] x0, y[0] \[Equal] y0}, {x[t], y[t]},
t]\)], "Input"],
Cell[TextData[{
"So no dice there. Clearly, we need to use NDSolve to find the solution. \
So, for example, if we want to plot a curve with initial condition ",
Cell[BoxData[
\(TraditionalForm\`\((1, 0)\)\)]],
" we could do:"
}], "Text"],
Cell[BoxData[
\(ParametricPlot[
Evaluate[{x[t], y[t]} /.
First[NDSolve[{\(x'\)[t] \[Equal] f[x[t], y[t]], \(y'\)[t] \[Equal]
g[x[t], y[t]], x[0] \[Equal] 1, y[0] \[Equal] 0}, {x[t],
y[t]}, {t, 0, 2}]]], {t, 0, 2}]\)], "Input"],
Cell[TextData[{
"This sounds good. But what about through the point ",
Cell[BoxData[
\(TraditionalForm\`\((1, \(-1\))\)\)]],
"?"
}], "Text"],
Cell[BoxData[
\(ParametricPlot[
Evaluate[{x[t], y[t]} /.
First[NDSolve[{\(x'\)[t] \[Equal] f[x[t], y[t]], \(y'\)[t] \[Equal]
g[x[t], y[t]], x[0] \[Equal] 1,
y[0] \[Equal] \(-1\)}, {x[t], y[t]}, {t, 0, 2}]]], {t, 0,
2}]\)], "Input"],
Cell[TextData[{
"Uh, I think we have a problem here. If you read the error messages, it is \
warning you that it can't compute a good solution all the way out to ",
Cell[BoxData[
\(TraditionalForm\`t = 2\)]],
" (it actually doesn't get much beyond ",
Cell[BoxData[
\(TraditionalForm\`t = 0.5\)]],
"), so it tells you it is extrapolating to get all the way to 2. As you \
can probably guess from looking at the numbers on the axes, this was pretty \
wildly unsuccessful. You could fix this by simply graphing from ",
Cell[BoxData[
\(TraditionalForm\`t = 0\)]],
" to ",
Cell[BoxData[
\(TraditionalForm\`t = 0.5\)]],
":"
}], "Text"],
Cell[BoxData[
\(ParametricPlot[
Evaluate[{x[t], y[t]} /.
First[NDSolve[{\(x'\)[t] \[Equal] f[x[t], y[t]], \(y'\)[t] \[Equal]
g[x[t], y[t]], x[0] \[Equal] 1,
y[0] \[Equal] \(-1\)}, {x[t], y[t]}, {t, 0, 0.5}]]], {t, 0,
0.5}]\)], "Input"],
Cell[TextData[{
"This works okay if you are only graphing a few curves, but it would really \
be nice to supply a list of points and then have ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" graph the solution curves for all of them, without having to set separate \
domains for individual curves. The best way to do this would be to just \
instruct ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" to not extrapolate approximate functions beyond their domain. \
Unfortunately, I haven't been able to figure out a way to do that (except on \
a case by case basis). \n\nSo, after spending way too much time on this over \
the weekend, I decided to just write my own function to take care of this. \
This was really annoying. The good news for you is that you don't really \
have to worry about it. You can just use the function I wrote (think \
cut-and-paste) and everything works out okay. (You owe me for this one...) \
Now, since I really didn't want to spend even more of my weekend making this \
really fancy, this means that you have to use this function ",
StyleBox["exactly",
FontSlant->"Italic"],
" as I show you; I haven't built in any fancy error-handling or usage notes \
or anything else. If you prefer not to use it, you can just plot individual \
curves (naming each one of them) and then combine them all together using the \
Show command. On the other hand, if you prefer to take advantage of my \
function, I list it below (don't forget to execute it); do ",
StyleBox["not",
FontSlant->"Italic"],
" change ",
StyleBox["anything",
FontSlant->"Italic"],
" in it (unless you really know ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" pretty well) or you might break it. (Unless you are just interested, it \
really isn't at all important that you understand how it works. Actually, if \
you need to, you could change the value for PlotPoints, but be careful.)"
}], "Text"],
Cell[BoxData[
\(TrajectoryPlot[system_,
initialPts_, {time_, tMin_, tMax_}, {xVar_, xMin_, xMax_}, {yVar_,
yMin_, yMax_}] :=
Module[{}, \[IndentingNewLine]Off[
NDSolve::ndsz]; \[IndentingNewLine]Off[
InterpolatingFunction::dmval]; \[IndentingNewLine]\
\[IndentingNewLine]getSolution[k_] :=
First[NDSolve[
Join[system, {x[0] \[Equal] initialPts[\([k, 1]\)],
y[0] \[Equal] initialPts[\([k, 2]\)]}], {x, y}, {t, tMin,
tMax}]]; \[IndentingNewLine]\[IndentingNewLine]getEstDomain[
funct_] := {t, \(Flatten[funct[\([1]\)]]\)[\([1]\)] +
0.2, \(Flatten[funct[\([1]\)]]\)[\([2]\)] -
0.2}; \[IndentingNewLine]\[IndentingNewLine]Block[{$\
DisplayFunction = Identity}, \[IndentingNewLine]\(solutionList =
Table[ParametricPlot[Evaluate[{x[t], y[t]} /. getSolution[k]],
Evaluate[getEstDomain[x /. getSolution[k]]],
PlotRange \[Rule] {{xMin, xMax}, {yMin, yMax}},
AspectRatio \[Rule] Automatic, PlotPoints \[Rule] 100], {k,
1, Length[initialPts]}];\)]; \[IndentingNewLine]On[
NDSolve::ndsz]; \[IndentingNewLine]On[
InterpolatingFunction::dmval]; \[IndentingNewLine]\
\[IndentingNewLine]Show[solutionList]\[IndentingNewLine]]\)], "Input"],
Cell[TextData[{
"So, how do you use it? Fortunately, it is very simple. You just have to \
supply it with the system (without any sort of initial conditions), a list of \
points you want curves through, the domain you want ",
StyleBox["t",
FontSlant->"Italic"],
" to vary over, and the domains for ",
StyleBox["x",
FontSlant->"Italic"],
" and ",
StyleBox["y",
FontSlant->"Italic"],
". Like this:"
}], "Text"],
Cell[BoxData[
\(points = {{1, 1}, {1, \(-1\)}, {\(-1\),
1}, {\(-1\), \(-1\)}}\)], "Input"],
Cell[BoxData[
\(solCurves =
TrajectoryPlot[{\(x'\)[t] == f[x[t], y[t]], \(y'\)[t] ==
g[x[t], y[t]]},
points, {t, \(-100\), 100}, {x, \(-2\), 2}, {y, \(-2\),
2}]\)], "Input"],
Cell["Now, let's put all our graphs together:", "Text"],
Cell[BoxData[
\(Show[vField, nullClinePlot, critPtsPlot, solCurves,
PlotRange \[Rule] {{\(-2\), 2}, {\(-2\), 2}}]\)], "Input"],
Cell["\<\
Let's skip the null-clines and add in some more points (you should think \
about what points to choose to make things clearer):\
\>", "Text"],
Cell[BoxData[
\(\(points = {{1, 1}, {1, \(-1\)}, {\(-1\),
1}, {\(-1\), \(-1\)}, {0, .25}, {\(-1\), \(-2\)}, {0, \(-2\)}, \
{0, .8}, {\(-1\), 1.5}, {1.5, 1.5}, {\(- .5\), 1}, { .5, \(-2\)}, {\(-1\),
1.2}, {0.5, 0}};\)\)], "Input"],
Cell[BoxData[
\(solCurves =
TrajectoryPlot[{\(x'\)[t] == f[x[t], y[t]], \(y'\)[t] ==
g[x[t], y[t]]},
points, {t, \(-100\), 100}, {x, \(-2\), 2}, {y, \(-2\),
2}]\)], "Input"],
Cell[BoxData[
\(Show[vField, critPtsPlot, solCurves,
PlotRange \[Rule] {{\(-2\), 2}, {\(-2\), 2}}]\)], "Input"],
Cell[TextData[{
"Or, if you're just lazy (like me), you could have ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" plug in a grid of points for you (and if you're lucky, it will look \
good):"
}], "Text"],
Cell[BoxData[
\(points =
Flatten[Table[{x, y}, {x, \(-2\), 2, 0.5}, {y, \(-2\), 2, 0.5}],
1]\)], "Input"],
Cell["\<\
(You need the Flatten command here to knock out some excess {}'s.)\
\>", "Text"],
Cell[BoxData[
\(solCurves =
TrajectoryPlot[{\(x'\)[t] == f[x[t], y[t]], \(y'\)[t] ==
g[x[t], y[t]]},
points, {t, \(-100\), 100}, {x, \(-2\), 2}, {y, \(-2\),
2}]\)], "Input"],
Cell[BoxData[
\(Show[vField, critPtsPlot, solCurves,
PlotRange \[Rule] {{\(-2\), 2}, {\(-2\), 2}}]\)], "Input"],
Cell[TextData[{
"Well, you get the idea. The combination of the direction field with the \
solution curves actually gives you a pretty decent idea of what's going on \
here. Of course, you still have to figure out what that ",
StyleBox["means, ",
FontSlant->"Italic"],
"but that depends on the specific application."
}], "Text"]
}, Closed]]
}, Closed]]
}, Open ]]
}, Open ]]
},
FrontEndVersion->"5.0 for Macintosh",
ScreenRectangle->{{0, 1024}, {0, 705}},
AutoGeneratedPackage->None,
WindowSize->{856, 533},
WindowMargins->{{61, Automatic}, {61, Automatic}},
CellLabelAutoDelete->True
]
(*******************************************************************
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[1754, 51, 73, 2, 46, "Input",
InitializationCell->True],
Cell[CellGroupData[{
Cell[1852, 57, 132, 5, 126, "Subtitle"],
Cell[CellGroupData[{
Cell[2009, 66, 48, 0, 121, "Section"],
Cell[CellGroupData[{
Cell[2082, 70, 40, 0, 62, "Subsection"],
Cell[2125, 72, 1058, 25, 203, "Text"],
Cell[CellGroupData[{
Cell[3208, 101, 42, 0, 61, "Subsubsection"],
Cell[3253, 103, 105, 2, 42, "Input"],
Cell[3361, 107, 1157, 23, 281, "Text"],
Cell[4521, 132, 53, 1, 42, "Input"],
Cell[4577, 135, 196, 4, 73, "Text"],
Cell[4776, 141, 91, 1, 42, "Input"],
Cell[4870, 144, 155, 4, 68, "Input"],
Cell[5028, 150, 582, 13, 181, "Text"],
Cell[5613, 165, 74, 1, 42, "Input"],
Cell[5690, 168, 240, 4, 73, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[5967, 177, 139, 3, 60, "Subsubsection"],
Cell[6109, 182, 1204, 23, 515, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[7350, 210, 47, 0, 38, "Subsubsection"],
Cell[7400, 212, 532, 12, 125, "Text"],
Cell[7935, 226, 59, 1, 42, "Input"],
Cell[7997, 229, 91, 1, 42, "Input"],
Cell[8091, 232, 216, 5, 120, "Input"],
Cell[8310, 239, 126, 4, 99, "Text"],
Cell[8439, 245, 64, 1, 42, "Input"],
Cell[8506, 248, 258, 5, 99, "Text"],
Cell[8767, 255, 46, 1, 42, "Input"],
Cell[8816, 258, 200, 4, 120, "Input"],
Cell[9019, 264, 117, 3, 42, "Input"],
Cell[9139, 269, 326, 8, 99, "Text"],
Cell[9468, 279, 46, 1, 42, "Input"],
Cell[9517, 282, 65, 1, 42, "Input"],
Cell[9585, 285, 104, 2, 42, "Input"],
Cell[9692, 289, 473, 12, 99, "Text"],
Cell[10168, 303, 46, 1, 42, "Input"],
Cell[10217, 306, 65, 1, 42, "Input"],
Cell[10285, 309, 104, 2, 42, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[10426, 316, 43, 0, 38, "Subsubsection"],
Cell[10472, 318, 267, 6, 73, "Text"],
Cell[10742, 326, 70, 1, 42, "Input"],
Cell[10815, 329, 79, 1, 42, "Input"],
Cell[10897, 332, 192, 5, 120, "Input"],
Cell[11092, 339, 379, 10, 73, "Text"],
Cell[11474, 351, 73, 1, 42, "Input"],
Cell[11550, 354, 74, 1, 42, "Input"],
Cell[11627, 357, 139, 5, 47, "Text"],
Cell[11769, 364, 83, 1, 42, "Input"],
Cell[11855, 367, 74, 1, 42, "Input"]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[11978, 374, 44, 0, 45, "Subsection"],
Cell[12025, 376, 1631, 39, 463, "Text"],
Cell[13659, 417, 50, 1, 42, "Input"],
Cell[13712, 420, 1402, 31, 172, "Input"],
Cell[15117, 453, 41, 1, 42, "Input"],
Cell[15161, 456, 58, 1, 42, "Input"],
Cell[15222, 459, 533, 8, 203, "Text"],
Cell[15758, 469, 1453, 32, 172, "Input"],
Cell[17214, 503, 121, 3, 73, "Text"],
Cell[17338, 508, 1701, 37, 198, "Input"],
Cell[19042, 547, 74, 1, 42, "Input"],
Cell[19119, 550, 558, 8, 229, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[19714, 563, 57, 0, 45, "Subsection"],
Cell[CellGroupData[{
Cell[19796, 567, 37, 0, 61, "Subsubsection"],
Cell[19836, 569, 468, 12, 107, "Text"],
Cell[20307, 583, 45, 1, 42, "Input"],
Cell[20355, 586, 117, 2, 64, "Input"],
Cell[20475, 590, 95, 2, 47, "Text"],
Cell[20573, 594, 136, 2, 95, "Input"],
Cell[20712, 598, 107, 3, 47, "Text"],
Cell[20822, 603, 227, 4, 122, "Input"],
Cell[21052, 609, 445, 11, 99, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[21534, 625, 76, 0, 38, "Subsubsection"],
Cell[21613, 627, 616, 13, 203, "Text"],
Cell[22232, 642, 75, 1, 42, "Input"],
Cell[22310, 645, 217, 5, 124, "Input"],
Cell[22530, 652, 133, 3, 68, "Input"],
Cell[22666, 657, 76, 0, 47, "Text"],
Cell[22745, 659, 64, 1, 42, "Input"],
Cell[22812, 662, 116, 3, 47, "Text"],
Cell[22931, 667, 121, 2, 42, "Input"]
}, Closed]]
}, Closed]]
}, Open ]],
Cell[CellGroupData[{
Cell[23113, 676, 52, 0, 121, "Section"],
Cell[23168, 678, 619, 15, 71, "Text"],
Cell[CellGroupData[{
Cell[23812, 697, 68, 0, 47, "Subsection"],
Cell[23883, 699, 262, 7, 73, "Text"],
Cell[24148, 708, 139, 2, 68, "Input"],
Cell[24290, 712, 847, 21, 229, "Text"],
Cell[25140, 735, 80, 1, 42, "Input"],
Cell[25223, 738, 619, 13, 177, "Text"],
Cell[25845, 753, 301, 10, 59, "Input"],
Cell[26149, 765, 19, 0, 47, "Text"],
Cell[26171, 767, 318, 11, 59, "Input"],
Cell[26492, 780, 52, 1, 42, "Input"],
Cell[26547, 783, 997, 21, 281, "Text"],
Cell[27547, 806, 91, 2, 42, "Input",
InitializationCell->True],
Cell[27641, 810, 555, 12, 192, "Text"],
Cell[28199, 824, 277, 8, 83, "Input"],
Cell[28479, 834, 297, 7, 151, "Text"],
Cell[28779, 843, 61, 1, 42, "Input"],
Cell[28843, 846, 46, 1, 42, "Input"],
Cell[28892, 849, 366, 11, 136, "Text"],
Cell[29261, 862, 79, 1, 42, "Input"],
Cell[29343, 865, 208, 6, 59, "Input"],
Cell[29554, 873, 72, 1, 42, "Input"],
Cell[29629, 876, 82, 1, 42, "Input"],
Cell[29714, 879, 92, 2, 47, "Text"],
Cell[29809, 883, 47, 1, 42, "Input"],
Cell[29859, 886, 111, 2, 42, "Input"],
Cell[29973, 890, 73, 0, 47, "Text"],
Cell[30049, 892, 83, 1, 42, "Input"],
Cell[30135, 895, 237, 7, 59, "Input"],
Cell[30375, 904, 77, 1, 42, "Input"],
Cell[30455, 907, 82, 1, 42, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[30574, 913, 60, 0, 31, "Subsection"],
Cell[30637, 915, 571, 11, 151, "Text"],
Cell[31211, 928, 81, 1, 42, "Input"],
Cell[31295, 931, 329, 9, 105, "Input"],
Cell[31627, 942, 72, 1, 42, "Input"],
Cell[31702, 945, 187, 4, 73, "Text"],
Cell[31892, 951, 63, 1, 42, "Input"],
Cell[31958, 954, 65, 1, 42, "Input"],
Cell[32026, 957, 373, 6, 125, "Text"],
Cell[32402, 965, 111, 3, 43, "Input"],
Cell[32516, 970, 111, 3, 43, "Input"],
Cell[32630, 975, 119, 3, 47, "Text"],
Cell[32752, 980, 262, 5, 95, "Input"],
Cell[33017, 987, 262, 5, 95, "Input"],
Cell[33282, 994, 181, 5, 99, "Text"],
Cell[33466, 1001, 113, 2, 42, "Input"],
Cell[33582, 1005, 59, 1, 42, "Input"],
Cell[33644, 1008, 55, 0, 47, "Text"],
Cell[33702, 1010, 51, 1, 42, "Input"],
Cell[33756, 1013, 47, 0, 47, "Text"],
Cell[33806, 1015, 82, 1, 42, "Input"],
Cell[33891, 1018, 302, 5, 99, "Text"],
Cell[34196, 1025, 101, 2, 42, "Input"],
Cell[34300, 1029, 622, 12, 151, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[34959, 1046, 109, 3, 31, "Subsection"],
Cell[35071, 1051, 4884, 143, 344, "Text"],
Cell[39958, 1196, 237, 7, 59, "Input"],
Cell[40198, 1205, 93, 2, 47, "Text"],
Cell[40294, 1209, 58, 1, 42, "Input"],
Cell[40355, 1212, 308, 8, 105, "Input"],
Cell[40666, 1222, 91, 1, 42, "Input"],
Cell[40760, 1225, 106, 3, 47, "Text"],
Cell[40869, 1230, 59, 1, 42, "Input"],
Cell[40931, 1233, 262, 6, 73, "Text"],
Cell[41196, 1241, 63, 1, 42, "Input"],
Cell[41262, 1244, 130, 5, 47, "Text"],
Cell[41395, 1251, 68, 1, 42, "Input"],
Cell[41466, 1254, 124, 3, 47, "Text"],
Cell[41593, 1259, 88, 1, 42, "Input"],
Cell[41684, 1262, 69, 0, 47, "Text"],
Cell[41756, 1264, 63, 1, 42, "Input"],
Cell[41822, 1267, 68, 1, 42, "Input"],
Cell[41893, 1270, 192, 4, 73, "Text"],
Cell[42088, 1276, 110, 2, 42, "Input"],
Cell[42201, 1280, 63, 1, 42, "Input"],
Cell[42267, 1283, 68, 1, 42, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[42372, 1289, 48, 0, 31, "Subsection"],
Cell[42423, 1291, 123, 3, 47, "Text"],
Cell[42549, 1296, 58, 1, 42, "Input"],
Cell[42610, 1299, 363, 9, 139, "Input"],
Cell[42976, 1310, 72, 1, 42, "Input"],
Cell[43051, 1313, 129, 3, 68, "Input"],
Cell[43183, 1318, 36, 0, 47, "Text"],
Cell[43222, 1320, 63, 1, 42, "Input"],
Cell[43288, 1323, 68, 1, 42, "Input"],
Cell[43359, 1326, 52, 0, 47, "Text"],
Cell[43414, 1328, 281, 8, 87, "Input"],
Cell[43698, 1338, 77, 1, 42, "Input"],
Cell[43778, 1341, 129, 3, 68, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[43944, 1349, 44, 0, 31, "Subsection"],
Cell[43991, 1351, 174, 5, 47, "Text"],
Cell[44168, 1358, 58, 1, 42, "Input"],
Cell[44229, 1361, 363, 9, 139, "Input"],
Cell[44595, 1372, 72, 1, 42, "Input"],
Cell[44670, 1375, 153, 4, 68, "Input"],
Cell[44826, 1381, 278, 7, 73, "Text"],
Cell[45107, 1390, 177, 4, 120, "Input"],
Cell[45287, 1396, 31, 0, 47, "Text"],
Cell[45321, 1398, 63, 1, 42, "Input"],
Cell[45387, 1401, 68, 1, 42, "Input"],
Cell[45458, 1404, 389, 12, 73, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[45884, 1421, 68, 0, 31, "Subsection"],
Cell[45955, 1423, 434, 8, 177, "Text"],
Cell[CellGroupData[{
Cell[46414, 1435, 38, 0, 61, "Subsubsection"],
Cell[46455, 1437, 153, 3, 73, "Text"],
Cell[46611, 1442, 45, 1, 42, "Input"],
Cell[46659, 1445, 79, 1, 42, "Input"],
Cell[46741, 1448, 208, 6, 59, "Input"],
Cell[46952, 1456, 72, 1, 42, "Input"],
Cell[47027, 1459, 123, 2, 68, "Input"],
Cell[47153, 1463, 209, 4, 73, "Text"],
Cell[47365, 1469, 218, 4, 94, "Input"],
Cell[47586, 1475, 376, 7, 99, "Text"],
Cell[47965, 1484, 136, 3, 68, "Input"],
Cell[48104, 1489, 139, 5, 47, "Text"],
Cell[48246, 1496, 52, 1, 42, "Input"],
Cell[48301, 1499, 89, 1, 42, "Input"],
Cell[48393, 1502, 42, 0, 47, "Text"],
Cell[48438, 1504, 50, 1, 42, "Input"],
Cell[48491, 1507, 89, 1, 42, "Input"],
Cell[48583, 1510, 157, 3, 68, "Input"],
Cell[48743, 1515, 124, 2, 68, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[48904, 1522, 55, 0, 38, "Subsubsection"],
Cell[48962, 1524, 189, 4, 73, "Text"],
Cell[49154, 1530, 79, 1, 42, "Input"],
Cell[49236, 1533, 208, 6, 59, "Input"],
Cell[49447, 1541, 72, 1, 42, "Input"],
Cell[49522, 1544, 136, 3, 68, "Input"],
Cell[49661, 1549, 52, 1, 42, "Input"],
Cell[49716, 1552, 82, 1, 42, "Input"],
Cell[49801, 1555, 136, 3, 68, "Input"],
Cell[49940, 1560, 138, 3, 68, "Input"],
Cell[50081, 1565, 59, 1, 42, "Input"],
Cell[50143, 1568, 199, 4, 73, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[50379, 1577, 79, 0, 38, "Subsubsection"],
Cell[50461, 1579, 49, 0, 47, "Text"],
Cell[50513, 1581, 83, 1, 42, "Input"],
Cell[50599, 1584, 245, 7, 59, "Input"],
Cell[50847, 1593, 77, 1, 42, "Input"],
Cell[50927, 1596, 639, 16, 151, "Text"],
Cell[51569, 1614, 247, 5, 94, "Input"],
Cell[51819, 1621, 553, 10, 177, "Text"],
Cell[52375, 1633, 136, 3, 68, "Input"],
Cell[52514, 1638, 139, 5, 47, "Text"],
Cell[52656, 1645, 53, 1, 42, "Input"],
Cell[52712, 1648, 89, 1, 42, "Input"],
Cell[52804, 1651, 637, 12, 203, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[53478, 1668, 100, 2, 38, "Subsubsection"],
Cell[53581, 1672, 1236, 30, 306, "Text"],
Cell[54820, 1704, 105, 2, 43, "Input"],
Cell[54928, 1708, 47, 0, 47, "Text"],
Cell[54978, 1710, 102, 2, 42, "Input"],
Cell[55083, 1714, 199, 4, 73, "Text"],
Cell[55285, 1720, 234, 5, 63, "Input"],
Cell[55522, 1727, 56, 1, 42, "Input"],
Cell[55581, 1730, 172, 3, 73, "Text"],
Cell[55756, 1735, 90, 1, 42, "Input"],
Cell[55849, 1738, 89, 1, 42, "Input"],
Cell[55941, 1741, 382, 6, 125, "Text"],
Cell[56326, 1749, 93, 1, 42, "Input"],
Cell[56422, 1752, 184, 4, 73, "Text"],
Cell[56609, 1758, 45, 1, 42, "Input"],
Cell[56657, 1761, 261, 5, 94, "Input"],
Cell[56921, 1768, 265, 6, 125, "Text"],
Cell[57189, 1776, 127, 3, 68, "Input"],
Cell[57319, 1781, 261, 6, 125, "Text"],
Cell[57583, 1789, 86, 1, 42, "Input"],
Cell[57672, 1792, 179, 3, 68, "Input"],
Cell[57854, 1797, 56, 0, 47, "Text"],
Cell[57913, 1799, 73, 1, 42, "Input"],
Cell[57989, 1802, 896, 14, 255, "Text"],
Cell[58888, 1818, 180, 3, 68, "Input"],
Cell[59071, 1823, 251, 6, 73, "Text"],
Cell[59325, 1831, 283, 5, 172, "Input"],
Cell[59611, 1838, 153, 5, 47, "Text"],
Cell[59767, 1845, 297, 6, 172, "Input"],
Cell[60067, 1853, 674, 17, 125, "Text"],
Cell[60744, 1872, 301, 6, 172, "Input"],
Cell[61048, 1880, 1957, 37, 463, "Text"],
Cell[63008, 1919, 1387, 23, 744, "Input"],
Cell[64398, 1944, 437, 13, 99, "Text"],
Cell[64838, 1959, 105, 2, 42, "Input"],
Cell[64946, 1963, 217, 5, 94, "Input"],
Cell[65166, 1970, 55, 0, 47, "Text"],
Cell[65224, 1972, 137, 2, 68, "Input"],
Cell[65364, 1976, 151, 3, 73, "Text"],
Cell[65518, 1981, 261, 4, 94, "Input"],
Cell[65782, 1987, 217, 5, 94, "Input"],
Cell[66002, 1994, 122, 2, 68, "Input"],
Cell[66127, 1998, 217, 6, 73, "Text"],
Cell[66347, 2006, 125, 3, 42, "Input"],
Cell[66475, 2011, 90, 2, 47, "Text"],
Cell[66568, 2015, 217, 5, 94, "Input"],
Cell[66788, 2022, 122, 2, 68, "Input"],
Cell[66913, 2026, 340, 7, 99, "Text"]
}, Closed]]
}, Closed]]
}, Open ]]
}, Open ]]
}
]
*)
(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)