# How to find an instance of this inequality?

I am trying to show that a given inequality is true for some triple but I am having trouble finding anything. Here is a tiny version of what I am working with: We define the following polynomials

T[1][x_, y_, z_] := x
T[2][x_, y_, z_] := y
T[3][x_, y_, z_] := z
T[4][x_, y_, z_] := x*z - y
T[5][x_, y_, z_] := y*z - x
T[6][x_, y_, z_] := x^2*z - x*y - z
T[7][x_, y_, z_] := x*y*z - x^2 - y^2 + 2
T[8][x_, y_, z_] := y^2*z - x*y - z
Test[x_, y_, z_] := x*y*z^2 - x^2*z - y^2*z + z
Rad[x_] = N[Abs[x + Sqrt[x^2 - 4]]/2]


Now I define a lists the functions T1 to T8 altered a little bit, mainly by throwing them into Rad and taking an appropriate root.

f1[x_, y_, z_] = Table[Rad[T[i ][x, y, z]], {i, 1, 2}];
f2[x_, y_, z_] = Table[(Rad[T[i ][x, y, z]])^(1/2), {i, 3, 3}];
f3[x_, y_, z_] = Table[(Rad[T[i ][x, y, z]])^(1/3), {i, 4, 5}];
f4[x_, y_, z_] = Table[(Rad[T[i ][x, y, z]])^(1/4), {i, 6, 8}];


So now $$f1$$ is just $$T1,T2$$ composed with $$Rad, f2$$ is $$T3$$ composed with Rad and then take the square root. $$f3$$ is $$T4,T5$$ composed with Rad and then you take the $$(1/3)$$ power. etc. In reality my list of polynomials is longer but they get modified in a similar manner.

Now I am looking for m,k,l where the following inequality holds:

Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l], ] < Rad[Test[m,k,l]]^(1/6)


This is equivalent to 8 inequalities of polynomials with various roots involved. I thought it would not be too hard for Mathematica to find an example but it is not handling it at least FindInstance/NSolve cannot do it. I played around with it and Mathematica clearly is having trouble with the square roots. I ran:

NSolve[Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l],
f5[m, k, l], f6[m, k, l]] < (Rad[Test[m, k, l]])^(1/6), {m, k, l}]


I left it running overnight and it could not find anything

Another idea I had was to look at a big mesh:

Do[If[
Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l],
f5[m, k, l],
f6[m, k,
l]] < (Abs[(Test[m, k, l] + Sqrt[Test[m, k, l]^2 - 4])/2 ])^(1/
6),
Print[{m, k, l}]], {m, 0., 1000., 10.}, {k, 0., 1000., 10.}, {l,
0., 1000., 10.}]]


However, for some reason the mesh takes a really really long time. Not that if found anything. How could i speed up the mesh? I know it can be quicker as the following code runs much faster It just checks random numbers to see if any of them satisfy the inequality:

Do[m = RandomReal[{0, 10000000}]; k = RandomReal[{m, 10000000}];
l = RandomReal[{0, 10000000}]; If[
Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l],  < (Rad[Test[m, k, l]])^(1/6),
Print[{m, k, l}]], 50000]]


MY question is, how do I speed FindInstance (if possible?) Why is my mesh so slow? I have a feeling everything comes down to square roots and how Methematica is trying to use exact algebra rather than approximations, but I am not sure. I tried changing WorkingPrecision or using N command to no avail. Any help/hints would be greatly appreciated.

EDIT: Thank you flinty for the suggestion. Changed the definitions of $$f1,...,f4$$ using your advice

• You could make your life a bit easier if you indexed T like this for example: T[5][x_,y_,z_]:=... instead of T5[x_,y_,z_]:=... and then you wouldn't need all of those ToExpression[StringJoin["T", ToString[i]] Oct 24, 2021 at 16:54
• @flinty I did not think of that... Thank you! I am very new to Mathematica. You should have seen my first attempts at coding this :P Oct 24, 2021 at 16:55
• @user64494 It is not a polynomial inequality per se, it has some square roots too. Oct 24, 2021 at 17:18

Somewhat changing your code to (?NumericQ doesn't hurt)

T[1][x_?NumericQ, y_?NumericQ, z_?NumericQ] := x;
T[2][x_?NumericQ, y_?NumericQ, z_?NumericQ] := y;
T[3][x_?NumericQ, y_?NumericQ, z_?NumericQ] := z;
T[4][x_?NumericQ, y_?NumericQ, z_?NumericQ] := x*z - y;
T[5][x_?NumericQ, y_?NumericQ, z_?NumericQ] := y*z - x;
T[6][x_?NumericQ, y_?NumericQ, z_?NumericQ] := x^2*z - x*y - z;
T[7][x_?NumericQ, y_?NumericQ, z_?NumericQ] := x*y*z - x^2 - y^2 + 2;
T[8][x_?NumericQ, y_?NumericQ, z_?NumericQ] := y^2*z - x*y - z;
Test[x_?NumericQ, y_?NumericQ, z_?NumericQ] :=  x*y*z^2 - x^2*z - y^2*z + z;
Rad[x_?NumericQ] = Abs[x + Sqrt[x^2 - 4]]/2;
f1[x_, y_, z_] = Table[Rad[T[i ][x, y, z]], {i, 1, 2}];
f2[x_, y_, z_] = Table[(Rad[T[i ][x, y, z]])^(1/2), {i, 3, 3}];
f3[x_, y_, z_] = Table[(Rad[T[i ][x, y, z]])^(1/3), {i, 4, 5}];
f4[x_, y_, z_] = Table[(Rad[T[i ][x, y, z]])^(1/4), {i, 6, 8}];


and making use of NMaximize, I obtain (superfluous , in Max is deleted)

NMaximize[-Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l]] +
Rad[Test[m, k, l]]^(1/6), {m, k, l}]


