symbolic solve not returning answer

Hello I am trying to solve a symbolic equation and get the result in terms of some other symbols but the solve is returning the following output:

Warning: Unable to find explicit solution. For options, see help.

> In sym/solve (line 317)

In matrix_mult (line 38)

ans =

Empty sym: 0-by-1

My code is :

syms x y0 yp0 y1 yp1 real;
%y0 = 27.4071; y1 = 19.9002; yp0 = -0.0009257929; yp1 = 0.000519017;
l = 2.948;
Q29g = 41.1925/l;
Q30g = -33.2563/l;
Q31g = -62.6938/l;
Q32g = 51.9106/l;
 
kQ29 = abs(Q29g)/x;
kQ30 = abs(Q30g)/x;
kQ31 = abs(Q31g)/x;
kQ32 = abs(Q32g)/x;
 
L1 = 9.416;
L2 = 0.541;
 
B6_l = 5;
B6_BL = 5.3371;
B = B6_BL/B6_l;
 
a = 0.030021188;
r0 = 499.6471242;
Q29 = [cosh(sqrt(kQ29)*l), (1/sqrt(kQ29))*sinh(sqrt(kQ29)*l); sqrt(kQ29)*sinh(sqrt(kQ29)*l), cosh(sqrt(kQ29)*l)];
Q30 = [cos(sqrt(kQ30)*l), (1/sqrt(kQ30))*sin(sqrt(kQ30)*l); -sqrt(kQ30)*sin(sqrt(kQ30)*l), cos(sqrt(kQ30)*l)];
Q31 = [cos(sqrt(kQ31)*l), (1/sqrt(kQ31))*sin(sqrt(kQ31)*l); -sqrt(kQ31)*sin(sqrt(kQ31)*l), cos(sqrt(kQ31)*l)];
Q32 = [cosh(sqrt(kQ32)*l), (1/sqrt(kQ32))*sinh(sqrt(kQ32)*l); sqrt(kQ32)*sinh(sqrt(kQ32)*l), cosh(sqrt(kQ32)*l)];
drift1 = [1,L1;0,1];
drift2 = [1,L2;0,1];
YY0 = [y0;yp0];
YY1 = [y1;yp1];
entryMat = drift2*(Q30*(drift1*(Q29*YY0)));
exitMat = drift1*(Q31*(drift2*(Q32*YY1)));
pretty(entryMat)
yentry = entryMat(1,1);
ypentry = entryMat(2,1);
 
yexit = exitMat(1,1);
ypexit = exitMat(2,1);
 
eq1 = 0 == x - B*(a^2*(r0 + yentry)^2/(a^2*r0 + a^2*yentry + 2*a*r0*ypentry + 2*a*yentry*ypentry + 2*yentry - 2*yexit))
solve(eq1,x,'Real',true)

I would like to have the solution for x in terms of y0, yp0, y1 and yp1. But symbolically it doesnt solve but if I enter the values for y0, yp0, y1 and yp1 it returns a numerical solution. How can i get the solution for x in terms of y0, yp0, y1 and yp1. Please let me know. I am fairly new to MATLAB and will really appreciate some help.

Thank you in advance,

Best Regards,

Dipanwita

Answer 1

[John D'Errico]

Do you know that an analytical solution exists at all, with those parameters as effectively unknowns? In fact, I'll claim quite confidently it does not. Have you looked at the mess it creates? Even if a solurion existed, it would surely be so nasty looking that you could not write it down in a large book anyway.

In fact, it is trivial to write equations far simpler than this that have provably no solution that can be formed. Congratulations! You wrote another. :)

If you provide values for the parameters, it will generate a solution.

Finally, be careful even if you use a numerical solver. If I substitute in the parameters you indicate, then plot the result. I also took out the unnecessary place where you set it to zero, since that is implicit in tools like solve and vpasolve.

