Why the below mentioned function plotting is unsuccessful?

17 次查看（过去 30 天）
Sarbhanu Saha
Sarbhanu Saha 2021-8-17
评论： Walter Roberson 2021-8-28
I want to create a plot for the below given function. But, I always get the same error. The code for the plotting is:
syms a m n b r s phi E D ri ro u;
ro = 80;
E = 210000;
s = 1;
u = 3./10;
D = (E.*s.^3)./(12.*(1-u.^2));
a = 0.005;
m = 1.1492;
n = 2;
b = 0.6158;
Pi = sym(pi);
assume(r >= ri);
assume(r <= ro);
ri = linspace(5, 15, 100);
z = exp(a.*r).*((r-ri).^2).*((r-ro).^2).*cos(n.*(tan(b).*log(ri./r) + phi))./(r.^m);
dzr = diff(z,r);
dzrr = diff(dzr,r);
dzp = diff(z,phi);
dzpp = diff(dzp,phi);
dzrp = diff(dzr,phi);
q2inside = ((1./r.^3).*dzr.*dzp).*r;
q2inint = vpaintegral(q2inside,phi,[0 2.*Pi]);
q2 = (-1).*vpaintegral(q2inint,r,[ri ro]);
q1inside = (D.*((dzrr + (1./r).*dzr + (1./r.^2).*dzpp).^2 + 2.*(1-u).*(((1./r).*dzrp - (1./r.^2).*dzp).^2 - dzrr.*((1./r).*dzr + (1./r.^2).*dzpp)))).*r;
q1inint = vpaintegral(q1inside,r,[ri ro]);
q1 = Pi.*vpaintegral(q1inint,phi,[0, 2.*Pi]);
Mkrit = q1./q2;
q = double(abs(Mkrit));
plot(ri, q)
The error that i am encountering is :
Operands to the logical and (&&) and or (||) operators must be convertible to logical scalar values.
Error in sym/vpaintegral (line 182)
(isinf(b) && isempty(symvar(b)) && ~isreal(b))
Error in Graph_Try (line 23)
q2 = (-1).*vpaintegral(q2inint,r,[ri ro]);

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2021-8-17
Your ri is a vector. You cannot use a vector as one of the bounds for vpaintegral, with the exception that you can supply a vector of two elements as the bounds instead of separate upper and lower boundaries.
Your expression to be integrated is a vector. Are you wanting to integrate each element of the vector with respect to each different lower bound, or are you wanting to integrate each one with respect to only its corresponding lower bound? Remember that you can arrayfun() calls to vpaintegral
  4 个评论
显示 2更早的评论
Sarbhanu Saha
Sarbhanu Saha 2021-8-28
Thank you for your response. I learned about integration with varying limits. Since ri is a vector in my previous code, i so decided to make the graph using for loop. But now it does not show any point in the graph. Here is the following code:
syms a m n b r s phi E D ri ro u x;
ro = 80;
E = 70000;
s = 1;
u = 3./10;
D = (E.*s.^3)./(12.*(1-u.^2));
a = 0.005;
m = 1.1492;
n = 2;
b = 0.6158;
Pi = sym(pi);
assume(r >= ri);
assume(r <= ro);
for x = 5:0.5:15
ri = x;
z = exp(a.*r).*((r-ri).^2).*((r-ro).^2).*cos(n.*(tan(b).*log(ri./r) + phi))./(r.^m);
dzr = diff(z,r);
dzrr = diff(dzr,r);
dzp = diff(z,phi);
dzpp = diff(dzp,phi);
dzrp = diff(dzr,phi);
q2inside = ((1./r.^3).*dzr.*dzp).*r;
q2inint = vpaintegral(q2inside,phi,[0 2.*Pi]);
q1inside = (D.*((dzrr + (1./r).*dzr + (1./r.^2).*dzpp).^2 + 2.*(1-u).*(((1./r).*dzrp - (1./r.^2).*dzp).^2 - dzrr.*((1./r).*dzr + (1./r.^2).*dzpp)))).*r;
q1inint = vpaintegral(q1inside,phi,[0 2.*Pi]);
q2 = (-1).*vpaintegral(q2inint,r,[ri ro]);
q1 = Pi.*vpaintegral(q1inint,r,[ri ro]);
Mkrit = q1./q2;
q = double(abs(Mkrit));
end
plot(ri, q)
Can you suggest me what to do so that the graph is visible?
Walter Roberson
Walter Roberson 2021-8-28
Ran in:
syms a m n b r s phi E D ri ro u x;
ro = 80;
E = 70000;
s = 1;
u = 3./10;
D = (E.*s.^3)./(12.*(1-u.^2));
a = 0.005;
m = 1.1492;
n = 2;
b = 0.6158;
Pi = sym(pi);
assume(r >= ri);
assume(r <= ro);
xvals = 5:0.5:15;
numx = length(xvals)
numx = 21
for xidx = 1 : numx
x = xvals(xidx);
ri = x;
z = exp(a.*r).*((r-ri).^2).*((r-ro).^2).*cos(n.*(tan(b).*log(ri./r) + phi))./(r.^m);
dzr = diff(z,r);
dzrr = diff(dzr,r);
dzp = diff(z,phi);
dzpp = diff(dzp,phi);
dzrp = diff(dzr,phi);
q2inside = ((1./r.^3).*dzr.*dzp).*r;
q2inint = vpaintegral(q2inside,phi,[0 2.*Pi]);
q1inside = (D.*((dzrr + (1./r).*dzr + (1./r.^2).*dzpp).^2 + 2.*(1-u).*(((1./r).*dzrp - (1./r.^2).*dzp).^2 - dzrr.*((1./r).*dzr + (1./r.^2).*dzpp)))).*r;
q1inint = vpaintegral(q1inside,phi,[0 2.*Pi]);
q2 = (-1).*vpaintegral(q2inint,r,[ri ro]);
q1 = Pi.*vpaintegral(q1inint,r,[ri ro]);
Mkrit = q1./q2;
q(xidx) = double(abs(Mkrit));
end
plot(xvals, q)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by