imaginary solution with solve

I am trying to solve for 2 variables R2 and theta. Since the systems are not linear, equationsToMatrix and linsolve don't work. So I tried solve.

At one point I was getting a symbolic solution with i in it (which I have since been unable to reproduce) and so I tried to initialize the variables as 'real':

E = 447000;
v = 849968100;
deltay = .00127;
q = 40;
t = .000211;
R = .0525;

R2 = sym('R2','real');
theta = sym('theta','real');
eqns = [deltay*(2*E*t)/(q*(1-v)) == (R2^2)*(1-cos(theta)), sin(theta) == R/R2];
vars = [R2 theta];
[solR2,soltheta] = solve(eqns, vars)
% R2 = double(R2)
% theta = double(theta)

And now, instead of getting an imaginary solution, I get an empty symbolic variable. In the end, I'd like to get a real numeric solution for R2 and theta.

TIA! :)

Answer 1

[Walter Roberson]

There are no real-valued solutions for R2 and theta for those constants.

The system of equations is a quadratic. The two forms are

R2 = 2*(-(1/4)/(((1/4)*R^2*(-1+v)*q+deltay*E*t)*(-1+v)*q))^(1/2)*t*deltay*E

theta = arctan(R/((-(1/4)/(((1/4)*R^2*(-1+v)*q+deltay*E*t)*(-1+v)*q))^(1/2)*t*deltay*E), (-2*deltay*E*t-R^2*(-1+v)*q)/(E*t*deltay))

and

R2 = -2*(-(1/4)/(((1/4)*R^2*(-1+v)*q+deltay*E*t)*(-1+v)*q))^(1/2)*t*deltay*E

theta = arctan(-R/((-(1/4)/(((1/4)*R^2*(-1+v)*q+deltay*E*t)*(-1+v)*q))^(1/2)*t*deltay*E), (-2*deltay*E*t-R^2*(-1+v)*q)/(E*t*deltay))

If you examine the ^(1/2) portions, you will see that they involve -1 times an expression. When all of the constants are positive and v is greater than 1, then the expression is strictly positive, leading to a negative inside the sqrt, which will be imaginary.

At least one of your constants would need to change sign, or your v would need to be less than 1, in order to be taking the square root of a positive value there.

[monarch]

OK that was very helpful because it helped me realize the equations for a few of my constants were wrong! I'm having to plug in a "reasonable" value for v and ignore my equations, but I'm getting closer! I no longer see i in the symbolic solution, but it still won't convert to numeric...

E = 447000;
v = .4;
deltay = .0127;
q = 40;
t = .0211;
R = .0525;

syms R2 theta
eqns = [(deltay*(2*E*t))/(q*(1-v)) == (R2^2)*(1-cos(theta)), sin(theta) == R/R2];
vars = [R2 theta];
[solR2,soltheta] = solve(eqns, vars)
% R2 = double(R2)
% theta = double(theta)

[Walter Roberson]

double(solR2), double(soltheta)

[monarch]

now we're cooking with gas! Thanks!

One final (I think) question:

is there a way to select the solution of the 2 that is for theta >0 && < pi/2 ? I'm not sure how I would use a for loop in this context if that is even the best way to do it...?

[Walter Roberson]

 dR2 = double(solR2)
 dtheta = double(theta);
 mask = dtheta > 0 & dtheta < pi/2;
 dR2(mask), dtheta(mask)

You might also consider

syms R2
syms theta positive

which might eliminate the negative branch solution.

[David Goodmanson]

Syms are great, but sometimes just solving stuff works as well. This is just a one-parameter problem in the known quantities. With a bit of algebra you can arrive at theta = acos(A) where parameter A is shown below.

E = 447000;   deltay = .00127;   q = 40;   t = .000211;   R = .0525;
v = .4
D = deltay*(2*E*t)/(q*(1-v));
A = R^2/D -1;
theta = acos(A);
R2 = R/sin(theta)       % definition of R2
% check other equation
(deltay*(2*E*t)/(q*(1-v))) - R2^2*(1-cos(theta))
ans = 5.2042e-18

acos gives

0 <= theta <= pi/2  for  0 <= A <= 1
pi/2 <= theta <= pi for -1 <= A <= 0
theta complex for 1 < abs( A)

so you just need appropriate values of v to make this work, assuming v is the constant that you are varying.

[monarch]

