0
$\begingroup$

Below I am trying to find the steady state of x1 (x2 is a func of x1) in terms of a variable q1. x1stst. Then I wish to compute q1 values giving "trace = 0" in equq1 (equq1 includes x1 steady states as a function of q1). Unfortunately NSolve doesn't work due to the hideous polynomial of equq1

(*Find q1*)
q2 = 46.2857;
q3 = 0.111;
q4 = 5.714285714285714`;
q5 = 0.14285714285714288`;
NSolve[q1 q2  (q3 + x1)/(1 + x1) x2 - (q4 x1)/(q5 + x1) == 0 && 
   1/(1 + x1^2) - x2 == 0, {x1, x2}];
x1stst = NSolve[
   q1 q2 (q3 + x1)/(1 + x1) 1/(1 + x1^2) - (q4 x1)/(q5 + x1) == 0, x1];
equq1 = -1 + (5.714285714285714` x1)/(0.14285714285714288` + x1)^2 - 
   5.714285714285714`/(0.14285714285714288` + x1) - (
   46.285714285714285` (0.111` + x1) x2 q1)/(1 + x1)^2 + (
   46.285714285714285` x2 q1)/(1 + x1);
NSolve[equq1 == 0 /. 
   x1stst[[1]] /. {x2 -> 1/(1 + x1stst[[1]]^2)}, q1, Reals]

I tried Together Collect etc. non of them worked. Looked similar questions but again couldn't find anything.

Any advice would be appreciated. Thank you.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ x1stst[[1]] is a replacement rule. An expression like: x2 -> 1/(1 + x1stst[[1]]^2) does not make sense. $\endgroup$ Commented Feb 11, 2021 at 13:20

2 Answers 2

Reset to default
2
$\begingroup$ 
Clear["Global`*"]

(*Find q1*)
q2 = 46.2857 // Rationalize;
q3 = 0.111 // Rationalize;
q4 = 5.714285714285714` // Rationalize;
q5 = 0.14285714285714288` // Rationalize;

Restrict x1 to Reals

x1stst = Solve[q1 q2 (q3 + x1)/(1 + x1) 1/(1 + x1^2) - (q4 x1)/(q5 + x1) == 0,
     x1, Reals] // Simplify;

equq1 = -1 + (5.714285714285714` x1)/(0.14285714285714288` + x1)^2 - 
    5.714285714285714`/(0.14285714285714288` + 
       x1) - (46.285714285714285` (0.111` + x1) x2 q1)/(1 + 
        x1)^2 + (46.285714285714285` x2 q1)/(1 + x1) // Rationalize // 
  FullSimplify

(* -1 - 40/(1 + 7 x1)^2 + (10287 q1 x2)/(250 (1 + x1)^2) *)

x1stst is a rule an cannot be used in an expression such as {x2 -> 1/(1 + x1stst[[1]]^2)}

f[q1_] = (equq1 /. x1stst[[1, 1]]) /. {x2 -> 
      1/(1 + (x1 /. x1stst[[1, 1]])^2)} // N[#, 20] &;

(sol = NSolve[f[q1] == 0, q1, Reals, WorkingPrecision -> 20]) // N

(* {{q1 -> 0.763707}} *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you! is these claims true:1-"Rationalize" makes the code faster due to memory issues 2- from the beginning one needs to restrict variables to Reals-not just the end solution-. x1stst[[1, 1]] is great, didn't know that. $\endgroup$ Commented Feb 11, 2021 at 14:46
  • $\begingroup$ In ...N[#, 20] &; why u need the &? When you write N[#, 20]doesn't guarantee WorkingPrecision -> 20 if not why? $\endgroup$ Commented Feb 11, 2021 at 14:49
  • 1
    $\begingroup$ (1) Rationalize generally makes the code slower; however, it may enable simplification which improves precision and any simplification would help the speed. (2) It is better to use any known constraints/assumptions throughout. (3) The & is needed to close a pure function using Slot (i.e., #). N[#, 20] is meaningless without the closing &. (4) N[#, 20]& will only guarantee WorkingPrecision -> 20 if all parts of the input have precision of at least 20. For example, compare Precision[N[0.5, 20]] with Precision[N[0.5//Rationalize, 20]] $\endgroup$ Commented Feb 11, 2021 at 15:22
1
$\begingroup$

If you assume x2 -> 1/(1 + x1^2) ContourPlot shows possible solutions

ContourPlot[Evaluate[0 == q1 q2 (q3 + x1)/(1 + x1) x2 - (q4 x1)/(q5 + x1) /.x2 -> 1/(1 + x1^2)], {x1, -10, 10}, {q1, -5, 5},FrameLabel -> {x1, q1}]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ That is so great :D thank you! $\endgroup$ Commented Feb 11, 2021 at 14:54
  • $\begingroup$ @confused You're welcome, hope it helps! $\endgroup$ Commented Feb 11, 2021 at 14:59

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.