## Real solutions of third and fourth degree equations

A few hours ago I "discovered" that if a third or fourth degree equation has distinct real solutions, it’s possible to calculate them avoiding complex numbers.

In particular, we have:

poly = (x - 1)(x - 2)(x - 3); {c, b, a} = CoefficientList[poly, x][[1 ;; 3]]; d = a^2/3 - b; e = 2 a^3/27 - a b/3 + c; f = ArcCos[-3 Sqrt[3] e/(2 d Sqrt[d])]; N[{-a/3 + 2 Sqrt[3 d] Cos[(f - 4 Pi)/3]/3,    -a/3 + 2 Sqrt[3 d] Cos[(f - 2 Pi)/3]/3,    -a/3 + 2 Sqrt[3 d] Cos[(f - 0 Pi)/3]/3}] 

{1., 2., 3.}

poly = (x - 1)(x - 2)(x - 3)(x - 4); {d, c, b, a} = CoefficientList[poly, x][[1 ;; 4]]; e = 3 a^2/4 - 2 b; f = 2 c - a b + a^3/4; g = b^2 + 12 d - 3 a c; h = 27 a^2 d - 9 a b c + 2 b^3 - 72 b d + 27 c^2; i = ArcCos[h/(2 g Sqrt[g])]; j = (e + 2 Sqrt[g] Cos[i/3])/3; N[{-a/4 - (Sqrt[j] + Sqrt[e - j + f/Sqrt[j]])/2,    -a/4 - (Sqrt[j] - Sqrt[e - j + f/Sqrt[j]])/2,    -a/4 + (Sqrt[j] - Sqrt[e - j - f/Sqrt[j]])/2,    -a/4 + (Sqrt[j] + Sqrt[e - j - f/Sqrt[j]])/2}] 

{1., 2., 3., 4.}

After a few moments of joy, however, I noticed that, for example, if I write (x - 10^-15) instead of (x - 1), I get respectively 6.66134*10^-16 and 8.88178*10^-16 instead of 10^-15 as the solution.

I intuitively believe that it’s a numerical problem and that the main cause is ArcCos, but I’m not too sure and also I’m not able to establish if something can be done to solve the issue, or if I have to give up and I must necessarily rely on it to the good Newton method.

Thank you!

## Fit data to a model of two-variable differential equations

I am trying to fit some data to the solution of a differential equation. I am using "ParametricNDSolve" to create a model which should fit the data using "FindFit".

Now, the code I propose is very similar to the examples I have seen to solve similar problems. However, the equation of this case is a little more complicated that the examples of Mathematica manuals (because is a function of two variables, and the initials conditions involves two numeric integrals).

The code keeps running for hours, and it does not appear to solve the problems at all:

data = {{21/100, 0.260882276365091}, {6/25, 0.330910580727009}, {27/ 100, 0.271601829283246}, {3/10, 0.294039749043066}, {33/100,  0.373363616981994}, {9/25, 0.467495450916}, {39/100,  0.50972503848512}, {21/50, 0.639300915026114}, {9/20,  0.679672806174314}, {12/25, 0.693446859556429}, {51/100,  0.697043731283207}, {27/50, 0.736147748563448}, {57/100,  0.706580484286456}, {3/5, 0.719128578284264}, {63/100,  0.762276173639189}, {33/50, 0.788753648395851}, {69/100,  0.802107836002146}, {18/25, 0.803992011056852}, {3/4,  0.796349018322968}, {39/50, 0.795845629380594}, {81/100,  0.808805024532776}, {21/25, 0.812210415525446}, {87/100,  0.811677706096654}, {9/10, 0.796614634731033}, {93/100,  0.802573366818874}, {24/25, 0.801851431134633}, {99/100,  0.799762250401537}, {51/50, 0.803234757269179}, {21/20,  0.812346679044962}, {27/25, 0.804974524275207}, {111/100,  0.813054479870107}, {57/50, 0.802494974488922}, {117/100,  0.811124462792484}, {6/5, 0.810259975308694}, {123/100,  0.810910967136978}, {63/50, 0.815891670404585}, {33/25,  0.803843274767398}, {69/50, 0.814778081357444}, {36/25,  0.809686831953762}, {3/2, 0.829532613157218}, {9/5,  0.815351190536115}, {21/10, 0.828143435644695}, {12/5,  0.791269047273677}, {27/10, 0.821321230878138}, {3,  0.81153212440728}, {33/10, 0.821798500231733}, {18/5,  0.819594302125574}, {39/10, 0.828629593518031}};  model = ParametricNDSolveValue[{ D[v[x, t], t] == d1*D[v[x, t], x, x], v[-5, t] == 0, v[10, t] == flux, v[x, 0] == flux*NIntegrate[Cos[y]^n1, {y, ArcTan[-b1*(x + a1)], Pi/2}]/NIntegrate[Cos[y]^n1, {y, -Pi/2, Pi/2}]}, v, {x, -5, 10}, {t, 0, 2}, {d1, flux, n1, b1, a1}]  fit = FindFit[data,  (c*model[d1, flux, n1, b1, a1][x, 2])/((c - 1)*model[d1, flux, n1, b1, a1][x, 2] + 1), {{d1, 0.001}, {flux, 0.57}, {c, 3.2}, {n1, 2}, {b1,1}, {a1, -0.4}}, x] 

