8
$\begingroup$

Hi I'm very new at Mathematica.

My real problems concerns a much more difficult equation with many parameters, but I want to check that the method is good, so I made up a test equation.

My test equation and roots:

equ = (5 x^3 + z + x z + x^2 (1 + z) == 0)  
sol = x /. Solve[equ, x]

Here is problem. I'm looking to the parameter z such that all the roots of the equation are less than zero (edit: real parts less than 0, where, z is a real number and the roots can be complex:

Reduce[Re[sol[[1]]] < 0 && Re[sol[[2]]] < 0 && Re[sol[[3]]] < 0 && Im[z] == 0, z]

The answer is (I think) z > 4, but the function Reduce doesn't calculate this in a reasonable time. How do I proceed?

EDIT after few days: In equation with more parameters: 

equ=-11.8193 + 12.0482 b - 59.0964 x + 12.0482 a x - 11.8193 d x + 12.0482 b d x + 0.320482 x^2 + 12.0482 c x^2 - 59.0964 d x^2 + 12.0482 a d x^2 + 1. x^3 + 0.320482 d x^3 + 1. d x^4 
Reduce[ForAll[x, equ == 0 \[And] Element[{a, b, c, d}, Reals], Re[x] < 0], {a, b, c, d}]

It is good code? I waiting and waiting for a response from the Mathematica.

Maybe it would be faster if the program could write only one/few numerical answer?

|improve this question
$\endgroup$
  • 3
    $\begingroup$ If the roots can be complex, what do you mean they are less than zero ? $\endgroup$ – b.gates.you.know.what Feb 11 '13 at 13:53
  • $\begingroup$ I mean that Real parts of all roots must be less than 0. $\endgroup$ – user5851 Feb 11 '13 at 14:10
  • $\begingroup$ Just to show how to formulate this with Reduce: Reduce[ForAll[x, 5 x^3 + z + x z + x^2 (1 + z) == 0, x < 0], z, Reals]. This is very fast. If you allow complex z and x, it takes longer than I cared to wait. (When you write these expressions, keep in mind that in Mathematica x<0 implicitly implies that x is real. This is important if you remove Reals from Reduce but keep x<0 instead of Re[x]<0) $\endgroup$ – Szabolcs Feb 11 '13 at 19:30
  • $\begingroup$ Another note: if Reduce seems to hang on inequalities, there's reason to believe that it won't finish in any reasonable time. I remember people commenting that it uses an algorithm called cylindrical algebraic decomposition, which has a "bad" complexity in the number of variables. Please google for details if you're interested, I don't want to give you incorrect information and I am not familiar with this kind of mathematics. $\endgroup$ – Szabolcs Feb 11 '13 at 19:33
  • 1
    $\begingroup$ Following @Szabolcs idea, Reduce[ForAll[x, 5 x^3 + z + x z + x^2 (1 + z) == 0 \[And] Im[z] == 0, Re[x] < 0], z, Complexes] gives the solution immediately. $\endgroup$ – Thies Heidecke Feb 11 '13 at 20:40
9
$\begingroup$

Finding the whole solution set

Inspired by @Szabolcs idea, you can let Reduce solve this problem with the help of existential quantifiers:

Reduce[ ForAll[x, 5 x^3 + z + x z + x^2 (1 + z) == 0 \[And] Im[z] == 0, Re[x] < 0], z]

It immediately affirms your conjecture:

(* z > 4 *)

This problem formulation via ForAll can be read as: Reduce the following statement with respect to z over the complex numbers: "For all x that satisfy the condition that 5 x^3 + z + x z + x^2 (1 + z) == 0 where z is restricted to the real number line, the real part of x is negative."

Finding a numeric solution

If you are not so much interested in all possible solutions but just one that meets certain criteria, you can use NMinimize to help solve that optimization problem. For that to work smoothly it's useful to set up some helper functions first:

PolyRoots[poly_, paramrules : {(_ -> _?NumericQ) ..}] :=
  x /. {ToRules@NRoots[poly /. paramrules, x]}

