How to properly write a code for this system of nonlinear equations involving hyperbolic and trigonometric functions

Hi! So I have the values for T1, T2, what im trying to solve are the values of gamma 1, gamma 2 and alpha0. I was able to write the following code but i believe the outputs are wrong:

T1 = 0.5780;
T2 = 0.1561;
syms y1 y2 alpha0
E1 = -T2/T1 + (y1/y2)*((y1^2+alpha0^2)/(y2^2+alpha0^2))^0.5;
E2 = 2+(2+(alpha0^4/(y2^2*(y2^2+alpha0^2))))*cos(y2)*cosh((y2^2+alpha0^2)^0.5)+((alpha0^2/(y2*(y2^2+alpha0^2)^0.5)))*sin(y2)*sinh((y2^2+alpha0^2)^0.5);
E3 = 2+(2+(alpha0^4/(y1^2*(y1^2+alpha0^2))))*cos(y1)*cosh((y1^2+alpha0^2)^0.5)+((alpha0^2/(y1*(y1^2+alpha0^2)^0.5)))*sin(y1)*sinh((y1^2+alpha0^2)^0.5);
eqs = [E1, E2, E3];
[y1,y2,alpha0] = vpasolve(eqs,[y1,y2,alpha0])

Any help is greatly appreciated. Thank you!

Answer 1

[Walter Roberson]

You are right, the solutions fail back-substitution

T1 = 0.5780;
T2 = 0.1561;
syms y1 y2 alpha0
E1 = -T2/T1 + (y1/y2)*((y1^2+alpha0^2)/(y2^2+alpha0^2))^0.5;
E2 = 2+(2+(alpha0^4/(y2^2*(y2^2+alpha0^2))))*cos(y2)*cosh((y2^2+alpha0^2)^0.5)+((alpha0^2/(y2*(y2^2+alpha0^2)^0.5)))*sin(y2)*sinh((y2^2+alpha0^2)^0.5);
E3 = 2+(2+(alpha0^4/(y1^2*(y1^2+alpha0^2))))*cos(y1)*cosh((y1^2+alpha0^2)^0.5)+((alpha0^2/(y1*(y1^2+alpha0^2)^0.5)))*sin(y1)*sinh((y1^2+alpha0^2)^0.5);
eqs = [E1, E2, E3];
[y1sol,y2sol,alpha0sol] = vpasolve(eqs,[y1,y2,alpha0])

y1sol = 

y2sol = 

alpha0sol = 

string(eqs)

ans = 1×3 string array
    "(y1*((alpha0^2 + y1^2)/(alpha0^2 + y2^2))^(1/2))/y2 - 1561/5780"    "cosh((alpha0^2 + y2^2)^(1/2))*cos(y2)*(alpha0^4/(y2^2*(alpha0^2 + y2^2)) + 2) + (alpha0^2*sinh((alpha0^2 + y2^2)^(1/2))*sin(y2))/(y2*(alpha0^2 + y2^2)^(1/2)) + 2"    "cosh((alpha0^2 + y1^2)^(1/2))*cos(y1)*(alpha0^4/(y1^2*(alpha0^2 + y1^2)) + 2) + (alpha0^2*sinh((alpha0^2 + y1^2)^(1/2))*sin(y1))/(y1*(alpha0^2 + y1^2)^(1/2)) + 2"

simplify(subs(eqs, {y1, y2, alpha0}, {y1sol, y2sol, alpha0sol}))

ans = 

[John Harvie delos Reyes]

Thanks! Any suggestion on how to approach the problem?

[Walter Roberson]

My tests suggest that the equations as coded are inconsistent.

... I have not checked yet to see whether they were coded correctly.

... they look to be coded correctly.

[Walter Roberson]

format long g
T1 = 0.5780;
T2 = 0.1561;
syms y1 y2 alpha0
E1 = -T2/T1 + (y1/y2)*((y1^2+alpha0^2)/(y2^2+alpha0^2))^0.5;
E2 = 2+(2+(alpha0^4/(y2^2*(y2^2+alpha0^2))))*cos(y2)*cosh((y2^2+alpha0^2)^0.5)+((alpha0^2/(y2*(y2^2+alpha0^2)^0.5)))*sin(y2)*sinh((y2^2+alpha0^2)^0.5);
E3 = 2+(2+(alpha0^4/(y1^2*(y1^2+alpha0^2))))*cos(y1)*cosh((y1^2+alpha0^2)^0.5)+((alpha0^2/(y1*(y1^2+alpha0^2)^0.5)))*sin(y1)*sinh((y1^2+alpha0^2)^0.5);
eqs = [E1, E2, E3];

