3
$\begingroup$

I wish to solve this function for PCOAL.

q == (0.0238849 PCOAL + 0.903548 PCOAL^1.86)

If an algebraic solution is not possible, I'd like to ask Mathamatica for an approximation, PCOAL = something_without_PCOAL / q

$\endgroup$
6
  • $\begingroup$ Have you looked up Solve in the documentation? $\endgroup$
    – Jens
    Sep 5, 2016 at 3:12
  • $\begingroup$ In principle, Solve is the right approach, but it yields 93 solutions. In most cases, all but one of these are complex. So, the question is, which of these is desired? $\endgroup$
    – bbgodfrey
    Sep 5, 2016 at 4:10
  • $\begingroup$ @Jens I have seen you voted to put this question on hold. Well, can be that it is formulated not in the best way, however, even in the current version it seems to be more deep than first impression shows. Note that OP asks for the approximation of the solution, not the solution itself. In fact these kind of equations can be solved analytically in terms of hypergeometric functions, see for instance this paper M.L.Glasser sciencedirect.com/science/article/pii/S0377042700002879 $\endgroup$
    – yarchik
    Sep 5, 2016 at 16:59
  • $\begingroup$ @yarchik There was ample time to clarify the question, I think. After all, I asked a question and didn't get a reply. $\endgroup$
    – Jens
    Sep 5, 2016 at 17:04
  • 1
    $\begingroup$ Thanks, yes Solve doesn't work well for this or more complicated instances, so an approximation might be best. $\endgroup$
    – r_al
    Sep 6, 2016 at 20:58

2 Answers 2

3
$\begingroup$

Let us write p instead of PCOAL for the sake of shortness, and first let us plot your equation at, say, q=0.5:

 Plot[{0.5, (0.0238849*p + 0.903548 p^1.86)}, {p, 0, 1}, 
 AxesLabel -> {Style["p", Italic, 16], Style["q", Italic, 16]}]

Here is the result: enter image description here

You get, thus, a single solution which is roughly close to q^1/2. It can be now obtained numerically using the FindRoot routine:

lst = Table[{q, 
   FindRoot[0.0238849*p + 0.903548 p^1.86 == q, {p, Sqrt[q] // N}][[1,
      2]]}, {q, 0, 1, 0.1}]

(*  {{0., 0.}, {0.1, 0.294456}, {0.2, 0.432039}, {0.3, 0.539891}, {0.4, 
  0.632052}, {0.5, 0.71406}, {0.6, 0.788786}, {0.7, 0.857954}, {0.8, 
  0.922694}, {0.9, 0.983794}, {1., 1.04183}}  *)

This solution may be then fitted to some analytical function entitled "model":

model = a*q^0.6 + b*q;
ff = FindFit[lst, model, {a, b}, q]


(*  {a -> 1.22995, b -> -0.190771}  *)

we can visually check the fit quality as follows:

Show[{
  ListPlot[lst, 
   AxesLabel -> {Style["q", Italic, 16], Style["p", Italic, 16]}],
  Plot[model /. ff, {q, 0, 1}, PlotStyle -> Red]
  }]

yielding

enter image description here

Done. Have fun!

$\endgroup$
1
  • $\begingroup$ This works too: s = p /. First@Solve[q == (0.0238849 p + 0.903548 p^1.86), p]; Plot[s, {q, 0, 1}, PerformanceGoal -> "Speed"]. $\endgroup$
    – bbgodfrey
    Sep 5, 2016 at 15:56
3
$\begingroup$

Let us make your equation little bit more general:

q == a x + b x^p

or

q - a x==b x^p

where in your case

a=0.0238849; 
b=0.903548;
p=1.86;

One can find a progression of approximations for $x$ (PCOAL) in your case by the following recursion:

z[0] = (q/b)^(1/p);
Table[z[i] = (q - b (1 - p) z[i - 1]^p)/(b p z[i - 1]^(p - 1) + a);
 z[i], {i, 3}]

This is kind of the Newton method. In most of the cases the zeroth order approximation is sufficient, i.e. x=z[0]. For more precision you can take x=z[1] which explicitly means x=(p q)/(a + p (b/q)^(1/p) q) . This is already quite accurate.

PCOAL=(1.86 q)/(0.0238849 + 1.76129 q^0.462366)

$\endgroup$
1
  • $\begingroup$ Thank you! sorry I went away for a day, these suggestions seem very responsive and I will check into them now! $\endgroup$
    – r_al
    Sep 6, 2016 at 20:53

Not the answer you're looking for? Browse other questions tagged or ask your own question.