Using Solve/FSolve for Multiple Trig Equations

[Jung]

Hello,

I am trying to solve for:
y
r
p

in the following equations where I, J, K are known.

I = siny * sinr + cosy * sinp * cosr
J = -siny * cosr + cosy * sinp * sinr
K = cosy * cosp

The code I have used is:

syms r p y i j k

S = solve(i == sin(y)*sin(r)+cos(y)*sin(p)*cos(r), j == -sin(y)*cos(r)+cos(y)*sin(p)*sin(r), k == cos(y)*cos(p))
S.r
S.p
S.y

but I cannot seem to get a result.

Any help would be appreciated.

Thank you,

[Nasser M. Abbasi]

it is normally very hard to obtain analytical solutions for
trig equations since these are nonlinear and involves inverse functions
with branch cuts as well. So to solve for 'y','p', and 'r'
about, you might want to try a numerical approach.

But I am not an expert in this.

--Nasser

[Bruno Luong]

"Jung" wrote in message <klrf49$bq1$1...@newscl01ah.mathworks.com>...

> Hello,
>
> I am trying to solve for:
> y
> r
> p
>
> in the following equations where I, J, K are known.
>
> I = siny * sinr + cosy * sinp * cosr
> J = -siny * cosr + cosy * sinp * sinr
> K = cosy * cosp
>

Your system looks like solving for three Euler angles knowing one column of the rotation matrix. See formulas before "Quaternion".

http://en.wikipedia.org/wiki/Euler_angles

(-J,I,K)' is the last column of [Z1*X2*Y3], where
theta1 = r;
theta2 = -p;
theta3 = y.

In this case meaning that you want to find the solid rotation, from knowing how end up only one unit vector after rotation.

I'm sorry to tell that you don't have unique solution [I wrote this affirmation somewhere in this newsgroup]. Thus your system is not solvable.

Bruno

[Alan_Weiss]

Your problem interested me. I believe that I have a proof that, at least
among real values for p, r, and y, there is no solution.

Look at the third equation
cosy * cosp = 0

Suppose that cos(y) = 0.
Then the first equation becomes
(+-1)*sin(r) + 0 = 0
So sin(r) = 0.

The second equation is
(+-1) * (+-1) + 0 = 0.
This is impossible.

So in the third equation we have cos(p) = 0. so sin(p) = (+-1)
Suppose sin(p) = 1.
Use the addition formulas for cos(a+b) and sin(a+b).
Then the first and second equations become
cos(r-y) = 0
sin(r-y) = 0
But this is impossible.

So suppose sin(p) = -1.
Then the first and second equations become
-cos(r+y) = 0
-sin(r+y) = 0
This is impossible as well.

So there is no solution to the equations, at least over the reals, and I
think even over complex, too.

Alan Weiss
MATLAB mathematical toolbox documentation

[Bruno Luong]

Alan_Weiss <awe...@mathworks.com> wrote in message
> [ snip]

>
> So there is no solution to the equations, at least over the reals, and I
> think even over complex, too.

There is solution(s) only if I^2 + J^2 + Z^2 = 1. And if there is one solution, there is an infinity of them, see my post.

Bruno

[Bruno Luong]

"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message <klrqr9$k14$1...@newscl01ah.mathworks.com>...

Here is a short code to illustrate that the problem of finding 3D rotation matrix R such that R*u = v has infinity solutions:

% Input vectors of norm 1
u = rand(3,1);
u = u/norm(u)
v = randn(3,1) ;
v = v/norm(v);

% (Random) Rotation solution of (R*u) = v
w = 0.5*(u + v);
q = cross(u,v);
r = randn(); % free parameter, any number will work
q = q + r*w;
q = q / norm(q);
d = dot(u, q);
c = q*d;

u1 = u-c;
u1 = u1/norm(u1);
v1 = v-c;
v1 = v1/norm(v1);
% Rodrigues's formula:
k = cross(u1, v1);
costheta = dot(u1,v1);
R =[ 0 -k(3) k(2);
k(3) 0 -k(1);
-k(2) k(1) 0];
R = costheta*eye(3) + R + k*k'*(1-costheta)/sum(k.^2);

% Check
disp(R)
disp(R*u) % close to v

% Bruno

[Bruno Luong]

> Here is a short code to illustrate that the problem of finding 3D rotation matrix R such that R*u = v has infinity solutions:
>
> % Input vectors of norm 1
> u = rand(3,1);
> u = u/norm(u)
> v = randn(3,1) ;
> v = v/norm(v);

> ...

To be clear, in relation of my post #3, the connection with the original problem is:

u = [0 0 1]';
v = [-J,I,K]';

The solution is given by:
r = theta1;
p = -theta2;
y = theta3.

where theta1, theta2, theta3 are 3 Euler's angles of (3x3) matrix R, when using the convention R = [Z1*X2*Y3].

Bruno

[Bruno Luong]

Here is the code using Euler's angles:

% Fake data
y = rand*2*pi;
r = rand*2*pi;
p = rand*2*pi;
I = sin(y) * sin(r) + cos(y) * sin(p) * cos(r);
J = -sin(y) * cos(r) + cos(y) * sin(p) * sin(r);
K = cos(y) * cos(p);
clear y r p;

% Solve

u = [0 0 1]';
v = [-J,I,K]';

v = v/norm(v);
% (Random) Rotation solution of (R*u) = v
w = 0.5*(u + v);
q = cross(u,v);

f = randn(); % free parameter, any number will work
q = q + f*w;

q = q / norm(q);
d = dot(u, q);
c = q*d;
u1 = u-c;
u1 = u1/norm(u1);
v1 = v-c;
v1 = v1/norm(v1);
% Rodrigues's formula:
k = cross(u1, v1);
costheta = dot(u1,v1);
R =[ 0 -k(3) k(2);
k(3) 0 -k(1);
-k(2) k(1) 0];
R = costheta*eye(3) + R + k*k'*(1-costheta)/sum(k.^2);

r = -atan2(R(1,2),R(2,2));
p = -asin(R(3,2));
y = -atan2(R(3,1),R(3,3));

% Check solution
II = (sin(y) * sin(r) + cos(y) * sin(p) * cos(r));
JJ = (-sin(y) * cos(r) + cos(y) * sin(p) * sin(r));
KK = (cos(y) * cos(p));
norm([I J K]-[II JJ KK]) % should be small

% Bruno