Solving 5 trygonometric equations with solve

Trying to solve 5 equations below and it takes ages (5min) for matlab, so perhaps something is wrong. When it finally completes the task it shows:

'Warning: Cannot solve symbolically. Returning a numeric approximation instead.' Please help. I'm not sure what I'm doing wrong. To get numerical values I was calling S.a etc.

gam = 28.024; % gyromagentic ratio, GHz/T
ThetaH = 0.1;
f = 14.27; %freq in GHz
tau = 0.3328; %relaxation in ns
H = 0.604; %Field in T
syms H1 H2 a Heff Theta
eqn1 = H1 == H*cos(ThetaH - Theta) + Heff*cos(Theta)^2;
eqn2 = H2 == H*cos(ThetaH - Theta) + Heff*cos(2*Theta);
eqn3 = f == gam * (H1 * H2)^1/2/(1 + a^2);
eqn4 = sin(2*Theta) == (2*H*sin(ThetaH - Theta))/Heff;
eqn5 = 1/tau == (a*gam*2*pi*(H1 + H2))/(2*(1 + a^2)^2);

eqns = [eqn1 eqn2 eqn3 eqn4 eqn5];
S = solve(eqns,[H1 H2 a Heff Theta]);

Answer 1

[John D'Errico]

You are not doing anything overtly wrong, except perhaps in the assumption that an analytical solution exists at all.

Of course MATLAB just does what it is asked to do, solve the system, or at least try to do so as hard as it can. So it grinds away, creating more and more complicated expressions, substituting one into the other. Could it have complex variable solutions? Of course, it does not care. JUST SOLVE THE THING.

Eventually, it decides to give up, and tries a numerical solver - vpasolve.

In fact, you would have been far better off to recognize in advance that the system surely has no analytical solution. And is there a reason you need super high precision accuracy in the solution anyway, when your parameters are defined with as little as 1 to 4 significant digits themselves?

So arguably, using fsolve would have been a reasonable choice here. Far faster. Or use the slower vpasolve. Skip the grind of telling solve to try to solve it though.

Note that no matter which numerical solver you try to use, fsolve or vpasolve, there will surely be multiple solutions. And the solution you get will depend on the starting values you choose. If you choose no starting values, then vpasolve will make its own arbitrary choice, but that does not mean this is the only solution. In fact, any time you have a complicated system of equations involving trig functions, it is almost certain that if at least one solution exists, then infinitely many solutions may also exist. In this case, I see that theta appears only inside calls to both sin and cos.

Now look at Theta, as it was found.

S.Theta
ans =
131.9546948417568318749204578266

Did it solve the problem?

vpa(subs(eqns,S2),5)
ans =
[ 7.7274 == 7.7274, 7.727 == 7.727, 14.27 == 14.27, 0.015606 == 0.015606, 3.0048 == 3.0048]

It seems to have done so. I only displayed 5 digits of course, but your parameters are known to only 3 digits, so really is there a good reason for more?

Next, I hope you realize, these are RADIANS, NOT degrees. That is a lot of radians. So, let me arbitrarily shoose a different solution.

S2 = S;
S2.Theta = S2.Theta - 2*pi;

vpa(subs(eqns,S2),5)
ans =
[ 7.7274 == 7.7274, 7.727 == 7.727, 14.27 == 14.27, 0.015606 == 0.015606, 3.0048 == 3.0048]

That seems to work just as well.

S2.Theta = rem(S2.Theta,2*pi);
vpa(S2.Theta)
ans =
0.0078033909855158594894357288618223

vpa(subs(eqns,S2),5)
ans =
[ 7.7274 == 7.7274, 7.727 == 7.727, 14.27 == 14.27, 0.015606 == 0.015606, 3.0048 == 3.0048]

As I said, infinitely many solutions. That was considering only Theta.

[Maciej Dabrowski]

Thanks for your answer. Ok, so how to impose constraints on solutions here, e.g. a o to 1, Theta (0 to 90)*pi/180 etc. That should solve the problem maybe?

[Maciej Dabrowski]

And I don't need an exact and high accuracy solution, approximated is fine.

[John D'Errico]

I would not be at all surpised if there are typically 2 independent solutions in the interval [0,2*pi], or viewed in degrees, [0,360].

However, solvers like vpasolve or even fsolve do not permit constraints on the solution. The solver that would most easily permist bounds on the oslution is lsqnonlin, from the optimization toolbox.

And, really, there is no reason at all to use symbolics here. You will want to solve for a zero solution to each equation. So move all terms to one side of the equality.

gam = 28.024; % gyromagentic ratio, GHz/T
ThetaH = 0.1;
f = 14.27; %freq in GHz
tau = 0.3328; %relaxation in ns
H = 0.604; %Field in T

% starting values (I have no idea what these variables mean so a random guess.)
v0 = ones(1,5);
LB = zeros(1,5);
UB = [100 100 100 100 2*pi];

[S,resnorm,resid,exitflgag] = lsqnonlin(@(vars) errfun(vars,gam,ThetaH,f,tau,H),V0,LB,UB);

Here is the objective function for lsqnonlin.

function err = errfun(vars,gam,ThetaH,f,tau,H)
  H1 = vars(1);
  H2 = vars(2);
  a = vars(3);
  Heff = vars(4);
  Theta = vars(5);
  err = [-H1 + H*cos(ThetaH - Theta) + Heff*cos(Theta)^2;
     -H2 + H*cos(ThetaH - Theta) + Heff*cos(2*Theta);
     -f + gam * (H1 * H2)^1/2/(1 + a^2);
     -sin(2*Theta) + (2*H*sin(ThetaH - Theta))/Heff;
     -1/tau + (a*gam*2*pi*(H1 + H2))/(2*(1 + a^2)^2)];

After the call to lsqnonlin, we see:

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

S
S =
       7.7274        7.727       7.5914       7.1264    0.0078185

resnorm
resnorm =
   1.0772e-09
   
resid
resid =
  -9.0835e-07
  -1.7455e-06
  -2.9155e-08
  -3.2761e-05
   9.2773e-08
   
rad2deg(S(5))
ans =
      0.44797

So all equations seem to be satisfied, and it found the effectively same solution from before. Only 0.45 degrees though. So VERY close to 0. However, if I try a different starting point, changing the starting value for Theta, we see this:

v0(5) = 3;

>> [S,resnorm,resid,exitflgag] = lsqnonlin(@(var) errfun(var,gam,ThetaH,f,tau,H),v0,LB,UB)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
S =
       7.7274       7.7269       7.5914       8.3284       3.1338
resnorm =
   2.5153e-23
resid =
   9.9476e-14
   9.9476e-14
  -1.3145e-12
  -7.1228e-15
   4.8379e-12
exitflgag =
     1

This does not lie in the interval you wanted to see however.

rad2deg(S(5))
ans =
       179.55