Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

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

Search the Archive

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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg102334] Re: [mg102317] Finding positive / non-complex solutions (newbie)
  • From: Bob Hanlon <hanlonr at cox.net>
  • Date: Thu, 6 Aug 2009 06:29:05 -0400 (EDT)
  • Reply-to: hanlonr at cox.net

The plot first tries to evaluate f symbolically. Restrict definition of f to only numeric input, i.e.,

f[rho_?NumericQ, t_?NumericQ] :=


Bob Hanlon

---- Dan <dan.halligan at gmail.com> 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}]




  • Prev by Date: Strange DSolve behavior for PDE solution - Bug?
  • Next by Date: Re: LinearModelFit regression estimated variance error
  • Previous by thread: Re: Finding positive / non-complex solutions (newbie)
  • Next by thread: Re: Finding positive / non-complex solutions (newbie)