% after substitution:
eq1
eq1 =
x + (110915138028460900352*((541*cos((737*(30233/(2680*x))^(1/2))/250)*((4269466172889345*cosh((737*(82385/(5896*x))^(1/2))/250))/4611686018427387904 - (1928603208551655*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/70368744177664))/1000 - cos((737*(30233/(2680*x))^(1/2))/250)*((240998723186793*cosh((737*(82385/(5896*x))^(1/2))/250))/8796093022208 - (4269466172889345*sinh((737*(82385/(5896*x))^(1/2))/250))/(4611686018427387904*(82385/(5896*x))^(1/2)) + (886705459556757*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/3435973836800) + (sin((737*(30233/(2680*x))^(1/2))/250)*((4269466172889345*cosh((737*(82385/(5896*x))^(1/2))/250))/4611686018427387904 - (1928603208551655*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/70368744177664))/(30233/(2680*x))^(1/2) + (541*sin((737*(30233/(2680*x))^(1/2))/250)*(30233/(2680*x))^(1/2)*((240998723186793*cosh((737*(82385/(5896*x))^(1/2))/250))/8796093022208 - (4269466172889345*sinh((737*(82385/(5896*x))^(1/2))/250))/(4611686018427387904*(82385/(5896*x))^(1/2)) + (886705459556757*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/3435973836800))/1000 - 2197471291370957/4398046511104)^2)/(2171181777475614225203200*cos((737*(2993001339026887/(140737488355328*x))^(1/2))/250)*((4787086884452255*cosh((737*(4956422973553657/(281474976710656*x))^(1/2))/250))/9223372036854775808 + (99501*sinh((737*(4956422973553657/(281474976710656*x))^(1/2))/250)*(4956422973553657/(281474976710656*x))^(1/2))/5000) - 230688210477147289600000*cos((737*(30233/(2680*x))^(1/2))/250)*((240998723186793*cosh((737*(82385/(5896*x))^(1/2))/250))/8796093022208 - (4269466172889345*sinh((737*(82385/(5896*x))^(1/2))/250))/(4611686018427387904*(82385/(5896*x))^(1/2)) + (886705459556757*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/3435973836800) + 3583566893165861404672000*cos((737*(30233/(2680*x))^(1/2))/250)*((4269466172889345*cosh((737*(82385/(5896*x))^(1/2))/250))/4611686018427387904 - (1928603208551655*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/70368744177664) + 230584300921369395200000*cos((737*(2993001339026887/(140737488355328*x))^(1/2))/250)*((6837753132798593*cosh((737*(4956422973553657/(281474976710656*x))^(1/2))/250))/343597383680000 + (4787086884452255*sinh((737*(4956422973553657/(281474976710656*x))^(1/2))/250))/(9223372036854775808*(4956422973553657/(281474976710656*x))^(1/2)) + (53830041*sinh((737*(4956422973553657/(281474976710656*x))^(1/2))/250)*(4956422973553657/(281474976710656*x))^(1/2))/5000000) - 115292150460684697600000*(cos((737*(30233/(2680*x))^(1/2))/250)*((4269466172889345*cosh((737*(82385/(5896*x))^(1/2))/250))/4611686018427387904 - (1928603208551655*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/70368744177664) + sin((737*(30233/(2680*x))^(1/2))/250)*(30233/(2680*x))^(1/2)*((240998723186793*cosh((737*(82385/(5896*x))^(1/2))/250))/8796093022208 - (4269466172889345*sinh((737*(82385/(5896*x))^(1/2))/250))/(4611686018427387904*(82385/(5896*x))^(1/2)) + (886705459556757*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/3435973836800))*((7314504539970061*cos((737*(30233/(2680*x))^(1/2))/250)*((4269466172889345*cosh((737*(82385/(5896*x))^(1/2))/250))/4611686018427387904 - (1928603208551655*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/70368744177664))/225179981368524800 - (8653018309761255*cos((737*(30233/(2680*x))^(1/2))/250)*((240998723186793*cosh((737*(82385/(5896*x))^(1/2))/250))/8796093022208 - (4269466172889345*sinh((737*(82385/(5896*x))^(1/2))/250))/(4611686018427387904*(82385/(5896*x))^(1/2)) + (886705459556757*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/3435973836800))/144115188075855872 + (8653018309761255*sin((737*(30233/(2680*x))^(1/2))/250)*((4269466172889345*cosh((737*(82385/(5896*x))^(1/2))/250))/4611686018427387904 - (1928603208551655*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/70368744177664))/(144115188075855872*(30233/(2680*x))^(1/2)) + (7314504539970061*sin((737*(30233/(2680*x))^(1/2))/250)*(30233/(2680*x))^(1/2)*((240998723186793*cosh((737*(82385/(5896*x))^(1/2))/250))/8796093022208 - (4269466172889345*sinh((737*(82385/(5896*x))^(1/2))/250))/(4611686018427387904*(82385/(5896*x))^(1/2)) + (886705459556757*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/3435973836800))/225179981368524800) + (230688210477147289600000*sin((737*(30233/(2680*x))^(1/2))/250)*((4269466172889345*cosh((737*(82385/(5896*x))^(1/2))/250))/4611686018427387904 - (1928603208551655*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/70368744177664))/(30233/(2680*x))^(1/2) + (230584300921369395200000*sin((737*(2993001339026887/(140737488355328*x))^(1/2))/250)*((4787086884452255*cosh((737*(4956422973553657/(281474976710656*x))^(1/2))/250))/9223372036854775808 + (99501*sinh((737*(4956422973553657/(281474976710656*x))^(1/2))/250)*(4956422973553657/(281474976710656*x))^(1/2))/5000))/(2993001339026887/(140737488355328*x))^(1/2) - 2171181777475614225203200*sin((737*(2993001339026887/(140737488355328*x))^(1/2))/250)*(2993001339026887/(140737488355328*x))^(1/2)*((6837753132798593*cosh((737*(4956422973553657/(281474976710656*x))^(1/2))/250))/343597383680000 + (4787086884452255*sinh((737*(4956422973553657/(281474976710656*x))^(1/2))/250))/(9223372036854775808*(4956422973553657/(281474976710656*x))^(1/2)) + (53830041*sinh((737*(4956422973553657/(281474976710656*x))^(1/2))/250)*(4956422973553657/(281474976710656*x))^(1/2))/5000000) + 3583566893165861404672000*sin((737*(30233/(2680*x))^(1/2))/250)*(30233/(2680*x))^(1/2)*((240998723186793*cosh((737*(82385/(5896*x))^(1/2))/250))/8796093022208 - (4269466172889345*sinh((737*(82385/(5896*x))^(1/2))/250))/(4611686018427387904*(82385/(5896*x))^(1/2)) + (886705459556757*sinh((737*(82385/(5896*x))^(1/2))/250)*(82385/(5896*x))^(1/2))/3435973836800) - 51918110721334201600000)

x0 = vpasolve(eq1)
x0 =
0.014727457616242243836166601226342

fplot(eq1)
yline(0);
hold on
plot(x0,0,'r*')

The root it found is in red. But those two vertical lines that look like singularities are exactly that. They are not zeros, but places where the function crosses from -inf to +inf. So while there is a root at those locations, it surely will not be the root you want to see.

But is that the real solution you wanted to find?

Lets take that plot, and expand it. Look very closely, as if with a magnifying glass. Zoom in. Now what do you see? Which root did vpasolve find? Look for the red *, as that indicates the root.

As you can see, this function you wish to solve seemes to have almost infinitely many roots.

Honestly, I think you are kidding yourself if you accept anything from this mess as a close approximation to reality.

[Dipanwita Banerjee]

Dear John,

Thanks a lot for your comments. So I did some corrections to the equation and now it looks like this :

p - (2.8841e-04*(y0*(cos(2.9480*(3.3819/p)^(1/2))*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2)) + 0.5410*sinh(2.9480*(4.1890/p)^(1/2))*cos(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2) - 0.5410*sin(2.9480*(3.3819/p)^(1/2))*(3.3819/p)^(1/2)*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2)) + (sinh(2.9480*(4.1890/p)^(1/2))*sin(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2))/(3.3819/p)^0.5000) + yp0*(0.5410*cos(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)) + cos(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000) + (sin(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)))/(3.3819/p)^0.5000 - 0.5410*sin(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000)*(3.3819/p)^(1/2)) + 499.6471)^2)/(2*yp1*(cos(2.9480*(6.3755/p)^(1/2))*(9.4160*cosh(2.9480*(5.2789/p)^(1/2)) + sinh(2.9480*(5.2789/p)^(1/2))/(5.2789/p)^0.5000) + 0.5410*cosh(2.9480*(5.2789/p)^(1/2))*cos(2.9480*(6.3755/p)^(1/2)) - 0.5410*sin(2.9480*(6.3755/p)^(1/2))*(6.3755/p)^(1/2)*(9.4160*cosh(2.9480*(5.2789/p)^(1/2)) + sinh(2.9480*(5.2789/p)^(1/2))/(5.2789/p)^0.5000) + (cosh(2.9480*(5.2789/p)^(1/2))*sin(2.9480*(6.3755/p)^(1/2)))/(6.3755/p)^0.5000) + 2.0009*y0*(cos(2.9480*(3.3819/p)^(1/2))*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2)) + 0.5410*sinh(2.9480*(4.1890/p)^(1/2))*cos(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2) - 0.5410*sin(2.9480*(3.3819/p)^(1/2))*(3.3819/p)^(1/2)*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2)) + (sinh(2.9480*(4.1890/p)^(1/2))*sin(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2))/(3.3819/p)^0.5000) + 30.0000*yp0*(cos(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)) - sin(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000)*(3.3819/p)^(1/2)) + (y0*(cos(2.9480*(3.3819/p)^(1/2))*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2)) + 0.5410*sinh(2.9480*(4.1890/p)^(1/2))*cos(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2) - 0.5410*sin(2.9480*(3.3819/p)^(1/2))*(3.3819/p)^(1/2)*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2)) + (sinh(2.9480*(4.1890/p)^(1/2))*sin(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2))/(3.3819/p)^0.5000) + yp0*(0.5410*cos(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)) + cos(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000) + (sin(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)))/(3.3819/p)^0.5000 - 0.5410*sin(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000)*(3.3819/p)^(1/2)))*(0.0600*yp0*(cos(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)) - sin(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000)*(3.3819/p)^(1/2)) + 0.0600*y0*(sinh(2.9480*(4.1890/p)^(1/2))*cos(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2) - sin(2.9480*(3.3819/p)^(1/2))*(3.3819/p)^(1/2)*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2)))) + 2.0009*yp0*(0.5410*cos(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)) + cos(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000) + (sin(2.9480*(3.3819/p)^(1/2))*cosh(2.9480*(4.1890/p)^(1/2)))/(3.3819/p)^0.5000 - 0.5410*sin(2.9480*(3.3819/p)^(1/2))*(9.4160*cosh(2.9480*(4.1890/p)^(1/2)) + sinh(2.9480*(4.1890/p)^(1/2))/(4.1890/p)^0.5000)*(3.3819/p)^(1/2)) - 2*y1*(cos(2.9480*(6.3755/p)^(1/2))*(cosh(2.9480*(5.2789/p)^(1/2)) + 9.4160*sinh(2.9480*(5.2789/p)^(1/2))*(5.2789/p)^(1/2)) + 0.5410*sinh(2.9480*(5.2789/p)^(1/2))*cos(2.9480*(6.3755/p)^(1/2))*(5.2789/p)^(1/2) - 0.5410*sin(2.9480*(6.3755/p)^(1/2))*(cosh(2.9480*(5.2789/p)^(1/2)) + 9.4160*sinh(2.9480*(5.2789/p)^(1/2))*(5.2789/p)^(1/2))*(6.3755/p)^(1/2) + (sinh(2.9480*(5.2789/p)^(1/2))*sin(2.9480*(6.3755/p)^(1/2))*(5.2789/p)^(1/2))/(6.3755/p)^0.5000) + 30.0000*y0*(sinh(2.9480*(4.1890/p)^(1/2))*cos(2.9480*(3.3819/p)^(1/2))*(4.1890/p)^(1/2) - sin(2.9480*(3.3819/p)^(1/2))*(3.3819/p)^(1/2)*(cosh(2.9480*(4.1890/p)^(1/2)) + 9.4160*sinh(2.9480*(4.1890/p)^(1/2))*(4.1890/p)^(1/2))) + 0.4503)

