2
$\begingroup$

Edit: Made question a lot easier to read.

Simple example. - Consider $p(x,y,z,t) = e^t(ax^2 + by^2 + cz^2)$ and I seek the values of $a,b,c$ which make $p\approx0$ at the point $(x=r\cos\theta,y=r\sin\theta,m,t)$ for all $t\ge0$ where $0 < \theta \le 2\pi$. A non-trivial solution will be $(a,b,c)=\{(m/r)^2,(m/r)^2,-1\}$. How could I have found this in Mathematica?

Actual problem description:

  • I have 3 real-valued functions $f(x,y,z,t)$, $g(x,y,z,t)$, $h(x,y,z,t)$ which contain some non-zero parameters/coefficients $a,b,c,d,e$ .

  • At carefully selected values of those parameters, the functions are each forced to be approximately zero for all the points parametrized by $(x=r\cos\theta,y=r\sin\theta,m,t)$, where $0 \le t$ with $0 < r \le R$ and $0 < \theta \le 2\pi$. How do I brute-force or solve for the parameters which drive those functions to zero?

  • In this question, by zero I mean some number $\epsilon$ such that $abs(\epsilon) > 10^{−5}$

Sample functions (the actual functions are much longer and complicated):

Find non-zero $a,b,d$ such that at all the parametrized points as defined above, then $f\approx0$ && $g\approx0$ && $\textrm{Norm}[h]\approx0$ where

f = 1/2 a^2 E^(-2 d^2 t) (E^(2 a x)+E^(2 a y)+E^(2 a z)+2 E^(a (y+z)) Cos[d x+a z] Sin[a x+d y]+2 E^(a (x+y)) Cos[a y+d z] Sin[d x+a z]+2 E^(a (x+z)) Cos[a x+d y] Sin[a y+d z]);
g = (a^2+b^2+ab)*Exp[2*(a^2+b^2+(a+b)^2)*t]*(Exp[a*(x-y)+b*(x-z)] + Exp[a*(y-z)+b*(y-x)] + Exp[a*(z-x)+b*(z-y)]);
h = {a E^(-d^2 t) (E^(a z) Cos[a x+d y]+E^(a x) Sin[a y+d z]),a E^(-d^2 t) (E^(a x) Cos[a y+d z]+E^(a y) Sin[d x+a z]),a E^(-d^2 t) (E^(a y) Cos[d x+a z]+E^(a z) Sin[a x+d y])};

