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! :)
采纳的回答
Walter Roberson
2017-12-1
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.
24 个评论
monarch
2017-12-1
Walter Roberson
2017-12-1
double(solR2), double(soltheta)
monarch
2017-12-1
Walter Roberson
2017-12-1
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
2017-12-1
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
2017-12-1
David Goodmanson
2017-12-2
编辑：David Goodmanson
2017-12-2
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
2017-12-2
monarch
2017-12-2
Walter Roberson
2017-12-2
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
2017-12-2
Walter Roberson
2017-12-2
编辑：Walter Roberson
2017-12-3
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
2017-12-3
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
2017-12-3
David Goodmanson
2017-12-3
编辑：David Goodmanson
2017-12-3
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
2017-12-3
Sorry, my earlier "Small v, small q" response was based upon misreading the Pi/2 part.
Walter Roberson
2017-12-3
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.
Walter Roberson
2017-12-3
monarch, which value ranges would you like me to explore?
monarch
2017-12-3
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
2017-12-3
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
2017-12-3
David Goodmanson
2017-12-3
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
2017-12-3
编辑：Walter Roberson
2017-12-3
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).
更多回答（1 个）
monarch
2017-12-3
0 个投票
3 个评论
David Goodmanson
2017-12-4
编辑：David Goodmanson
2017-12-4
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
2017-12-4
David Goodmanson
2017-12-4
编辑：David Goodmanson
2017-12-7
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Number Theory 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