This takes a polynomial and a list of rules which specify the parameters and their values and returns a list of all roots. The big obfuscated looking pattern match has the purpose of making sure that the function definition fires only when the right hand side of all rules (a -> 123) assumes numeric values.

To make it easier for NMinimize to get an idea how bad the current set of parameter values is, we introduce a function that assigns a penalty value to a set of parameter values based on the sum of the squares of the real parts of all roots if they are positive and zero if they are negative. That allows NMinimze to "sense" it's way to better solutions until it hits the promised land, where solutions have penalty value zero.

Penalty[poly_, paramrules : {(_ -> _?NumericQ) ..}] :=
  Total[HeavisideTheta[#] #^2 & @ Thread @ Re @ PolyRoots[poly, paramrules]]

The final step is setting up the polynomial (including a cosmetic change of variable names so that the variables in our final solution have the desired names) and letting NMinimize do the work:

Module[{poly},
  poly = -11.8193 + 12.0482 b - 59.0964 x + 12.0482 a x - 
    11.8193 d x + 12.0482 b d x + 0.320482 x^2 + 12.0482 c x^2 - 
    59.0964 d x^2 + 12.0482 a d x^2 + 1. x^3 + 0.320482 d x^3 + 
    1. d x^4 /. {a -> aa, b -> bb, c -> cc, d -> dd};
  {penalty, solution} = 
  NMinimize[
    {Penalty[poly, {aa -> a, bb -> b, cc -> c, dd -> d}], a > 0, b > 0, c > 0, d > 0},
    {a, b, c, d}]
 ]

(* {0., {a -> 13.2852, b -> 7.3045, c -> 2.38894, d -> 0.171584}} *)
|improve this answer
$\endgroup$
  • $\begingroup$ You don't even need to specify the domain as Complexes, as that's the default. Nice and +1! $\endgroup$ – Szabolcs Feb 11 '13 at 21:03
  • $\begingroup$ Ah nice, that simplifies the solution even more! $\endgroup$ – Thies Heidecke Feb 11 '13 at 21:10
  • $\begingroup$ This gets a vote. I'm kinda' embarrassed now by the things I tried. $\endgroup$ – Daniel Lichtblau Feb 12 '13 at 3:09
  • $\begingroup$ @Thies Heidecke What do you think about this: equ=-11.8193 + 12.0482 b - 59.0964 x + 12.0482 a x - 11.8193 d x + 12.0482 b d x + 0.320482 x^2 + 12.0482 c x^2 - 59.0964 d x^2 + 12.0482 a d x^2 + 1. x^3 + 0.320482 d x^3 + 1. d x^4 Reduce[ForAll[x, equ == 0 [And] Element[{a, b, c, d}, Reals], Re[x] < 0], {a, b, c, d}] ? $\endgroup$ – user5851 Feb 16 '13 at 14:45
  • 1
    $\begingroup$ Thies Heidecke, I editing your answer to simplify the pattern in the definitions. Also, your chosen line break point (starting a line with :=) is invalid in copy&paste so I corrected it. I hope you don't mind. $\endgroup$ – Mr.Wizard♦ Feb 18 '13 at 10:49
4
$\begingroup$ 
RegionPlot[ And @@ Thread[(Re[sol /. z -> (a + b I)]) < 0], {a, 3, 10}, {b, -10, 10}]

Mathematica graphics

Edit

As per your edit, z is Real (really bad choice for a name), so:

Plot[Re[sol], {z, -1, 10}, Evaluated -> True, PlotRange -> {-2, 2}]

Mathematica graphics

You can see that from z==4 onwards, the two complex conjugate solutions are negative, as well as the real one.

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for the answer, but I think the plot will not work, if I have for example 3 parameters in the equation. $\endgroup$ – user5851 Feb 11 '13 at 14:46
  • 1
    $\begingroup$ @user5851 I think there is no universal answer/method to your problem, as it may end up with almost any set of inequalities in the general case. This is just an example for your toy problem. It may be useful if you try to describe your actual equations. $\endgroup$ – Dr. belisarius Feb 11 '13 at 14:53

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.