What I have tried:

  1. ReplaceAll for the variables $(x,y)=(r\cos\theta,r\sin\theta)$, enforce constraints on $z,\theta,r$ and $t$ then use Minimize. This doesn't work, Mathematica simply returns the input command. I have also asked this on another question. e.g.

    points = {x^2+y^2->r^2,x->r*Cos[\[Theta]],y->r*Sin[\[Theta]]}; Minimize[{f /. points, 0 < z <= 200 && 0 <= \[Theta] < 2\[Pi] && t > 0}, {a,d}]

  2. ReplaceAll for the variables, and then SolveAlways for $z,r,\theta$ and $t$ in the constraints. I also tried Reduce. Neither worked e.g. for SolveAlways,

    fpoints = Simplify[f, {x^2 + y^2 -> r^2, x->r*Cos[theta], y->r*Sin[theta]}]; gpoints = Simplify[g, {x^2 + y^2 -> r^2, x->r*Cos[theta], y->r*Sin[theta]}]; hpoints = Simplify[h, {x^2 + y^2 -> r^2, x->r*Cos[theta], y->r*Sin[theta]}]; SolveAlways[fpoints==0 && gpoints == 0 && hpoints == 0 && t >= 0 && 0 < theta <= 2*Pi && m >= z > 0,{theta,r,z,t}]

  3. Manual brute-force process. I tried fixing some parameters (my random guesses), evaluating the functions and using FindMinValue and FindMaxValue for all theta, t and z e.g. fpoints = Simplify[f, {x^2 + y^2 -> r^2, x->r*Cos[theta], y->r*Sin[theta]}]; FindMinValue[{fpoints /.{r->1,a->1,b->1,c->1/10,d->-1/100,e->25}, {0 <= z < m && 0 <= theta < 2*Pi}}, {theta,t,z} FindMaxValue[{fpoints /.{r->1,a->1,b->1,c->1/10,d->-1/100,e->25}, {0 <= z < m && 0 <= theta < 2*Pi}}, {theta,t,z}

    (For ℎ I did the FindMinValue on Norm[hpoints] after the ReplaceAll). By randomly changing the parameters manually, I can occasionally obtain cases in which FindMinValue and FindMaxValue give me numerical zeros e.g. FindMinValue gives −6.07768∗10−7 and FindMaxValue gives me 5.52429∗10−8. Then I move on to the other functions and see if I get numerical zero for $g$ and $h$ as well at those same parameters. This isn't always working out so far. Since I am changing parameters manually I am sure I am definitely missing the sweet spot in between the parameter values.

$\endgroup$
17
  • 2
    $\begingroup$ A smaller but complete example that shows the issue would be quite helpful. $\endgroup$ Commented Mar 27, 2020 at 21:33
  • $\begingroup$ @DanielLichtblau Thanks! Does the smaller example I included help? $\endgroup$
    – Cogicero
    Commented Mar 27, 2020 at 22:23
  • 1
    $\begingroup$ It was not clear to me that m and r are (apparently) adjustable rather than free parameters. So expanding as a series in theta and equating coefficients is perhaps a viable way forward. $\endgroup$ Commented Mar 30, 2020 at 13:57
  • 1
    $\begingroup$ If there is a solution that makes the expression zero then it should work for any order. $\endgroup$ Commented Mar 30, 2020 at 14:13
  • 1
    $\begingroup$ Also, are you expecting an analytical solution or perhaps it can only be solved numerically? $\endgroup$
    – Chris K
    Commented Mar 30, 2020 at 14:19

3 Answers 3

1
$\begingroup$

I have no idea on your actual problem, but your simple example can be solved by choosing a few arbitrary values of θ:

FullSimplify[
  Solve[0 == E^t (a x^2 + b y^2 + c z^2)
    /. {x -> r Cos[θ], y -> r Sin[θ], z -> m} /. θ -> {0, 1, 2}, {a, b}]]
(* {{a -> -((c m^2)/r^2), b -> -((c m^2)/r^2)}} *)
$\endgroup$
2
  • $\begingroup$ Thanks. Please you post this answer on the other bounty question (which this fully answers): mathematica.stackexchange.com/questions/217238 so I can award it. $\endgroup$
    – Cogicero
    Commented Mar 30, 2020 at 14:56
  • $\begingroup$ Thanks, I didn't see that one. $\endgroup$
    – Chris K
    Commented Mar 30, 2020 at 15:17
1
$\begingroup$

It is best to start with visualizations. This hard to visualize.

There are two types of functions in this question: Exp and Sin and Cos.

For f:

a dominate the Exp so has to be taken arbitrarily from the Exp condition. If a is big, then for large {x,y,z} the conditions can not be met at all. The Sin and Cos terms do not cancel out the exponential divergence at all. For smaller a the are then arbitrarily many potential minima. So the x,y,z-ranges have to be constrained to find minimum.

A meaningful approach are visualizations like:

With[{z = 0}, 
 Manipulate[
  Plot3D[E^(2 a x) + E^(2 a y) + E^(2 a z) + 
    2 E^(a (y + z)) Cos[d x + a z] Sin[a x + d y] + 
    2 E^(a (x + y)) Cos[a y + d z] Sin[d x + a z] + 
    2 E^(a (x + z)) Cos[a x + d y] Sin[a y + d z], {x, -6, 6}, {y, -6,
     6}], {a, 0, 2}, {b, 0, 2}, {d, 0, 2}]]

If did put in sliders for a,b and d just set the other coordinates zero.

This is a z=0 cut through the four-dimensional function.

It did put t=0 first since this is only an amplitude to the combination of trigonometrical function with exponential coefficients.

For most triples of a,b and d the function f looks like an egg tray. The minima are a zig-zag. The maxima grow like the exponential functions suggest.

The use of Series is a nice help. To zeros order this results in a nice complicated polynomial in a,b and d.

Since everything does not really work of the shelve from Mathematica good old knowledge must be applied. The only method are Lagrange multipliers as additional dimensions. One for each f,g and h and for the minimization constrains: $f\approx0$ && $g\approx0$ && $\textrm{Norm}[h]\approx0$.

Simpler version in Mathematica are presented in many questions here on Mathematica.stackexchange. These are easily generalized and coded for this case. But take under consideration the is five-dimensional from the definition of the function with Exp and trigonometrics and there are there conditions to be met.

Take for example this: How can I implement the method of Lagrange multipliers to find constrained extrema?

This question goes a bit further in visualization: 4D Table, Slider limitation and plotting ListSliceContourPlot3D

It is hard to decide whether the solutions will be maxima or minima without a visual idea. The only point is the this will be extrema in the complete space. Since no values for a, b or d are given numerical simulation will not be easily find something.

The conditions are only given in a proximal form and not as equation is the main flaw.

Hope that solves the question so far.

$\endgroup$
1
+50
$\begingroup$

Firstly, one needs to appreciate that there is no unique solution to the problem at hand, so one cannot ask Mathematica to find the expected answer. At best, it can find the relation between the dependent variables ($a,b,c$ in the simplified example) in terms of independent variables ($r,t,\theta,m$). Fortunately, this is quite doable.

The main function that we will be using is Reduce, which will give all possible cases for which given function is zero. We will then eliminate the possibilities which require independent variables to take specific values (such as $r=0$), as these cases are isolated solutions whereas we are interested in generalized solution. Finally, we will assume that the required inequalities for the solution is satisfied (such as $r\ne0$, which is actually a condition for OP's preferred answer as it involves $a=m^2/r^2$).

The code to do these is as follows:

ClearAll[solve];
Options[solve] = Options[FullSimplify];
solve[independentParameters_List, opts : OptionsPattern[]] := 
  Module[{condition, replace},
   condition[a_, b_] := Table[FreeQ[Equal[a, b], i], {i, 
   Subsets[Alternatives @@ independentParameters, {Length[independentParameters] - 1}]}];
   replace = Equal[a_, b_] :> False /; (Or @@ condition[a, b]);
   FullSimplify[Reduce[# == 0] /. replace /. Unequal[a_, b_] -> True,
   Assumptions -> opts]
  ] &;

We can see it in action as follows:

p[a_, b_, c_][t_, x_, y_, z_] := Exp[t] (a x^2 + b y^2 + c z^2);
p[a, b, c][t, r Sin[\[Theta]], r Cos[\[Theta]], m] // solve[{r, \[Theta], m, t}]
(* a + b Cot[\[Theta]]^2 + (c m^2 Csc[\[Theta]]^2)/r^2 == 0 *)

which gives the general solution:

$$a+b \cot ^2(\theta )+\frac{c m^2 \csc ^2(\theta )}{r^2}=0$$

We can now fix any solution we like; for example, we can get OP's result back as follows:

a + b Cot[\[Theta]]^2 + (c m^2 Csc[\[Theta]]^2)/r^2 == 0 /. {c -> -1, b -> m^2/r^2} // FullSimplify
(* a == m^2/r^2 *)

The code as written is quite general and should work with other input. In particular, we did not make use of anything specific to the function $p(t,x,y,z)$. As an example, consider a similar yet modified function:

p2[a_, b_, c_][t_, x_, y_, z_] := Exp[2 t] (a x^4 + b y^4 + c z^4);
p2[a, b, c][t, r Sin[\[Theta]], r Cos[\[Theta]], m] // solve[{r, \[Theta], m, t}]
(* a + b Cot[\[Theta]]^4 + (c m^4 Csc[\[Theta]]^4)/r^4 == 0 *)

for which fixing $$c=-1\;,b=\frac{m^4}{r^4}$$ gives us the answer $$a=\frac{m^4 \left(\cot ^4(\theta )+\csc ^4(\theta )\right)}{r^4}$$

We can of course use the code for functions with other number of variables. For example:

p3[a_, b_, c_, d_][t_, x_, y_, z_, u_] := Exp[t] (a x^2 + b y^2 + c z^2 + d u^2);
p3[a, b, c, d][t, r Sin[\[Theta]], r Cos[\[Theta]] Sin[\[Phi]], r Cos[\[Theta]] Cos[\[Phi]], m] // solve[{r, \[Theta], \[Phi], m, t}]
(* a + (d m^2 Csc[\[Theta]]^2)/r^2 + Cot[\[Theta]]^2 (c Cos[\[Phi]]^2 + b Sin[\[Phi]]^2) == 0 *)

for which fixing $$b=c=\frac{m^2}{r^2}\;, d=-1$$ fixes $$a=\frac{m^2}{r^2}$$

$\endgroup$
5
  • $\begingroup$ Thanks! This looks great. Hmmm... I've been running it for the past 20 hours on my actual object function, and it's still running, so it's hard to tell if it's working out yet. Is there any way to break this down into parts to ensure it hasn't hung? $\endgroup$
    – Cogicero
    Commented Mar 31, 2020 at 17:22
  • $\begingroup$ ANSWER PART 1: Reduce is extremely slow when it comes to complicated expressions, and should really not be used straightforwardly. Technically, it works and produce results but just takes forever (To see that the bottleneck is Reduce, forget about my code and just run Reduce[f==0] where f is your relevant function). $\endgroup$ Commented Mar 31, 2020 at 17:33
  • $\begingroup$ PART 2: I wrote the answer above for sufficiently small expressions;for larger ones,one should either come up with a problem-specific solutions(for example, in your simple example, the result can be chosen $\theta$ independent due to spherical symmetry and one can impose it in the code) or one should decompose the large expression into smaller ones and apply the code to small ones separately. For the latter, Factor can be handy: you can even write a small code, which factors a given expression into smaller ones, apply my function to each case, and return the results with an overall Or. $\endgroup$ Commented Mar 31, 2020 at 17:34
  • $\begingroup$ You're right. Reduce [f == 0, {a,b,c}] does run forever also (I had previously tried that). You'e obviously put in a lot of thought and I wish I could award this. I'm starting to think a brute-force search across a grid might be the only way out here. Such as this mathematica.stackexchange.com/questions/218391 $\endgroup$
    – Cogicero
    Commented Mar 31, 2020 at 21:44
  • $\begingroup$ Well, thanks for your assistance. I will award the bounty even though this didn't work for me. $\endgroup$
    – Cogicero
    Commented Apr 1, 2020 at 7:34

Your Answer

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