I guess that this code is computationally inefficient, but I need help to understand why.

Thank you all very much

Ps: english is obviously not my first lenguage, so I apologize for the grammar mistakes.

## Pass equations into shaders to define graphics – HLSL or other shaders

Is it possible with HLSL (or other popular shader languages) to pass instead of an image, an equation that would define the pixel color / position output by the shader? This would allow for more dynamic drawing of smooth shapes that could be sampled with infinite detail. They would also by able to change and morph in real time.

Not sure how this is possible without passing a string with the eq to the shader and having it parse it. Seems like that might lose all the performance gained by putting it onto the GPU.

## How to minimize a problem with differential equations in it?

My problem is to minimize a target function with differential equations in it. The target function is $$\min \mathrm{J}\left( p \right) =\sum_{i=1}^5{|x_i\left( p \right) -y_i|}\,\,$$ where $$x_i(p)$$ are differential equations like this:

stateFunctionV[v1_, v2_, v3_, v4_, v5_, v6_, v7_, v8_] := { \[Mu] = (67/100)*(x2[t]/(x2[t] + 28/100))*(1 - x3[t]/940)*(1 - x4[t]/1050)^3*(1 - x5[t]/361); x1'[t] == \[Mu]*x1[t], x2'[t] == -(v1 + \[Mu]/v5)* x1[t], x3'[t] == (v2 + \[Mu]*v6)*x1[t], x4'[t] == (v3 + \[Mu]*v7)* x1[t], x5'[t] == (v4 + \[Mu]*v8)* x1[t], x1[0] == 0.2025, x2[0] == 441.337, x3[0] == 0, x4[0] == 0,  x5[0] == 0} 

and $$y_i$$ is the experiment data.

The code:

PIP1[value_] := { Sum[Sum[Abs[ NDSolveValue[{ [Mu] = (67/100)*(x2[t]/(x2[t] + 28/100))*(1 - x3[t]/940)*(1 - x4[t]/1050)^3*(1 - x5[t]/361); x1'[t] == \[Mu]*x1[t], x2'[t] == -(value[[1]] + \[Mu]/value[[5]])* x1[t], x3'[t] == (value[[2]] + \[Mu]*value[[6]])*x1[t], x4'[t] == (value[[3]] + \[Mu]*value[[7]])* x1[t], x5'[t] == (value[[4]] + \[Mu]*value[[8]])* x1[t], x1[0] == 0.2025, x2[0] == 441.337, x3[0] == 0, x4[0] == 0,  x5[0] == 0}, {x1, x2, x3, x4, x5}, {t, 0, 6.5}][[i]] [data2[[i, j, 1]]] - data2[[i, j, 2]] ], {j, 1, 7}], {i, 1, 5}]} NMinimize[PIP1[{v1, v2, v3, v4, v5, v6, v7, v8}], {v1, v2, v3, v4, v5, v6, v7, v8}] 

doesn’t work. I also tried to take the differential functions out as a new function and it also failed.

So I wonder if there’s a way to minimize this problem with NMinimize or not? I would be grateful for your advice on how to couple NMinimize and NDSolveValue.

The experiment data are below:

data2 = {{{0, 0.2025}, {2, 0.76}, {3, 1.28}, {4, 2.2}, {5, 2.75}, {6,       2.95}, {6.5, 3.055}}, {{0, 441.337}, {2, 400.435}, {3,       300.326}, {4, 219.348}, {5, 120.109}, {6, 50}, {6.5, 25}}, {{0,       0}, {2, 42.2895}, {3, 78.6579}, {4, 117.026}, {5, 152.158}, {6,       229.553}, {6.5, 224.1205}}, {{0, 0}, {2, 1.66667}, {3,       13.55}, {4, 16.6667}, {5, 22.7333}, {6, 34.4333}, {6.5,       71.2333}}, {{0, 0}, {2, 6.73913}, {3, 7.86957}, {4, 9.3913}, {5,       8.34783}, {6, 10.9565}, {6.5, 11.8261}}}; 

## Equations with Tensor product and Ket in Mathematica:

I tried to express this equation in Mathematica:

I defined necessary things:

P = {{Ket[0], Ket[2], Ket[1], Ket[3], Ket[5], Ket[4], Ket[6], Ket[8],     Ket[7]}, {Ket[2], Ket[1], Ket[0], Ket[5], Ket[4], Ket[3], Ket[8],     Ket[7], Ket[6]}, {Ket[1], Ket[0], Ket[2], Ket[4], Ket[3], Ket[5],     Ket[7], Ket[6], Ket[8]}, {Ket[6], Ket[8], Ket[7], Ket[0], Ket[2],     Ket[1], Ket[3], Ket[5], Ket[4]}, {Ket[8], Ket[7], Ket[6], Ket[2],     Ket[1], Ket[0], Ket[5], Ket[4], Ket[3]}, {Ket[7], Ket[6], Ket[8],     Ket[1], Ket[0], Ket[2], Ket[4], Ket[3], Ket[5]}, {Ket[a], Ket[c],     Ket[b], Ket[6], Ket[8], Ket[7], Ket[\[Alpha]], Ket[\[Gamma]],     Ket[\[Beta]]}, {Ket[c], Ket[b], Ket[a], Ket[8], Ket[7], Ket[6],     Ket[\[Gamma]], Ket[\[Beta]], Ket[\[Alpha]]}, {Ket[b], Ket[a],     Ket[c], Ket[7], Ket[6], Ket[8], Ket[\[Beta]], Ket[\[Alpha]],     Ket[\[Gamma]]}}  hadamand[\[Omega]_, \[Omega]2_] := {   {1, 1, 1, 1, 1, 1, 1, 1, 1},   {1, \[Omega], \[Omega]2, 1, \[Omega], \[Omega]2,     1, \[Omega], \[Omega]2},   {1, \[Omega]2, \[Omega], 1, \[Omega]2, \[Omega],     1, \[Omega]2, \[Omega]},   {1, 1, 1, \[Omega], \[Omega], \[Omega], \[Omega]2, \[Omega]2, \ \[Omega]2},   {1, 1, 1, \[Omega], \[Omega], \[Omega], \[Omega]2, \[Omega]2, \ \[Omega]2},   {1, \[Omega]2, \[Omega], \[Omega],     1, \[Omega]2, \[Omega]2, \[Omega], 1},   {1, 1, 1, \[Omega]2, \[Omega]2, \[Omega]2, \[Omega], \[Omega], \ \[Omega]},   {1, \[Omega], \[Omega]2, \[Omega]2,     1, \[Omega], \[Omega], \[Omega]2, 1},   {1, \[Omega]2, \[Omega], \[Omega]2, \[Omega], 1, \[Omega],     1, \[Omega]2}   } H = hadamand[\[Omega], \[Omega]^2]  

My attempt is:

A[u_, j_][n_, P_, Had_] :=   1/Sqrt[n] Sum[    TensorProduct[Ket[k], Part[P, k, j]]* Bra[k].Had. Ket[u], {k, 0,      n - 1}]  B[1, 2][9, P, H] 

I obtain:

Where is issues?

## Specifying statistical models / Equations between Random Variables

I’m new to Mathematica and confused about how random variables work.

Say I have a standard normal random variable $$X\sim \mathcal N(0,1)$$ and set $$Y := 2X$$. Then $$\operatorname{Cov}(X,Y) = 2$$. Now if I generate $$N$$-iid. draws for $$X$$, the sample covariance of $$X_1,\dots, X_N$$ and $$Y_1,\dots, Y_N$$ with $$Y_n := 2X_n$$ should be close to $$2$$.

My question is: how can I compute statistics (e.g. sample covariance) in Mathematica? Naively, I tried

x:=RandomVariate[NormalDistribution[], 5000] y:=2*x Covariance[x,y] 

but the result is clearly incorrect (usually not close to $$2$$), because in actuality (I presume) under this model $$Y$$ is independent of $$X$$.

What is the most convenient way of achieving the correct result?

## Numerically solving 3D Maxwell equations with NDEigensystem

I am trying to get electric $$\vec{E}$$ and magnetic $$\vec{B}$$ fields in a cylindrical cavity with a dielectric, as in the following Figure. .

Both the cavity (pink) and dielectric (blue, with dielectric constant \epsilon_r) are cylindrical and share the main axis. The cavity is assumed to be conducting, such as the field at its surface has to be null.

By resorting to the vector potential $$\vec{A}$$, and using the generalized Coulomb Gauge transformation such that:

$$\vec{E}(t,\vec{r})=-\partial_t \vec{A}(t,\vec{r}),$$

$$\vec{B}(t,\vec{r})=\vec{\nabla} \times \vec{A}(t,\vec{r}),$$

$$\vec{\nabla} \cdot \left[ \epsilon_r(\vec{r}) \vec{A}(t,\vec{r}) \right]= 0,$$

I got the system of differential equations that $$\vec{A}$$ has to satisfy: $$$$\vec{\nabla}^2 \vec{A}(t,\vec{r}) – \frac{1 + \epsilon_r(\vec{r})}{c^2}\partial_t^2 \vec{A}(t,\vec{r})=0.$$$$ Here, $$\vec{r} = (x,y,z)$$.

I would like to have the eigenfrequencies and (spatial) eigenfunctions of this last operator. As far as I have understood, NDEigensystem with DirichletCondition is what I have to use; with the first output being the frequencies of the modes and the second output the spatial envelopes of the eigenfunctions.

I tried that, and failed (quite miserably). Not only the eigenfrequencies are imaginary and the spatial envelopes of the eigenfunctions do have an imaginary part$$^{*}$$. If I plug the outcome of NDEigensystem into the differential equations, I find that these equations are not even satisfied. I am sure that it is me failing somewhere, but after long time spent trying, I am starting being really frustrated. My code is the following:

{vals, funs} = NDEigensystem[{EqSt[t, x, y, z, z1, z2, e, c] == {0, 0, 0}, BndCnd}, {Ax[t, x, y, z], Ay[t, x, y, z], Az[t, x, y, z]}, t, {x, y, z} \[Element] Cylinder[{{0, 0, 0}, {0, 0, d}}, r], 16, Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}]; 

where EqSt[t, x, y, z, z1, z2, e, c] is the system of differential equations

EqSt[t_, x_, y_, z_, z1_, z2_, er_, c_] := Laplacian[{Ax[t, x, y, z], Ay[t, x, y, z], Az[t, x, y, z]}, {x, y, z}] - (1 + fer[z, z1, z2, er])/c^2 {D[Ax[t, x, y, z], {t, 2}], D[Ay[t, x, y, z], {t, 2}], D[Az[t, x, y, z], {t, 2}]}; 

and BndCnd the boundary conditions

BndCnd = {DirichletCondition[Ax[t, x, y, z] == 0, True], DirichletCondition[Ay[t, x, y, z] == ,True], DirichletCondition[Az[t, x, y, z] == 0, True]}; 

Finally, the dielectric function fer[z, z1, z2, er] that I use to mimic the dielectric is the following:

fer[z_, z1_, z2_, e_] := e (HeavisideTheta[z - z1] - HeavisideTheta[z - z2]); 

which is a step function in the z-th coordinate (= the axis of both cylinders) having value e betweem z1 and z2 and zero elsewhere.

I have tried different methods and different values of "MaxCellMeasure" (there is little improvement, but the error which I obtain by plugging the solution into the differential equations is above one!!!). Do you have any idea what is wrong here?

I have few "bonus questions" here.

-First, I tried increasing the mesh resolution for getting better outcomes, but all the times the eigenfunctions look quite bad, not smooth at all. I guess that a 3D mesh is quite demanding, but I know that the most difficult part is going to be at the dielectric. How can I tell NDEigensystem to have a finer mesh there?

-Second, there is a way to change the normalization of the eigenfunctions? I have read that NDEigensystem spits out eigenfunctions $$\vec{\phi}_i$$ s.t.: $$\int \vec{\phi}_i^* \cdot \vec{\phi}_j d \vec{r} = \delta_i,j$$. I would change that to $$\int \epsilon_r(\vec{r}) \vec{\phi}_i^* \cdot \vec{\phi}_j d \vec{r} = \delta_i,j$$, would this be possible?

-Finally, I guess that it might be easier for NDEigensolver to remove the time dependence from the differential equations. This can be easily done by assuming that $$\vec{A}(t,\vec{r}) \rightarrow \vec{A}(\vec{r})e^{-i \omega t})$$. In this case, one can rewrite the operator equation above as: $$$$\vec{\nabla}^2 \vec{A}(\vec{r}) + \frac{\omega^{2}}{c^2}[1 + \epsilon_r(\vec{r})] \vec{A}(\vec{r})=0.$$$$ However, I do not know how to tell NDEigensystem that $$\omega$$ is not a parameter, but should be found with the system constraints… There is a way to remove the time dependence from the equation to be given in NDEigensystem?

$$^*$$I guess this is a stupid problem; if I change the sign of the Laplacian I get real eigenfrequencies and eigenfunctions… But I am quite confident that the sign in the above equation is correct. There is something that I am missing here?

## Full line trajectories plot for the solution of Second Order nonlinear coupled differential equations

I wanted to plot a phase plane containing the trajectories of the solutions found by using ‘NDSolve’ using the initial conditions for x[0], y[0], x'[0] and y'[0]. The equations are: x”[t] – 2 y'[t] == -x[t] + y[t]^2; y”[t] + 2 x'[t] == x[t] + y[t] + x[t]*y[t]

The equilibrium point for the system is (0,0). I have plotted the stream plot for the system but unable to plot a phase portrait that would give me the full line trajectories of the system for different initial conditions. I am also looking for any periodic solution if present in it. The stream plot I got is given below and I would take initial conditions from it.

I get this by using the Parametric Plot of the NDSolve solution:

Kindly help in this capacity. Thanks in advance.

## System of equations and distinct result sets

Question on the below system of equations:

First off, in section 1 below, I had to manually play with that 4th parameter (300) so I can finally see that there are actually 96 sets of solutions. Is there a way that I do not provide that and Mathematica tells me there are 96? If I remove the parameter I only get 1 set. As I increase it, I eventually see that beyond 96, no matter how high I go it stops at 96.

Second, in section 2 below, I see the CountDistinct gives me the total number of sets. But if you pay attention there are only 4 unique sets (3,21,27,29), (9,11,27,33), (7,13,29,31), (11,13,19,37) – But there are 24 permutations for each set, since the 4 vars are interchangeable, hence 24*4=96 total). How can I ask Mathematica to also list these 4 unique sets?

eqn = FullSimplify[{w + x + y + z == 80, w^2 + x^2 + y^2 + z^2 == 2020}] 
Table[FindInstance[eqn, {w, x, y, z}, Integers, 300] ]  {{w -> 3, x -> 21, y -> 27, z -> 29}, {w -> 3, x -> 21, y -> 29,    z -> 27}, {w -> 3, x -> 27, y -> 21, z -> 29}, {w -> 3, x -> 27,    y -> 29, z -> 21}, {w -> 3, x -> 29, y -> 21, z -> 27}, {w -> 3,    x -> 29, y -> 27, z -> 21}, {w -> 7, x -> 13, y -> 29,    z -> 31}, {w -> 7, x -> 13, y -> 31, z -> 29}, {w -> 7, x -> 29,    y -> 13, z -> 31}, {w -> 7, x -> 29, y -> 31, z -> 13}, {w -> 7,    x -> 31, y -> 13, z -> 29}, {w -> 7, x -> 31, y -> 29,    z -> 13}, {w -> 9, x -> 11, y -> 27, z -> 33}, {w -> 9, x -> 11,    y -> 33, z -> 27}, {w -> 9, x -> 27, y -> 11, z -> 33}, {w -> 9,    x -> 27, y -> 33, z -> 11}, ..... ...(snip)....   {w -> 33, x -> 27, y -> 9, z -> 11}, {w -> 33, x -> 27, y -> 11,    z -> 9}, {w -> 37, x -> 11, y -> 13, z -> 19}, {w -> 37, x -> 11,    y -> 19, z -> 13}, {w -> 37, x -> 13, y -> 11, z -> 19}, {w -> 37,    x -> 13, y -> 19, z -> 11}, {w -> 37, x -> 19, y -> 11,    z -> 13}, {w -> 37, x -> 19, y -> 13, z -> 11}} 
CountDistinct[Table[FindInstance[eqn, {w, x, y, z}, Integers, 300] ]]  {96, 4} 

-Thanks

## Solving trigonometric equations with two variables in fixed range?

I am trying to use ‘Solve’ function to solve two trigonometric equations with two variables a1 and a2 (with range 0,2Pi). I wrote the code in the following:

  Vets={Cos[2a1](I+Cos[2a2])+Sin[2a1]Sin[2a2], (I-Cos[2a2])Sin[2a1]+Cos[2a1]Sin[2a2]};   Solve[Vets[[1]] == 1 && Vets[[2]] == 0 && 0<=a1<=2Pi && 0<=a2<=2Pi, {a1,a2}]  

but it gives errors:

Solve::nsmet: This system cannot be solved with the methods available to Solve.

I looked up the documents for Solve function and it should be no problem with the above code. But it doesn’t do anything and give errors, I don’t understand why?

It will be super cool if anyone could answer this question, thank you very much in advance!