Re: Something really wrong with FindRoot or NDSolve in Mathematica 5.0
- To: mathgroup at smc.vnet.net
- Subject: [mg44017] Re: Something really wrong with FindRoot or NDSolve in Mathematica 5.0
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Fri, 17 Oct 2003 05:14:51 -0400 (EDT)
- Organization: The University of Western Australia
- References: <bmlk3m$6fm$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In article <bmlk3m$6fm$1 at smc.vnet.net>, Arturas Acus <acus at itpa.lt>
wrote:
> Here is the example which worked really fine with 4.1
> and completely refuses to work with 5.0:
>
> teisingasKampas =
> FindRoot[
> (u[\[Xi]] /. First[sprendinys1 =
> NDSolve[{Derivative[2][\[Psi]][\[Xi]] == (\[Xi]^2 - energija)*
> \[Psi][\[Xi]], Derivative[1][u][\[Xi]] == \[Psi][\[Xi]]^2,
> \[Psi][\[Xi]min] == 0,
> Derivative[1][\[Psi]][\[Xi]min] == kampas, u[\[Xi]min] == 0},
> {\[Psi][\[Xi]], u[\[Xi]]},
> {\[Xi], \[Xi]min, \[Xi]max}]] /. \[Xi] -> \[Xi]max) == 1,
> {kampas, -0.04, 0.04},
> MaxIterations -> 30]
>
> It uses FindRoot to find shooting angle for derivative. In 5.0 Mathematica
> complains about undefined symbol value of kampas. I need to modify
> modify this code as little as possible to make it work.
The key is that you only want the argument of NDSolve to evaluate when
all parameters are numeric. The easiest (and cleanest) way to obtain
this behaviour is to break the computation down into small (easily
maintainable) chunks. Here is a Notebook at does this.
(************** Content-type: application/mathematica **************
CreatedBy='Mathematica 5.0'
*******************************************************************)
(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[ 4473, 129]*)
(*NotebookOutlinePosition[ 5111, 151]*)
(* CellTagsIndexPosition[ 5067, 147]*)
(*WindowFrame->Normal*)
Notebook[{
Cell["Break the code down into three routines:", "Text"],
Cell[BoxData[
\(TraditionalForm\`Clear[f, g, teisingasKampas]\)], "Input"],
Cell[BoxData[
FormBox[
RowBox[{\(g[energija_, kampas_?NumberQ, \[Xi]min_, \[Xi]max_]\),
":=",
RowBox[{"First", "[",
RowBox[{"NDSolve", "[",
RowBox[{
RowBox[{"{",
RowBox[{
RowBox[{
RowBox[{
SuperscriptBox["\[Psi]", "\[Prime]\[Prime]",
MultilineFunction->None], "(", "\[Xi]", ")"}],
"\[Equal]", \(\((\[Xi]\^2 -
energija)\)\ \(\[Psi](\[Xi])\)\)}], ",",
RowBox[{
RowBox[{
SuperscriptBox["u", "\[Prime]",
MultilineFunction->None], "(", "\[Xi]", ")"}],
"\[Equal]", \(\(\[Psi](\[Xi])\)\^2\)}],
",", \(\[Psi](\[Xi]min) \[Equal] 0\), ",",
RowBox[{
RowBox[{
SuperscriptBox["\[Psi]", "\[Prime]",
MultilineFunction->None], "(", "\[Xi]min",
")"}],
"\[Equal]", "kampas"}], ",", \(u(\[Xi]min) \[Equal]
0\)}],
"}"}], ",", \({\[Psi], u}\),
",", \({\[Xi], \[Xi]min, \[Xi]max}\)}], "]"}], "]"}]}],
TraditionalForm]], "Input"],
Cell[BoxData[
\(TraditionalForm\`f[energija_, kampas_?NumberQ, \[Xi]min_,
\[Xi]max_] :=
u(\[Xi]max) /. \[InvisibleSpace]g[energija,
kampas, \[Xi]min, \[Xi]max]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`teisingasKampas[energija_, \[Xi]min_, \[Xi]max_]
:=
FindRoot[f[energija, kampas, \[Xi]min, \[Xi]max] \[Equal]
1, {kampas, \(-0.04\), 0.04}, MaxIterations \[Rule] 30]\)],
"Input"],
Cell[TextData[{
"Note that ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" require ",
Cell[BoxData[
FormBox[
StyleBox["kampas",
"Input"], TraditionalForm]]],
" to be numeric before they will evaluate. This function appears to
work \
fine,"
}], "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(TraditionalForm\`teisingasKampas[\(-1\), 0, 2]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`{kampas \[Rule] 0.292393448767806`}\)], "Output"]
}, Open ]],
Cell["and a plot shows that the required condition is satisfied.",
"Text"],
Cell[BoxData[
\(TraditionalForm\`\(Plot[
Evaluate[
u(\[Xi]) /.
g[\(-1\), kampas /. teisingasKampas[\(-1\), 0, 2], 0,
2]], {\[Xi], 0, 2}, PlotRange \[Rule] All];\)\)], "Input"]
},
FrontEndVersion->"5.0 for Macintosh",
ScreenRectangle->{{4, 1280}, {0, 832}},
WindowSize->{1081, 684},
WindowMargins->{{20, Automatic}, {Automatic, 44}}
]
(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)
--
Paul Abbott Phone: +61 8 9380 2734
School of Physics, M013 Fax: +61 8 9380 1014
The University of Western Australia (CRICOS Provider No 00126G)
35 Stirling Highway
Crawley WA 6009 mailto:paul at physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul