Solve nonlinear systems of trigonometric equations

Question

I face the challenge to solve systems of nonlinear equations including trigonometric functions like this one:

$$ y_1 = x_1 \\ y_2 = x_2 \\ \dot{y}_1 = u_1 \cos(x_3) \\ \dot{y}_2 = u_1 \sin(x_3) \\ \ddot{y}_1 = \dot{u_1} \cos(x_3) - u_1^2 \tan(u_2) \sin(x_3) \\ \ddot{y}_2 = \dot{u_1} \sin(x_3) + u_1^2 \tan(u_2) \cos(x_3) \\ $$

I need to solve this system to $x_{1,2,3},u_{1,2},\dot{u}_1$, so the derivative of $u_1$ is allowed to remain.

I know a solution is

$$ x_1 = y_1 \\ x_2 = y_2 \\ x_3 = \arctan(\dot{y}_2/\dot{y}_1) \\ u_1 = \frac{\dot{y}_1}{\cos(\arctan(\dot{y}_2/\dot{y}_1))} \\ u_2 = \Psi_1(y_1,y_2,\dot{y}_1,\dot{y}_2,\ddot{y}_1,\ddot{y}_2) \\ \dot{u}_1 = \Psi_2(y_1,y_2,\dot{y}_1,\dot{y}_2,\ddot{y}_1,\ddot{y}_2) $$

as shown in this paper. I'm particular interested in the $\Psi_{1,2}$ equations not shown.

How can one solve systems of equations like that using CAS tools like MuPad, Maple and others? I tried using the solve() function of MuPad but didn't get it. Although a by-hand solution might be possible as well, I got similar but larger equations which also have to be solved.

I should add that in my case a approximate numerical solution (if possible) would be also fine (no symbolic notation needed).

Thanks in advance, Clemens

Answer 1

In Maple the solve and eliminate command can be useful.

Below I'll lprint the solutions, for legibility here.

restart;

eq1 := y1=x1:
eq2 := y2=x2:
eq3 := d1y1=u1*cos(x3):
eq4 := d1y2=u1*sin(x3):
eq5 := d2y1=d1u1*cos(x3)-u1^2*tan(u2)*sin(x3):
eq6 := d2y2=d1u1*sin(x3)+u1^2*tan(u2)*cos(x3):

sols := simplify(solve({eq1,eq2,eq3,eq4,eq5,eq6},
                       [x1,x2,x3,u1,u2,d1u1], explicit),
                 size) assuming real:

map(lprint,sols[1]):
  x1 = y1
  x2 = y2
  x3 = arctan(d1y2, d1y1)
  u1 = (d1y1^2+d1y2^2)^(1/2)
  u2 = arctan((d1y1*d2y2-d1y2*d2y1)/(d1y1^2+d1y2^2)^(3/2))
  d1u1 = (d1y1*d2y1+d1y2*d2y2)/(d1y1^2+d1y2^2)^(1/2)

map(lprint,sols[2]):
  x1 = y1
  x2 = y2
  x3 = arctan(-d1y2, -d1y1)
  u1 = -(d1y1^2+d1y2^2)^(1/2)
  u2 = -arctan((d1y1*d2y2-d1y2*d2y1)/(d1y1^2+d1y2^2)^(3/2))
  d1u1 = -(d1y1*d2y1+d1y2*d2y2)/(d1y1^2+d1y2^2)^(1/2)

Esol := eliminate({eq1,eq2,eq3,eq4,eq5,eq6},
                  [x1,x2,x3,u1,u2,d1u1]):

Esol[2];
                           {}

A := map(simplify,simplify(Esol[1]),size) assuming real:

seq(lprint(s=eval(s,A)), s in [x1,x2,x3,u1,u2,d1u1]):
  x1 = y1
  x2 = y2
  x3 = arctan(d1y2/d1y1)
  u1 = signum(d1y1)*(d1y1^2+d1y2^2)^(1/2)
  u2 = arctan((d1y1*d2y2-d1y2*d2y1)*signum(d1y1)/(d1y1^2+d1y2^2)^(3/2))
  d1u1 = (d1y1*d2y1+d1y2*d2y2)*signum(d1y1)/(d1y1^2+d1y2^2)^(1/2)

We can resubstitute and check those answers.

simplify( eval( map( (rhs-lhs),
                     {eq1,eq2,eq3,eq4,eq5,eq6} ),
                A ) ) assuming real;

                          {0}

simplify( eval( map( (rhs-lhs),
                     {eq1,eq2,eq3,eq4,eq5,eq6} ),
                sols[1] ) ) assuming real;

                          {0}

 simplify( eval( map( (rhs-lhs),
                      {eq1,eq2,eq3,eq4,eq5,eq6} ),
                 sols[2] ) ) assuming real;

                          {0}