Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Finding positive / non-complex solutions (newbie)

• To: mathgroup at smc.vnet.net
• Subject: [mg102336] Re: [mg102317] Finding positive / non-complex solutions (newbie)
• From: Jaebum Jung <jaebum at wolfram.com>
• Date: Thu, 6 Aug 2009 06:29:27 -0400 (EDT)
• References: <200908050944.FAA18226@smc.vnet.net>

f[rho, 0.1] is evaluated first without plugin any value at rho and this cause res = {}.

One way avoid this evaluation is giving addition arguments in f.  For example,

f[rho_?NumericQ, t_] := 
  Module[{Pu, Pv, pi, d, eqn1, eqn2, res, res2, pi2}, 
   pi = (u Pu + (1 - u) Pv)/(2 s);
   d = (Pu - Pv)/(2 s);
   eqn1 = -t pi (1 - pi) + t u (1 - u) d^2;
   eqn2 = d (rho + u + (2 pi - 1) t - (u - (1 - u)) (1 + t d)) + pi;
   res = NSolve[{eqn1 == 0, eqn2 == 0}, {Pu, Pv}];
   res2 = 
    Cases[{Pu, Pv} /. 
       res, {(a_Real)?Positive, (b_Real)?Positive} :> {a, b}][[1]];
   pi2 = (u res2[[1]] + (1 - u) res2[[2]])/(2 s);
   (1 - pi2)/u];

i.e., f will not be evaluated when rho is not numeric.

- Jaebum




Dan wrote:
> Hi,
>
> This is probably easily solved because I'm very much a newbie.  I have
> two equations (eqn1 and eqn2) with two variables and want to find the
> positive (non-zero, no-complex solutions) for when they both equal 0,
> as a function of two other parameters, then plot the results.  So far
> (after several days work), I've found a kind of solution (it plots
> correctly), but produces an error.  The error, as far as I can tell,
> is produced by the "Cases" statement, which is the function I'm using
> to find the positive solutions.  This will work when the parameters
> (theta and rho) are defined, but will fall over when undefined, which
> seems to be what happens when I call Plot.
>
> Any help on this much appreciated (along with any pointers for where
> my code looks obscure / overly complex -- as I say, I'm not used to
> Mathematica).
>
> Cheers
> Dan
>
> In[1]:=
> f[rho_, t_] := Module [{Pu, Pv, pi, d , eqn1, eqn2, res, res2, pi2},
>    pi = (u Pu + (1 - u) Pv)/(2 s);
>    d = (Pu - Pv)/(2 s);
>    eqn1 = -t pi (1 - pi) + t u (1 - u) d^2;
>    eqn2 = d (rho + u + (2 pi - 1) t - (u - (1 - u)) (1 + t d)) + pi;
>    res = NSolve[{eqn1 == 0, eqn2 == 0}, {Pu, Pv}];
>    res2 =
>     Cases[{Pu, Pv} /.
>        res, {(a_Real)?Positive, (b_Real)?Positive} :> {a, b}][[1]];
>    pi2 = (u res2[[1]] + (1 - u) res2[[2]])/(2 s);
>    (1 - pi2)/u
>    ];
>
> In[2]:= u = 0.004;
> s = 0.02;
>
> In[4]:= f[1, 0.1]
>
> Out[4]= 0.237377
>
> In[5]:= f[10, 0.1]
>
> Out[5]= 0.0081024
>
> In[6]:= LogLogPlot[{f[rho, 0.1], f[rho, 10]}, {rho, 0.01, 100}]
>
>
>
>   

• References:
  • Finding positive / non-complex solutions (newbie)
    • From: Dan <dan.halligan@gmail.com>