E = simplify(sum(eqs.^2));
Efun = matlabFunction(E, 'vars', {[y1, y2, alpha0]});

Efun([.1,.2,.3])

ans = 
          193.713495338566

[bestpos, fval] = fmincon(Efun, [1 -1 0])

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

bestpos = 1×3

1.8740 -4.6941 0.0000

fval = 
         0.184510551313022

[Y1, Y2, AL] = meshgrid(linspace(-50,50,500));
R = Efun([Y1(:), Y2(:), AL(:)]);
[fval, idx] = min(R(:))

fval = 
         0.156681145608294

idx = 
    61620227

besty1 = Y1(idx)

besty1 = 
         -1.90380761523046

besty2 = Y2(idx)

besty2 = 
         -4.70941883767535

bestalpha0 = AL(idx)

bestalpha0 = 
        -0.701402805611222

[Walter Roberson]

format long g
T1 = 0.5780;
T2 = 0.1561;
syms y1 y2 alpha0
E1 = -T2/T1 + (y1/y2)*((y1^2+alpha0^2)/(y2^2+alpha0^2))^0.5;
E2 = 2+(2+(alpha0^4/(y2^2*(y2^2+alpha0^2))))*cos(y2)*cosh((y2^2+alpha0^2)^0.5)+((alpha0^2/(y2*(y2^2+alpha0^2)^0.5)))*sin(y2)*sinh((y2^2+alpha0^2)^0.5);
E3 = 2+(2+(alpha0^4/(y1^2*(y1^2+alpha0^2))))*cos(y1)*cosh((y1^2+alpha0^2)^0.5)+((alpha0^2/(y1*(y1^2+alpha0^2)^0.5)))*sin(y1)*sinh((y1^2+alpha0^2)^0.5);
eqs = [E1, E2, E3];

E = simplify(sum(eqs.^2));
Efun = matlabFunction(E, 'vars', {[y1, y2, alpha0]});

[Y1, Y2, AL] = meshgrid(linspace(-2,-1,300), linspace(-5,-4,300), linspace(-1,0,300));
R = Efun([Y1(:), Y2(:), AL(:)]);
[fval, idx] = min(R(:))

fval = 
       0.00850646301047734

idx = 
     1538486

besty1 = Y1(idx)

besty1 = 
         -1.90635451505017

besty2 = Y2(idx)

besty2 = 
         -4.71571906354515

bestalpha0 = AL(idx)

bestalpha0 = 
         -0.94314381270903

[bestpos, fval] = fmincon(Efun, [besty1, besty2, bestalpha0])

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

bestpos = 1×3

-1.9341 -4.8704 -3.1463

fval = 
      0.000293263727026955

vpa(subs(eqs, [y1, y2, alpha0], bestpos))

ans = 

If we were truly successful then the back-substitution should give all 0.

[John Harvie delos Reyes]

Thank you, good Sir!

If values obtained via back substitution are not consistent with the original equation even using the "best" possible initial test values, could this mean that the problem itself is wrong? I could try other method since im dealing with hyperbolic and trigonometric functions or simply my codes for the equations are wrong.

[Alex Sha]

How about the results below:

y1: -1.90447748745904
y2: -4.93345791290994
alpha0: -4.02843940412624

Fevl:
-8.32667268468867E-16
-3.32534000335727E-12
1.13686837721616E-12

[John Harvie delos Reyes]

Hello, Alex!

Solution found indeed satisfies the problem mathematically but I have to check for the negative values since alpha corresponds to a physical quantity in structural dynamics and y1 and y2 are eigenvalues which when found negative translates to an unstable stiffness matrix (again related to structural dynamics).

Im trying to solve this problem using fsolve currently. May I know the process you took to arrive to those values?

Thank you very much!

Edit:

apparently, postive values of the solutions you found also satisfies the problem

a0 =

4.0241

y1 =

1.9046

y2 =

4.9332

Check:

E1 =

0

E2 =

-8.5265e-14

E3 =

0

Im suspecting the function is very sensitive and could grow exponentially out of control. In any case, did you use some algorithm to arrive at a local maxima/minima?