(***********************************************************************
Mathematica-Compatible Notebook
This notebook can be used on any computer system with Mathematica 3.0,
MathReader 3.0, or any compatible application. The data for the notebook
starts with the line of 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[ 5667, 167]*)
(*NotebookOutlinePosition[ 6380, 192]*)
(* CellTagsIndexPosition[ 6336, 188]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[CellGroupData[{
Cell["Numerical solution of boundary value problems", "Title"],
Cell[CellGroupData[{
Cell["Linear, second order", "Section"],
Cell[TextData[{
" Although the standard Mathematica routine NDSolve is supposed to work for \
this class of problems, it is remarkable delicate, and fails often even if \
\[Epsilon] is not particularly small. It uses the idea of 'shooting'. Using \
for ex. the left boundary condition (BC) and guessing on a second left \
condition (say, for the first derivative) two different solutions are found \
numerically. Typically neither of these two will satisfy the right BC. \
However, because of linearity, one can explicitly find a linear conbination \
of the two trial solutions will do this.\n\t\t\nWe use below ",
StyleBox["Mathematica",
FontSlant->"Italic"],
"'s built-in solver for the example in B&O page 431. "
}], "Text"],
Cell[BoxData[
\(ode =
\[Epsilon]\ \(\(y'\)'\)[x] + \((1 + x)\)\ \(y'\)[x] + y[x] == 0\)],
"Input"],
Cell[BoxData[
\(sol =
NDSolve[{ode\ /. \[Epsilon] -> 0.1, y[0] == 1, y[1] == 1},
y[x], {x, 0, 1}]\)], "Input"],
Cell["\<\
The output of NDSolve takes the form of an InterpolatingFunction - \
essentially a table of intermediate values from which the function values at \
arbitrary points can be obtained by standard interpolation. We plot the \
result as follows\
\>", "Text"],
Cell[BoxData[
\(Plot[Evaluate[y[x]\ /. \ sol], {x, 0, 1}]\)], "Input"],
Cell[TextData[
"Mathematica fails for this problem if we try \[Epsilon]=0.05 - still not a \
very small value."], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Nonlinear, or higher order equations", "Section"],
Cell[TextData[{
"We can now implement basically the same principle ourselves. Instead of \
just using two tries, we have to allow several tries, and use Mathematica's \
root finding algorithm (for general nonlinear functions) in order to \
converge in on the solution. \n\nAs an example of this, we consider a third \
order ODE from B&O page 446 (linear, but third order; ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" will not accept the direct approach above):"
}], "Text"],
Cell[BoxData[
\(ode\ = \ \[Epsilon]\ \(\(\(y'\)'\)'\)[x] - \(y'\)[x] + x\ y[x] == 0\)],
"Input"],
Cell["\<\
with the BC y[0]=y'[0]=y[1]=1. Since two conditions are already given at \
x=0, we decide to there add a third; y''[0]=a. The problem has now become to \
choose the value of a so that y[1]=1. Given any value of a, we can compute \
the corresponding value of y[1] through the following function \
\>", "Text"],
Cell[BoxData[
\(r[a_] := \
\((sol = NDSolve[{ode\ /. \ \[Epsilon] -> 0.01, y[0] == 1,
\(y'\)[0] == 1, \(\(y'\)'\)[0] == a}, y[x], {x, 0, 1}]; \t\t
Evaluate[y[x] /. sol]\ /. x -> 1)\)\ [\([1]\)]\)], "Input"],
Cell["\<\
We could for example try a few values of the variable a 'manually' to see if \
we can get the result to become close to 1:\
\>", "Text"],
Cell[BoxData[
\(r[0]\)], "Input"],
Cell[BoxData[
\(r[\(-10\)]\)], "Input"],
Cell[BoxData[
\(r[\(-9\)]\)], "Input"],
Cell[BoxData[
\(r[\(-8.9\)]\)], "Input"],
Cell["\<\
Mathematica has a function FindRoot that does this kind of search for us. \
After specifying the equation to solve, FindRoot requires one or two starting \
guesses, and then conducts a local search. With one guess, it tries a \
Newton-type method - and fails in this type of cases. With two initial \
guesses (which normally don't have to be very good at all) it successfully \
tries a secant method. So\
\>", "Text"],
Cell[BoxData[
\(FindRoot[r[a] == 1, {a, 0, 1}]\)], "Input"],
Cell["\<\
FindRoot called r[a] several times, and we can still access the full solution \
generated by its last ODE evaluation. So we plot this\
\>", "Text"],
Cell[BoxData[
\(Plot[Evaluate[y[x]\ /. \ sol], {x, 0, 1}]\)], "Input"]
}, Closed]]
}, Open ]]
},
FrontEndVersion->"Microsoft Windows 3.0",
ScreenRectangle->{{0, 1024}, {0, 712}},
WindowSize->{706, 607},
WindowMargins->{{10, Automatic}, {Automatic, 15}},
PrintingCopies->1,
PrintingPageRange->{Automatic, Automatic}
]
(***********************************************************************
Cached data follows. If you edit this Notebook file directly, not using
Mathematica, you must remove the line containing CacheID at the top of
the file. The cache data will then be recreated when you save this file
from within Mathematica.
***********************************************************************)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[1731, 51, 62, 0, 150, "Title"],
Cell[CellGroupData[{
Cell[1818, 55, 39, 0, 53, "Section"],
Cell[1860, 57, 737, 12, 128, "Text"],
Cell[2600, 71, 112, 3, 30, "Input"],
Cell[2715, 76, 131, 3, 30, "Input"],
Cell[2849, 81, 263, 5, 52, "Text"],
Cell[3115, 88, 75, 1, 30, "Input"],
Cell[3193, 91, 121, 2, 33, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[3351, 98, 55, 0, 33, "Section"],
Cell[3409, 100, 491, 9, 109, "Text"],
Cell[3903, 111, 104, 2, 30, "Input"],
Cell[4010, 115, 319, 5, 71, "Text"],
Cell[4332, 122, 246, 4, 70, "Input"],
Cell[4581, 128, 146, 3, 33, "Text"],
Cell[4730, 133, 37, 1, 30, "Input"],
Cell[4770, 136, 43, 1, 30, "Input"],
Cell[4816, 139, 42, 1, 30, "Input"],
Cell[4861, 142, 44, 1, 30, "Input"],
Cell[4908, 145, 427, 7, 71, "Text"],
Cell[5338, 154, 63, 1, 30, "Input"],
Cell[5404, 157, 157, 3, 33, "Text"],
Cell[5564, 162, 75, 1, 30, "Input"]
}, Closed]]
}, Open ]]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)