Re: again: roots of transcendental function...
- To: mathgroup at smc.vnet.net
- Subject: [mg65428] Re: again: roots of transcendental function...
- From: "David W. Cantrell" <DWCantrell at sigmaxi.org>
- Date: Fri, 31 Mar 2006 06:09:27 -0500 (EST)
- References: <e0gcvl$i78$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Dule <dule23 at gmx.de> wrote:
> Thank y'all for the very helpful suggestions to my previous post!
> I'm wondering that i didn't find "Union" by myself in the Help
> Browser...;-) The "RootSearch" Package is very useful for my application
> as well!
>
> Anyway, i didn't fix my problem yet! I need to fit this model which
> includes the roots of this transcendental function Cot[x] == x/a -
> a/(4*x) to experimental data. I'm using FindMinimum and that works fine
> with other models. The problem now is that i need to define the function
> which needs to be minimized with the mentioned model. The parameter "a"
> is one of the fitting parameters. In previous cases where the objective
> function contained more steps and numerical operations, I defined a
> module including all steps.
> Doing this with FindRoot (or RootSearch) made me a lot of difficulties.
> For every iteration step in FindMinimum the roots for the predicted
> parameter "a" have to be determined and then used for the calculation of
> the model. Therefore i defined a module in which a list of the roots for
> a predicted "a" should be construct. The values of this list should than
> be used for calculating the model. At present i get error messages
> because there is no numerical value for the parameter "a" given.
> In case this sounds confusing, i've attached the code (in this case with
> FindRoot).
> i hope somebody can give me a hint...
>
> time = Range[0.2, 30, 0.2];
>
> model[a_, Ï?_, t_] := Module[{δ, x, gl},
> δ = x /. Union[Table[FindRoot[Cot[x] == x/a - a/(4*x), {x, i}], {i,
> 1, 50}]];
> gl = Exp[a/2 (1 - (t/
> Ï?)/2)] Sum[(δ[[
> i]] (a Sin[δ[[i]]] + 2 δ[[i]] Cos[
> δ[[i]]]) Exp[-(δ[[i]])^2*(t/Ï?)/a])/((δ[[
> i]])^2 + (a/2)^2 + a), {i, 1, Length[δ]}];
> gl
> ]
>
> hju = Map[model[a, Ï?, t] /. {Pe -> 50, Ï? -> 5} /. {t -> #} &, time]
I won't be answering your question here. Sorry. But I think you must to be
warned about something.
You had mentioned earlier that 0 < a < 200, and so parameter a might be
0.01, for example. But if you use, as you do above,
Union[Table[FindRoot[Cot[x] == x/a - a/(4*x), {x, i}], {i, 1, 50}]]
you will get _only three_ out of the 16 positive roots less than 50.
In[1]:=
a = 0.01; Union[Table[FindRoot[Cot[x] == x/a - a/(4*x), {x, i}], {i, 1, 50}]]
Out[1]=
{{x -> 0.0999584}, {x -> 0.0999584}, {x -> 0.0999584}, {x -> 0.0999584},
{x -> 0.0999584}, {x -> 0.0999584}, {x -> 0.0999584}, {x -> 0.0999584},
{x -> 0.0999584}, {x -> 0.0999584}, {x -> 6.28478}, {x -> 21.9916}}
As I mentioned in your previous thread, it is much nicer to use simple
facts about the nature of the roots. We can get _all_ the roots in this
case very easily using just a slight modification of your code (and Union
then becomes superfluous).
In[2]:=
a = 0.01; Table[FindRoot[Cot[x] == x/a - a/(4*x), {x, i}], {i, 10^-6, 50, Pi}]
Out[2]=
{{x -> 0.0999584}, {x -> 3.14477}, {x -> 6.28478}, {x -> 9.42584},
{x -> 12.5672}, {x -> 15.7086}, {x -> 18.8501}, {x -> 21.9916},
{x -> 25.1331}, {x -> 28.2747}, {x -> 31.4162}, {x -> 34.5578},
{x -> 37.6994}, {x -> 40.8409}, {x -> 43.9825}, {x -> 47.1241}}
Alternatively, one could use RootSearch effectively, if its options are set
correctly. Please see Ted Ersek's recent posting "RootSearch options".
David