{9.78323*10^101, {m -> 2.9198*10^101, k -> 2.76565*10^101, l -> -7.2135*10^204}}

It should be noticed the result of

Do[If[Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l],
f5[m, k, l],   f6[m, k,l]] < (Abs[(Test[m, k, l] + Sqrt[Test[m, k, l]^2 - 4])/2])^(1/
6), Print[{m, k, l}]], {m, 0, 1000, 50}, {k, 0, 1000,  50}, {l, -10000, 1000, 50}] // AbsoluteTiming


{236.963, Null}

NMaximize[-Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l]] +  Rad[Test[m, k, l]]^(1/6), {m, k, l}, Method -> "DifferentialEvolution", WorkingPrecision -> 300]


3.74935038801842110353464276516933302348545630212011818458831915811903\ 0793344098122744324893952095295694788455503975127389929059962968704043\ 6881370442253308770250177442044423272044187593686089430206852739192174\ 0595925055750106491697039559913784579335096631505494767041517240115471\ 973857388296811256455, {m -> 1.98013793454019207592242312881283322169808200261541215105513767193\ 0870106777592772985933795485465318119483695895009204189302294808515099\ 9779073309670576359259470504535021376122781762263918698523122418519664\ 8195244174523245204781896173050152749166932705747720831083863681068816\ 520365302792008495405802, k -> 1.9999999999999999999999999999999999999999999999999999999999999\ 9999999999999999999999999999999999999999999999999999999999999999999999\ 9999999999999999999275328811206522664422270216568018771487510297390814\ 4933888201872371628499589283248580142572079374639911411784811573519263\ 4993954318339160502713223156, l -> -52.\ 9652209312619898091777901435838542259450834506662795135424544498398647\ 0551375776535710779069976335053268620868482857525272051942906312477045\ 4441665248094428463493580799352815503503749245996384759173664636995545\ 2753587315953377756535466737332160322349594182435936783874876095419393\ 835831676486330137}}

NMaximize[-Max[f1[m,k,l],f2[m,k,l],f3[m,k,l],f4[m,k,l]]+Rad[Test[m,k,l]]^(1/6),{m,k,l},Method->"NelderMead",WorkingPrecision->100,MaxIterations->1000]


{0, {m -> \ -0.8791191610602191314205455796582100447267293930053710937500000000000\ 000000000000000000000000000000000, k -> -0.270441708311808068441244494924898828995341267717094183822178\ 3554940992632775892932476234742236225088, l -> -0.130332026868731814972253192545810358650959045014917095028118\ 4911145868479857313294263504614289782781}}

RegionPlot3D[-Max[f1[m, k, l], f2[m, k, l], f3[m, k, l],
f4[m, k, l]] + Rad[Test[m, k, l]]^(1/6) >= 0, {m, -10^30,10^30}, {k, -10^30, 10^30}, {l, -10^30, 10^30}, WorkingPrecision -> 35, PlotPoints -> 50]


RegionPlot3D[-Max[f1[m, k, l], f2[m, k, l], f3[m, k, l],
f4[m, k, l]] + Rad[Test[m, k, l]]^(1/6) > 0, {m, -10^30,10^30}, {k, -10^30, 10^30}, {l, -10^30, 10^30},  WorkingPrecision -> 100 , PlotPoints -> 50]


produces an empty plot.

RegionPlot3D[-Max[f1[m, k, l], f2[m, k, l], f3[m, k, l],

• Hey! Thank you for your help. I have a few questions. First of all, what is the function of ?NumericQ ? I am also very confused by NMaximize. FindInstance cannot find $(x,y,z)$ where the value is non negative but NMaximize can find where it is very positive? What is the reason for that? Lastly, the plots at the end make very little sense to me. The second plot should be contained in the first, so why is the first one empty? Oct 25, 2021 at 18:52
• @2132123: (i) Indeed, NMaximize[-Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l]] + Rad[Test[m, k, l]]^(1/6), {m, k, l}, Method -> "DifferentialEvolution", WorkingPrecision -> 300] works without ?NumericQ. I prefer to blow on cold water. (ii) FindInstance and NMinimize are quite different commands. Maybe, FindInstance has problems with complex numbers in f3[x_, y_, z_] though x,y,z are real numbers. The use of NMinimize to find an instance of an inequality is a known method. Oct 26, 2021 at 6:38
• (iii) I think RegionPlot3D[-Max[f1[m, k, l], f2[m, k, l], f3[m, k, l], f4[m, k, l]] + Rad[Test[m, k, l]]^(1/6) > 0, {m, -10^30,10^30}, {k, -10^30, 10^30}, {l, -10^30, 10^30}, WorkingPrecision -> 100 , PlotPoints -> 50] is empty because of huge ranges. Every command has its limitations. Oct 26, 2021 at 6:39