0
$\begingroup$

I'm trying to solve the following system of equations. The first two equations depend on the variable TeV

RxA[TeV_] :=   Exp[E1/TeV]/ngas*(2*Pi*(me*mp/mp)*qe*TeV/hp^2)^(3/2);(*kT in eV*)
RxB[TeV_] :=  Exp[E2/TeV]/ngas*(2*Pi*(mp^2/(2*mp))*qe*TeV/hp^2)^(3/2);(*kT in eV*)
eqpart = NSolve[{2*y^2/(x - y) == RxA,4*(x - y)^2/(1 - x) == RxB}, {x, y}];
Plot[Evaluate[{x, y} /. eqpart], {TeV, 0.1, 10}, PlotRange -> All]

I don't understand why there is no output on the plot. Is the syntax ok?

Thanks in advance! bye!

$\endgroup$
1
$\begingroup$
Clear["Global`*"]

You need to provide values for all of the parameters. I have arbitrarily set them all to 1

me = mp = qe = hp = E1 = E2 = ngas = 1;

RxA[TeV_] := Exp[E1/TeV]/ngas*(2*Pi*(me*mp/mp)*qe*TeV/hp^2)^(3/2);(*kT in eV*)
RxB[TeV_] := 
 Exp[E2/TeV]/ngas*(2*Pi*(mp^2/(2*mp))*qe*TeV/hp^2)^(3/2);(*kT in eV*)

You must include arguments for RxA and RxB in the equations, and since eqpart uses a numeric technique it should only be called when TeV has a numeric value.

eqpart[TeV_?NumericQ] := 
  NSolve[{2*y^2/(x - y) == RxA[TeV], 4*(x - y)^2/(1 - x) == RxB[TeV]}, {x, y}];

Show[
 Plot[#[[1]] /. eqpart[TeV], {TeV, 0.1, 10},
    PlotPoints -> 50,
    PlotStyle -> ColorData[97][#[[2]]]] & /@
  {{x, 1}, {y, 2}},
 Frame -> True,
 Axes -> False,
 FrameLabel -> (Style[#, 14, Bold] & /@ {"TeV", "x, y"}),
 PlotRange -> All]

enter image description here

Or, if you want a ParametricPlot

ParametricPlot[{x, y} /. eqpart[TeV],
 {TeV, 0.1, 10}, AspectRatio -> 1]

enter image description here

EDIT: To use Solve rather than NSolve (and again setting all parameters to 1)

eqns = {2*y^2/(x - y) == RxA[TeV],
  4*(x - y)^2/(1 - x) == RxB[TeV]}

(* {(2 y^2)/(x - y) == 2 Sqrt[2] E^(1/TeV) π^(3/2) TeV^(3/2), (4 (x - y)^2)/(
  1 - x) == E^(1/TeV) π^(3/2) TeV^(3/2)} *)

By assuming that x != y && x != 1 the equations can be put in a more manageable form.

eqns2 = eqns // Simplify[#, x != y && x != 1] &

(* {Sqrt[2] E^(1/TeV) π^(3/2) TeV^(3/2) (x - y) == y^2, 
 E^(1/TeV) π^(3/2) TeV^(3/2) (-1 + x) + 4 (x - y)^2 == 0} *)

The solutions are each a ConditionalExpression with the condition that TeV > 0. Simplifying with this assumption eliminates the condition. The solutions are expressed as Root objects.

sol = Solve[eqns2, {x, y}, Reals] //
  Simplify[#, TeV > 0] &

(* {{x -> Root[
    E^(4/TeV) π^6 TeV^6 - 
      64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 - 
         32 Sqrt[2] E^(4/TeV) π^6 TeV^6 + 
         128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 + 
         96 Sqrt[2] E^(4/TeV) π^6 TeV^6 - 
         64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(4/TeV) π^6 TeV^6 - 
         96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
          2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 + 
         32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
          2/TeV) π^3 TeV^3 - 
         512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
          2/TeV) π^3 TeV^3 + 512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 
      256 #1^8 &, 1], 
  y -> Root[-2 E^(
       2/TeV) π^3 TeV^3 Root[
        E^(4/TeV) π^6 TeV^6 - 
          64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 - 
             32 Sqrt[2] E^(4/TeV) π^6 TeV^6 + 
             128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 + 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6 - 
             64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
              4/TeV) π^6 TeV^6 - 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
              2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 + 
             32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
              2/TeV) π^3 TeV^3 - 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
              2/TeV) π^3 TeV^3 + 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 1]^2 + 
      4 E^(2/TeV) π^3 TeV^3 Root[
        E^(4/TeV) π^6 TeV^6 - 
          64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 - 
             32 Sqrt[2] E^(4/TeV) π^6 TeV^6 + 
             128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 + 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6 - 
             64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
              4/TeV) π^6 TeV^6 - 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
              2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 + 
             32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
              2/TeV) π^3 TeV^3 - 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
              2/TeV) π^3 TeV^3 + 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 1] #1 - 
      2 E^(2/TeV) π^3 TeV^3 #1^2 + #1^4 &, 2]}, {x -> 
   Root[E^(4/TeV) π^6 TeV^6 - 
      64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 - 
         32 Sqrt[2] E^(4/TeV) π^6 TeV^6 + 
         128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 + 
         96 Sqrt[2] E^(4/TeV) π^6 TeV^6 - 
         64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(4/TeV) π^6 TeV^6 - 
         96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
          2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 + 
         32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
          2/TeV) π^3 TeV^3 - 
         512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
          2/TeV) π^3 TeV^3 + 512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 
      256 #1^8 &, 2], 
  y -> Root[-2 E^(
       2/TeV) π^3 TeV^3 Root[
        E^(4/TeV) π^6 TeV^6 - 
          64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 - 
             32 Sqrt[2] E^(4/TeV) π^6 TeV^6 + 
             128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 + 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6 - 
             64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
              4/TeV) π^6 TeV^6 - 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
              2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 + 
             32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
              2/TeV) π^3 TeV^3 - 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
              2/TeV) π^3 TeV^3 + 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 2]^2 + 
      4 E^(2/TeV) π^3 TeV^3 Root[
        E^(4/TeV) π^6 TeV^6 - 
          64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 - 
             32 Sqrt[2] E^(4/TeV) π^6 TeV^6 + 
             128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 + 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6 - 
             64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
              4/TeV) π^6 TeV^6 - 
             96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
              2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 + 

             32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
              2/TeV) π^3 TeV^3 - 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
              2/TeV) π^3 TeV^3 + 
             512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 2] #1 - 
      2 E^(2/TeV) π^3 TeV^3 #1^2 + #1^4 &, 2]}} *)

The Plots for the two solution sets are

Column[
 Plot[Evaluate[{x, y} /. sol[[#[[1]]]]], {TeV, 0.1, 10},
    PlotRange -> #[[2]],
    Frame -> True,
    FrameLabel -> (Style[#, 14, Bold] & /@ {"TeV", ""}),
    ImageSize -> 360,
    PlotLegends -> Placed[{x, y}, {.5, .25}]] & /@
  {{1, Automatic}, {2, All}}]

enter image description here

Your actual parameter values and the context of your problem would presumably tell you which solution set(s) is(are) of interest.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much for your help! I confirm it works well now. So my problem was mostly linked to the correct syntax to use (as I suspected). If I move from NSolve to Solve, how should I modify your code? Thanks again! $\endgroup$ – Riccardo Oct 10 '19 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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