its still pretty nasty but basically the equation is a function of p, y0, yp0, y1 and yp1 and I would like to solve for p. If I now replace the values for y0, yp0 , y1 and yp1 as:

y0 = 0.0274071+(14.623*-0.0009257929);
y1 = 0.0199002-(13.623*0.000519017);
yp0 = -0.0009257929;
yp1 = 0.000519017;

and ask for the solution I get the following:

p0 = vpasolve(eq1)
fplot(eq1)
yline(0);
hold on
plot(p0,0,'r*')

p0 =
172.3262

which is the correct expected value.

If I now zoom into the red dot solution it looks like this :

So now I am wondering if I give an assumption for p (e.g > 40 && < 300) why do I not see the solution in terms of y0, y1, yp0 and yp1. These values are always real. Please let me know what you think... Thank you again in advance.

[John D'Errico]

Again, if y0, y1, yp0 and yp1 are left as essentially symbolic unknowns, then no solution will be found. You can want for that to be different, but mathematics rarely gives out X-mas presents. Those parameters being always real is not relevant. Magical solutions don't happen, and now matter what you have changes, this is still one incredibly messy function.

In fact, you should see there are actually MULTIPLE solutions. You chose ONE solution. Yet I see several solutions near 0, just looking at that plot. So why choose the one you did? A solution is a solution.

I know. You still want it to be different. Hey, there are a lot of things in this world I wish were different.

Seriously, if you want to do better than this, then you need to prove that what you have is actually valid. It seems like you made a change, and everything got vastly different. What assurance do we have that this version is now accurate? What argument is there to show this particular solution really is a valid one, and that it tells you something meaningful?

So I would spend some SIGNIFICANT effort validating what you have done. Validating that what you have is consistent with reality. It is far too easy to just generate a number and believe the formula you have written down. After all, mathematics cannot lie, now can it? And computers are NEVER wrong, now are they? Of course then there is always the old phrase - garbage in, garbage out.

Once you have done that, AND you can show that this is really the golden solution, then you MIGHT be able to do something, perhaps in terms of a low order Taylor series expansion in mutliple dimensions. Don't bet the farm on it though.

Answer 2

[Walter Roberson]

I am trying to solve a symbolic equation and get the result in terms of some other symbols

l = 2.948;

So, does that mean that l is 2948/1000 exactly? If so, why didn't you write it that way?

fprintf('%.999g\n', 2.948)

2.947999999999999953814722175593487918376922607421875

Is that what you meant by 2.948, the value 6638305850744111 / 2251799813685248 exactly ? If so, why didn't you write it that way?

In terms of notations in science, l = 2.948 means that l is some value that is not precisely known, other than being somewhere between 29475/10000 (inclusive) and 29485/10000 (exclusive). But does it make sense to ask for an indefinitely precise closed-form solution for p when the value of your constants are not precisely known?

Using solve() with equations that have floating point numbers is almost always a category mistake. If you have floating point numbers, vpasolve() might be appropriate, since it is only expected to find approximate solutions.