I want to find omega value in the equation.
can i anybody tell me how to solve this equation.

 채택된 답변

Walter Roberson
Walter Roberson 2021년 6월 21일

0 개 추천

The search turns out to take over 20 seconds for each eload value, and even then it seems to fail.
Notice I ended at 0.002, which is the most I could do with this online version without timing out.
tic
syms omega
a=3.666;
k=0.467;
s=100;
r=1.4;
ef=0.08483;
N=(a*k)/(1+(a*k));
eloadvals = 0:0.001:0.002;
for eloadidx = 1 : length(eloadvals)
eload = eloadvals(eloadidx);
c1=(1-N)*(((r*a)/2)+(r/2)-1)-(r-1);
c2=(N*(r/(2*a))+(r/2)-1);
Ac2=((r*a)*(1-N))/2;
Ah1=(r*N)/(2*a);
A=-8*(1-N)*c2*sin(pi-atan(s/((1-N)*omega)))+8*N*c1*sin(2*(pi-atan(1/(N*omega))))+(32*N*ef+32*N*eload);
B=16*(1-N)*(-1*(((1-N)*omega)/(2*s))*(((Ah1*c1*sin(pi-atan(1/(N*omega)))*cos(pi-atan(1/(N*omega))))+((((N*s)/((1-N)*omega))+(N*(r/2-1)))*c2*cos((atan((((1-N)*((r/2)-1))/(s))*omega))+(pi-atan(s/((1-N)*omega))))*cos(pi-atan(s/((1-N)*omega)))))))*c2*sin(2*(pi-atan(s/((1-N)*omega))))+(16*(((N*omega)/2)*((((((1-N)/(N*omega))+((1-N)*(((r/2)-1)-(r-1))))*c1*cos((atan((((1-N)*((r/2)-1)-(r-1))/(1-N))*(N*omega)))+(pi-atan(1/(N*omega))))*cos(atan((((1-N)*((r/2)-1)-(r-1))/(1-N))*(N*omega))))+(Ac2*c2*sin(pi-atan(s/((1-N)*omega)))*cos(pi-atan(s/((1-N)*omega)))))))+2*c1)*N*c1*sin(2*(pi-atan(1/(N*omega))))+N*c1*c1*sin(4*(pi-atan(1/(N*omega))));
thissol = vpasolve(A==B);
if length(thissol) >= 1
sol(eloadidx) = thissol(1);
else
sol(eloadidx) = nan;
end
end
toc
Elapsed time is 54.511198 seconds.
plot(eloadvals, sol)
Nothing visibile in the plot because vpasolve() failed.

댓글 수: 7

