Why the below mentioned function plotting is unsuccessful?

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]);

Answer 1

[Walter Roberson]

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

[Sarbhanu Saha]

Thank you for your answer. I want to integrate the whole function for each value of ri ( but a fixed upperbound ro value and variable lowerbound ri value) so that i can get the corresponding value of Mkrit and so will be able to plot a graph between Mkrit and ri.

It is my 1st time seeing arrayfun(), so i will look into it.

[Walter Roberson]

Here is an efficiency trick:

Instead of integrating the vector q2inint for ri(K) to ro, you can sort the ri and integrate in pieces and sum to get the individual parts.

Assuming all ri are <= ro

[sorted_bounds, sortidx] = sort(ri);
nint = length(sorted_bounds);
sorted_bounds(end+1) = ro;
partials = zeros(1, nint, 'sym');
for K = 1 : nint
  partials(K) = vpaintegral(q2inint, r, sorted_bounds(K:K+1));
end
totals = cumsum(partials);
q2 = -1 .* totals(sortidx);

... of course if you know that your ri is in sorted order already, you can skip some of these steps.

[Sarbhanu Saha]

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]

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)