0
$\begingroup$

I need to solve an polynimial equation, but when i try to use Solve or NSolve in mathematica, its cant find the solution in appropriate time (i interrupted calculation after 5 hour left). The equations is:

Remove["Global`*"];
g = 9.8;
C1 = 3.6;
C2 = 0.25;
C3 = 666;
C4 = 2;
c = -36;
e = 0.86;
P0 = 24.44;
h0 = 0.1;

h = h0*(1 + C3*P^C4);
NSolve[{P/P0 == ((1 - e*(1 - h/e))/(1 - e))^((2 + c)/3)/(1 - h/e)^4.65}, {P}]

Maybe my computer is too weak or I'm using the wrong method?

Improve this question
$\endgroup$
3
  • $\begingroup$ I should check the code before publication, I just tried to move only part that was a problem, I corrected value for x0 I do not quite understand "Also using the dependent variable dP in the RHS of h in the wrong way." - I can express from the first equation the variable P or h and make a single equation, but the result is the same. I apologize if I do not understand some obvious things I have enough experience in dealing with mathematics $\endgroup$
    – user1058051
    Commented Jan 8, 2014 at 4:37
  • $\begingroup$ @Nasser: NDsolve? Unless the post changed, that's not being used. OP: Nasser's points are spot-on - if this is some kind of HW problem, please specify such, and perhaps add some context of the over-arching problem. The MM solvers are not magic machines, you might need to 'help' them. Hint: try plotting your function with P around 0.02368. $\endgroup$
    – ciao
    Commented Jan 8, 2014 at 4:51
  • $\begingroup$ @Nasser Wow, thanks very much, It works. $\endgroup$
    – user1058051
    Commented Jan 8, 2014 at 14:17

1 Answer 1

Reset to default
3
$\begingroup$

Ok, try this to see if works for you. Made a function in P (wish you used lowerLetterCase) and plotted it. For some values of P, the function is complex valued.

One can see the roots at around 0 and .1. Hence use FindRoot to pick them up.

Clear[g, C1, C2, C3, C4, c, e, P0, h0, exp, P]
parms = {g -> 9.8, C1 -> 3.6, C2 -> 0.25, C3 -> 666, C4 -> 2, c -> -36, e -> 0.86, 
        P0 -> 24.44, h0 -> 0.1, exp -> 4.65};
h = h0*(1 + C3*P^C4);
f[P_?NumericQ] := P/P0 - ((1 - e*(1 - h/e))/(1 - e))^((2 + c)/3)/(1 - h/e)^exp;
(*verify*)
Simplify[f[P] /. parms]

Mathematica graphics

Not a nice looking function in P. If you can simply it a little it will help. You can see the numerator on the second term is almost zero. The limit of the second term as P->0 is 0.0039, so may be you can simplify this whole function to a linear function (straight line) in P around zero.

Here is a plot of f[P]

Plot[Re[f[P] /. parms], {P, -.25, .25}, Exclusions -> None]

Mathematica graphics

Plot[Im[f[P] /. parms], {P, -.25, .25}, Exclusions -> None]

Mathematica graphics

use FindRoot 

FindRoot[(f[P] /. parms) == 0, {P, .01}]
(* {P -> 0.0236849} *)

FindRoot[(f[P] /. parms) == 0, {P, .11}]
(* {P -> -0.0256356 + 0.0197882 I} *)

FindRoot[(f[P] /. parms) == 0, {P, -.11}]
(* {P -> 0.0236849 + 1.34501*10^-27 I}  *)

may be someone will see a better way to handle this.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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