I like the idea, but that gives me a different value for R2 than the 2 I got without trying to do human algebra (which I haven't done in 20+ years!) I get R2 = -2.2342 or 2.2342 the original way...and R2 = 0.0761 this new way. Both of them are potentially plausible. The main thing I need to be sure is that I'm only getting solutions for which theta is > 0 and < pi/2 (this is basically one quadrant of a sphere I'm modeling here).

[David Goodmanson]

I assume that the other calculation is also for v = .4, is that correct? Then it looks like theta would have been .0235 radians, is that also the case?

[monarch]

Here's what I'm currently working with...I changed some of the constants.

v=.47 
E = 447000;  
deltay = 0.0170
q = 9806.65; 
t = .00682; 
R = .0525;

syms R2 
syms theta
eqns = [deltay == (q*(R2^2)*(1-v)*(1-cos(theta)))/(2*E*t), sin(theta) == R/R2];
vars = [R2 theta];
[solR2,soltheta] = solve(eqns, vars)
R2 = double(solR2);
theta = double(soltheta);

This gives me R2 =

   -0.1035
    0.1035

theta =

    9.9568
    2.6096

[monarch]

v=.47 
E = 447000;  
deltay = 0.0170
q = 9806.65; 
t = .00682; 
R = .0525;

D = deltay*(2*E*t)/(q*(1-v));
A = R^2/D -1;
theta = acos(A);
R2 = R/sin(theta)       % definition of R2
% check other equation
check=(deltay*(2*E*t)/(q*(1-v))) - R2^2*(1-cos(theta))

gives me theta =

    2.6096

R2 =

    0.1035

check =

   1.3878e-17

[Walter Roberson]

In R2017b using that code gives me

         -2.60958116291844
          2.60958116291844

for theta. The two solutions are exact negatives of each other (same is true for solR2, the solutions are +/- a number.)

I have checked with a different software package and confirmed that the answers are necessarily negatives of each other and that there is a closed form solution.

[monarch]

Thanks for the check. I'm still disappointed that I have yet to get an answer for theta that is less than pi/2 even with all my futzing with constants. I'm satisfied though that it's not due to me coding improperly...perhaps I'm violating the assumptions of the model to too great a degree!

[Walter Roberson]

Small v, small q. For example

syms v E deltay q t R R2 theta
eqns = [deltay == (q*(R2^2)*(1-v)*(1-cos(theta)))/(2*E*t), sin(theta) == R/R2];
vars = [R2 theta];
[solR2,soltheta] = solve(eqns, vars);
subs(soltheta, {'E','R','deltay','t', 'q', 'v'}, {447000, 0.0525, 0.01, 0.01, 40, 1/2})

gets you a theta about 3.106

[David Goodmanson]

Hi monarch

It looks like you are fine with changing the value of q, so try

% first four constants same as in original question
E = 447000;  deltay = .00127;  t = .000211;  R = .0525;
v = 1/2;
q = 280;

D = deltay*(2*E*t)/(q*(1-v));
A = R^2/D -1
theta = acos(A)
R2 = R/sin(theta)

% check that the original equations are satisfied

(deltay*(2*E*t)/(q*(1-v))) - R2^2*(1-cos(theta))
     ans = 2.1684e-19
sin(theta) - R/R2
     ans =  -1.1102e-16

The result of this was

theta = 0.9138      R2 = 0.0663

and theta is definitely less than pi/2. I used a different means of solution but the equations are verified for equality.

As I mentioned in a previous comment, you can use any combination of constants that you like, and as long as the constant A above satisfies 0 < A < 1 you will get an angle between 0 and pi/2. If you don't satisfy that condition for A, you wont get theta less than pi/2, simple as that. In this case

A = 0.6107

[monarch]

If I vary q, I also have to vary delta...and then theta is no longer less than pi/2. Here it is with a combo closer to the 280 you suggested. the v, t, and q in my original question are flat wrong (I realized after I asked)...but delta y and q vary together. This odd value of theta is quite possibly a result of me violating assumptions of the model....except that the ratios of R to R2 are coming out as I expect.

v=.47; 
E = 447000;  
deltay = .0048
q = 239.980263158 
t =.00715 
R = .0525; 

syms R2 
syms theta
eqns = [deltay == (q*(R2^2)*(1-v)*(1-cos(theta)))/(2*E*t), sin(theta) == R/R2];
vars = [R2 theta];
[solR2,soltheta] = solve(eqns, vars)
R2 = double(solR2);
theta = double(soltheta);

final_theta = theta(find(theta > 0 & theta <= 4))
final_R2 = R/(sin(final_theta))

final_theta =

    2.9903

final_R2 =

    0.3483

[David Goodmanson]

It would help if you could explain what this code is actually modeling. That would make some of the relationships of the variables clearer. For example how is q related to deltay? Can you vary E or t or R if you want, independent of the others? Are there restrictions on the sizes of some of these variables? If there is any independence at all, then it is not hard to get back to theta less than pi/2 again.

[Walter Roberson]

Sorry, my earlier "Small v, small q" response was based upon misreading the Pi/2 part.

[Walter Roberson]

If you use

{'E','R','deltay','t'}, {447000, 0.0525, 0.01, 0.01}

with

eqns = [deltay == (q*(R2^2)*(1-v)*(1-cos(theta)))/(2*E*t), sin(theta) == R/R2];

then as v approaches 1 from below, q can become arbitrarily large and as that happens the solutions can be less than Pi/2

For v = 0.4, you can get arbitrarily close to 0 as q approaches 47680000/441 from below.

[monarch]

The equations are for a thin-walled spherical shell with uniform internal or external pressure, q force=unit area; tangential edge support (Roark's equations for stress and strain Table 13.1 3a). All values should be in SI units. One underlying assumption is that R2 is more than 10x t and I haven't violated that yet. There aren't any equations for the relationship between q and deltay (other than the above) but I have different experimental conditions for which those are obtained and they vary together. V is nu (poisson's ratio for the material) which is the same under all conditions - I had my equations wrong to obtain it originally). E is elastic modulus and we are making the assumption that it is constant for all conditions, however, I have some wiggle room to vary it within a range of +/- 97kPa around the 447kPa (the mean). t has very little wiggle room: basically one set of conditions where it is 6.82–7.48mm and another where it is 5.76–6.93mm. R has a range of 5 to 7cm that is independent of all other variables.

[Walter Roberson]

monarch, which value ranges would you like me to explore?

[monarch]

There's still something wrong other than what constants I'm inputing.

My thetas are not matching up when I check them with the definition of R2.

v=.47; 
E = 447000;  
deltay = .0048 
q = 239.980263158 
t =.00715 
R = .0525; 

D = deltay*(2*E*t)/(q*(1-v));
A = R^2/D -1;
theta = acos(A)
R2 = R/sin(theta)       % definition of R2
% check other equation
check=(deltay*(2*E*t)/(q*(1-v))) - R2^2*(1-cos(theta))

final_asum_grtr_than_10 = R2/t;
check_R2 = R/(sin(theta));
check_theta = asin(R/R2);
check_deltay =(q*(R2.^2)*(1-v)*(1-cos(theta)))/(2*E*t);

[monarch]

v=.47;
E = 447000;  
deltay = .0048 
q = 239.980263158 
t =.00715 
R = .0525; 

syms R2 
syms theta
eqns = [deltay == (q*(R2.^2)*(1-v)*(1-cos(theta)))/(2*E*t), sin(theta) == R/R2];
vars = [R2 theta];
[solR2,soltheta] = solve(eqns, vars)
R2 = double(solR2);
theta = double(soltheta);

final_theta = theta(find(theta > 0 & theta <= 4))
final_R2 = R/(sin(final_theta))
final_theta_d = radtodeg(final_theta)
final_asum_grtr_than_10 = final_R2/t;
final_theta_draw = pi - final_theta;
check_R2 = R/(sin(theta(2)));
check_theta = asin(R/R2);
check_deltay =(q*(R2.^2)*(1-v)*(1-cos(final_theta_draw)))/(2*E*t);

[monarch]

my check_theta and theta should at least have one value in common (using solve function (below), or the algebraic method (above))....in the 2nd block of code, final_theta_draw is pi minus theta (which looks more like the theta I was expecting) except when I plug that back in to the equation for deltay, I don't get the same deltay.

[David Goodmanson]

If R/R2 is positive (and less than one) then asin outputs into the first quadrant, so you get check_theta = 0.1513. But since sin(theta) = sin(pi-theta), (pi-theta) is also a valid solution for asin. Indeed if you compare check_theta with theta from the acos calculation, then if you check whether (pi - check_theta) = theta, that checks out.

Geometrically, is check_theta really the angle that counts? It's certainly in the range you wanted.

You did not yet mention what deltay is. Is it a deflection?

[monarch]

Yes - deflection. "deltay is change in the height dimension y; y is length of cylindrical or conical shell and is also used as a vertical position coordinate, positive upward."

I agree that check_theta is the one that makes most sense to me geometrically. However, plugging it back in to the deltay equation gives me a different deltay than the fixed one I input from the beginning.

The is the plot I get for theta on x axis, deltay on y axis for the above set of conditions (where deltay is supposed to be .0048).

Answer 2

[monarch]

For the curious...the mystery has been solved. For this material, E varies under these conditions (but we can't measure it in this experiment) so I solved for E under various conditions and found that it's nowhere near the constant for E I was using (that was derived from a different experiment). When I use the derived E and feed everything back in, the thetas match. HOORAY!

[David Goodmanson]

well, I guess ... I would be pretty careful about using a derived E unless you have support from experiment, which it sounds like you don't. In fact there is probably every reason to keep to the value of E that does have experimental support and look at other possibilites. For one thing, I went and looked at Roark's formula. I will add to this comment in a few hours but for now want to ask since q is positive then the shell was pressurized on the inside, is that true?

[monarch]

Yes that is true the shell was pressurized on the inside...and I understand your concern with using the derived E. It does indeed problems. I was just happy to see the theta matching...but yes, it is problematic and most likely pointing to some serious assumption violating going on!

[David Goodmanson]

Could you describe the experimental setup, for example is it a hemisphere, and how you are measuring y, and how is your angle theta measured or defined?

I'm not sure what you don't like about the plot that you posted since it looks like it peaks out at .0048 just like it should.