Walter Roberson
Walter Roberson 2021년 6월 21일
My tests show that if the question is to find the omega for which A == B, then in order for there to be solutions, eload must be negative, below approximately -.014
You can derive this by taking symbolic eload, and calculate A-B and solving for eload, and start plotting. You will see a discontinuity, and that it never crosses 0; the closest to zero is at approximately omega = 1.42.
Ali
Ali 2021년 6월 21일
right sir can we find eload (positve or negative) at which omega is only one root?
Walter Roberson
Walter Roberson 2021년 6월 21일
If I recall correctly, the plot of eload in terms of omega has a discontinuity at 0, and rises from there, then falls to a minima at about omega = 1.42, and then increases again. The places where omega has only one root would then be the point of the minima (at about 1.42), together with the points that are above the local maxima that it rises to after the discontinuity at 0... as a quick mental memory, it was somewhere around omega = 3
I would need to recheck what happened on the other side of the discontinuity.
We could probably identify that singular point, and the range where it goes above the local maxima. But I have to ask why you would want to know those ??
I cannot quite get three plots without timing out, so the plots will be split between two posts.
syms omega
Q = @(v) sym(v);
a = Q(3.666);
k = Q(0.467);
s = Q(100);
r = Q(1.4);
ef = Q(0.08483);
N = (a*k)/(1+(a*k));
syms eload
Pi = Q(pi);
c1 = (1-N)*(((r*a)/2)+(r/2)-1)-(r-1);
c2 = (N*(r/(2*a))+(r/2)-1);
Ac2 = ((r*a)*(1-N))/2;
Ah1 = (r*N)/(2*a);
A = -8*(1-N)*c2*sin(Pi-atan(s/((1-N)*omega)))+8*N*c1*sin(2*(Pi-atan(1/(N*omega))))+(32*N*ef+32*N*eload);
B = 16*(1-N)*(-1*(((1-N)*omega)/(2*s))*(((Ah1*c1*sin(Pi-atan(1/(N*omega)))*cos(Pi-atan(1/(N*omega))))+((((N*s)/((1-N)*omega))+(N*(r/2-1)))*c2*cos((atan((((1-N)*((r/2)-1))/(s))*omega))+(Pi-atan(s/((1-N)*omega))))*cos(Pi-atan(s/((1-N)*omega)))))))*c2*sin(2*(Pi-atan(s/((1-N)*omega))))+(16*(((N*omega)/2)*((((((1-N)/(N*omega))+((1-N)*(((r/2)-1)-(r-1))))*c1*cos((atan((((1-N)*((r/2)-1)-(r-1))/(1-N))*(N*omega)))+(Pi-atan(1/(N*omega))))*cos(atan((((1-N)*((r/2)-1)-(r-1))/(1-N))*(N*omega))))+(Ac2*c2*sin(Pi-atan(s/((1-N)*omega)))*cos(Pi-atan(s/((1-N)*omega)))))))+2*c1)*N*c1*sin(2*(Pi-atan(1/(N*omega))))+N*c1*c1*sin(4*(Pi-atan(1/(N*omega))));
sol = solve(A-B, eload)
sol = 
double(limit(sol, omega, -inf))
ans = -0.0848
double(limit(sol, omega, inf))
ans = -0.0848
fplot(sol, [-10 10]); xlabel('omega'); ylabel('eload')
%fplot(sol, [1.35 1.5]); xlabel('omega'); ylabel('eload')
%fplot(sol, [-1.5 -1.35]); xlabel('omega'); ylabel('eload')
syms omega
Q = @(v) sym(v);
a = Q(3.666);
k = Q(0.467);
s = Q(100);
r = Q(1.4);
ef = Q(0.08483);
N = (a*k)/(1+(a*k));
syms eload
Pi = Q(pi);
c1 = (1-N)*(((r*a)/2)+(r/2)-1)-(r-1);
c2 = (N*(r/(2*a))+(r/2)-1);
Ac2 = ((r*a)*(1-N))/2;
Ah1 = (r*N)/(2*a);
A = -8*(1-N)*c2*sin(Pi-atan(s/((1-N)*omega)))+8*N*c1*sin(2*(Pi-atan(1/(N*omega))))+(32*N*ef+32*N*eload);
B = 16*(1-N)*(-1*(((1-N)*omega)/(2*s))*(((Ah1*c1*sin(Pi-atan(1/(N*omega)))*cos(Pi-atan(1/(N*omega))))+((((N*s)/((1-N)*omega))+(N*(r/2-1)))*c2*cos((atan((((1-N)*((r/2)-1))/(s))*omega))+(Pi-atan(s/((1-N)*omega))))*cos(Pi-atan(s/((1-N)*omega)))))))*c2*sin(2*(Pi-atan(s/((1-N)*omega))))+(16*(((N*omega)/2)*((((((1-N)/(N*omega))+((1-N)*(((r/2)-1)-(r-1))))*c1*cos((atan((((1-N)*((r/2)-1)-(r-1))/(1-N))*(N*omega)))+(Pi-atan(1/(N*omega))))*cos(atan((((1-N)*((r/2)-1)-(r-1))/(1-N))*(N*omega))))+(Ac2*c2*sin(Pi-atan(s/((1-N)*omega)))*cos(Pi-atan(s/((1-N)*omega)))))))+2*c1)*N*c1*sin(2*(Pi-atan(1/(N*omega))))+N*c1*c1*sin(4*(Pi-atan(1/(N*omega))));
sol = solve(A-B, eload)
sol = 
double(limit(sol, omega, -inf))
ans = -0.0848
double(limit(sol, omega, inf))
ans = -0.0848
%fplot(sol, [-10 10]); xlabel('omega'); ylabel('eload')
fplot(sol, [1.35 1.5]); xlabel('omega'); ylabel('eload')
fplot(sol, [-1.5 -1.35]); xlabel('omega'); ylabel('eload')
Walter Roberson
Walter Roberson 2021년 6월 21일
So look at the top plot. The left and right halfs are symmetric around eload = -0.0848
You can see near the discontinuity at 0 that the left side close to 0 shares values with the right side close to 0, so that section is not unique.
If we look down the discontinuity on the left side close to 0, to the point where the one on the right close to 0 gives out, and then mentally follow that level near -0.1 to the left, we can see that we are not unique, there is the other side of the valley that rises up again. How high does the left valley rise? Well it rises to double(limit(sol, omega, -inf)) = -0.0848, and by following across towards the discontinuity at 0, we can see that that entire area gets us to a matching edge just to the left of 0.
What point does not have a duplicate? Only the lowest part of the valley, at approximately omega = -1.42 corresponding to eload of about -0.15342 : if you follow that one point left and right, there is no corresponding point on either side of the discontinuity.
Likewise, when we look on the right side of the discontinuity, we have a balance point at about 1.42 (same distance away from 0 as the other one), corresponding to eload of about -0.01403 .
Those are the only two omega that have unique solutions: +/- 1.42 (approximately) corresponding to eload of about -0.01403 and eload of about -0.15342 .
If there are any stable points, it would be only those two.
If omega were being controlled by some force pushing towards a balance, then if you were to disturb omega a little, it should rebalance itself. But if not... well, if there is no restoring pressure on omega, then if the system gets knocked slightly so that omega varies, then you would start getting two eload solutions.
Ali
Ali 2021년 6월 22일
Thanks a lot sir

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Reshma Nerella
Reshma Nerella 2021년 6월 21일

0 개 추천

Hi,
You can use solve function from symbolic math toolbox to solve equations and obtain the values of variables.
Refer to the documentation page for more informataion and examples: Solve.

카테고리

태그

질문:

Ali
2021년 6월 18일

편집:

2024년 7월 11일

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by