how to plot a function where there is inside 1x1 sym?

I want to plot T(OMEGA):

T=real(1i*n*alfak*(besselj(np+1,alfak))^2*(conj(A3)*B3-A3*conj(B3)))

where:

• np is known
• n is known
• alfak is known
• besselj(np+1, alfak) is the bessel function of the fisrt kind of order np+1

But A3 and B3 depends on OMEGA and they appears as 1x1 sym in the workspace.

I used this code:

fplot(T)
grid on

But this message appears:

Error in fplot>vectorizeFplot (line 191)

hObj = cellfun(@(f) singleFplot(cax,{f},limits,extraOpts,args),fn{1},'UniformOutput',false);

Error in fplot (line 161)

hObj = vectorizeFplot(cax,fn,limits,extraOpts,args);

Error in calcolo_grafico_coppia (line 100)

fplot(T)

I post the entire code if you want to try (the code is for sure correct, only the function plot is of interest for the question):

clear all
clc 


a=10e-3;
b=10e-3;
c=5e-3;    
d=5e-3;
e=8e-3;

z1=a;
z2=z1+b;
z3=z2+c; 
z4=z3+d;
z5=z4+e;

Br=1.25;  %T

R1=25e-3;
R2=65e-3;
R3=90e-3;

u0=4*pi*1e-7;
urb=1000;

sigmab=7e6;  %M=1000000
sigmac=57e6;

syms OMEGA

alfa=0.9; 

p=4;

%TORQUE T
A=(1:2:10); %vector in order to take n odd
summeT=0.0;
for n=A(1):A(length(A))  
    for k=1:10
    np=n*p;
    kind=1;
    alfaks = (besselzero(np, k, kind))/R3;  
    alfak=alfaks(k);
    
    gamma4=sqrt((alfak^2)+1i*np*sigmac*u0*OMEGA);
    gamma5=sqrt((alfak^2)+1i*np*sigmab*urb*u0*OMEGA); %devo moltiplicare u0*urb?
    
    syms r; %devo definire la variabile r
    Mnksym=((8*Br)/(n*pi*u0*R3^2*(besselj(np+1,alfak*R3))^2))*sin(n*alfa*pi/2)*int(r*besselj(np,alfak*r),r,R1,R2);
    Mnk=double(Mnksym);  %in order to not  have 1x1 sym 
   
    syms A1 B1 A2 B2 A3 B3 A4 B4 A5 B5 A6 B6 A7 B7 A8 B8 A9 B9 A10 B10;
    eqn1 = A1-B1==0; 
    eqn2 = A1*exp(alfak*z1)+B1*exp(-alfak*z1)==A2*exp(alfak*z1)+B2*exp(-alfak*z1);
    eqn3 = A1*exp(alfak*z1)-B1*exp(-alfak*z1)==(1/urb)*(A2*exp(alfak*z1)+B2*exp(-alfak*z1)+(Mnk/alfak));
    eqn4 = A2*exp(alfak*z2)+B2*exp(-alfak*z2)==A3*exp(alfak*z2)+B3*exp(-alfak*z2);
    eqn5 = A2*exp(alfak*z2)-B2*exp(-alfak*z2)==A3*exp(alfak*z2)-B3*exp(-alfak*z2)+(Mnk/alfak);
    eqn6 = A3*exp(alfak*z3)+B3*exp(-alfak*z3)==(gamma4/(n^2*p^2*sigmac*u0*OMEGA))*(A4*exp(gamma4*z3)-B4*exp(-gamma4*z3));
    eqn7 = A3*exp(alfak*z3)-B3*exp(-alfak*z3)==(alfak/(n^2*p^2*sigmac*u0*OMEGA))*(A4*exp(gamma4*z3)+B4*exp(-gamma4*z3));
    eqn8 = A4*exp(gamma4*z4)+B4*exp(-gamma4*z4)==(sigmac/sigmab)*(A5*exp(gamma5*z4)+B5*exp(-gamma5*z4));
    eqn9 = gamma4*A4*exp(gamma4*z4)-gamma4*B4*exp(-gamma4*z4)==(sigmac/(sigmab*urb))*(gamma5*A5*exp(gamma5*z4)-gamma5*B5*exp(-gamma5*z4));
    eqn10 = A5*exp(gamma5*z5)+B5*exp(-gamma5*z5)==0;
    
    sol = solve([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10], [A1 B1 A2 B2 A3 B3 A4 B4 A5 B5]);
    
    summeT=summeT+(1i*n*alfak*(besselj(np+1,alfak*R3))^2*(conj(A3)*B3-A3*conj(B3)));
    end

end


T=(pi/2)*u0*R3^2*p*real(summeT);


fplot(T)
grid on

[Walter Roberson]

Are you using https://www.mathworks.com/matlabcentral/fileexchange/48403-bessel-zero-solver for besselzero() ?

[Luigi Stragapede]

yes, the only problem is the plot of T.

All the previous code has been controlled and it is ok!

Answer 1

[Walter Roberson]

    summeT=summeT+(1i*n*alfak*(besselj(np+1,alfak*R3))^2*(conj(A3)*B3-A3*conj(B3)));

That line is not using A3 and B3 from the solution, which are sol.A3 and sol.B3

[Walter Roberson]

Note: the code is now Very Slow. Very slow.

[Walter Roberson]

A=(1:2:10); %vector in order to take n odd
summeT=0.0;
for n=A(1):A(length(A))  

What is the point of creating A that way, to contain [1 3 5 7 9] ? The only place you use the content of A is in the |for n=A(1):A(length(A)) which is 1:9 without skipping any values -- in which case why not just make A = [1 9] ?

I suspect you should be using

for n=A

to have n iterate over the values stored in A.

[Luigi Stragapede]

I want to consider only n odd integer from 1 to 9. In fact A(1)=1, A(2)=3, A(3)=5, A(4)=7 and A(5)=A(lenght(A))=9.

Do you think it is wrong?

I am running the code, but it is very very long. Is it normal?

[Walter Roberson]

With some adjustments to try to simplify the generated equations:

Q = @(v) sym(v);
V16 = @(v) vpa(v, 16);

a=Q(10e-3);
b=Q(10e-3);
c=Q(5e-3);
d=Q(5e-3);
e=Q(8e-3);

z1=a;
z2=z1+b;
z3=z2+c;
z4=z3+d;
z5=z4+e;

Br=Q(1.25);  %T

R1=Q(25e-3);
R2=Q(65e-3);
R3=Q(90e-3);

u0=Q(4)*Q(pi)*Q(1e-7);
urb=Q(1000);

sigmab=Q(7e6);  %M=1000000
sigmac=Q(57e6);

syms OMEGA real

alfa=Q(0.9);

p=(4);

%TORQUE T
A=(1:2:10); %vector in order to take n odd
summeT=Q(0.0);
for n=A
    for k=1:10
        np=n*p;
        kind=1;
        alfaks = (besselzero(np, k, kind))/R3;
        alfak=alfaks(k);
        
        gamma4=sqrt((alfak^2)+1i*np*sigmac*u0*OMEGA);
        gamma5=sqrt((alfak^2)+1i*np*sigmab*urb*u0*OMEGA); %devo moltiplicare u0*urb?
        
        syms r; %devo definire la variabile r
        Mnksym=((8*Br)/(n*Q(pi)*u0*R3^2*(besselj(np+1,alfak*R3))^2))*sin(n*alfa*Q(pi)/2)*int(r*besselj(np,alfak*r),r,R1,R2);
        Mnk=double(Mnksym);  %in order to not  have 1x1 sym
        
        syms A1 B1 A2 B2 A3 B3 A4 B4 A5 B5 A6 B6 A7 B7 A8 B8 A9 B9 A10 B10;
        eqn1 = A1-B1==0;
        eqn2 = A1*exp(alfak*z1)+B1*exp(-alfak*z1)==A2*exp(alfak*z1)+B2*exp(-alfak*z1);
        eqn3 = A1*exp(alfak*z1)-B1*exp(-alfak*z1)==(1/urb)*(A2*exp(alfak*z1)+B2*exp(-alfak*z1)+(Mnk/alfak));
        eqn4 = A2*exp(alfak*z2)+B2*exp(-alfak*z2)==A3*exp(alfak*z2)+B3*exp(-alfak*z2);
        eqn5 = A2*exp(alfak*z2)-B2*exp(-alfak*z2)==A3*exp(alfak*z2)-B3*exp(-alfak*z2)+(Mnk/alfak);
        eqn6 = A3*exp(alfak*z3)+B3*exp(-alfak*z3)==(gamma4/(n^2*p^2*sigmac*u0*OMEGA))*(A4*exp(gamma4*z3)-B4*exp(-gamma4*z3));
        eqn7 = A3*exp(alfak*z3)-B3*exp(-alfak*z3)==(alfak/(n^2*p^2*sigmac*u0*OMEGA))*(A4*exp(gamma4*z3)+B4*exp(-gamma4*z3));
        eqn8 = A4*exp(gamma4*z4)+B4*exp(-gamma4*z4)==(sigmac/sigmab)*(A5*exp(gamma5*z4)+B5*exp(-gamma5*z4));
        eqn9 = gamma4*A4*exp(gamma4*z4)-gamma4*B4*exp(-gamma4*z4)==(sigmac/(sigmab*urb))*(gamma5*A5*exp(gamma5*z4)-gamma5*B5*exp(-gamma5*z4));
        eqn10 = A5*exp(gamma5*z5)+B5*exp(-gamma5*z5)==0;
        
        sol = solve([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10], [A1 B1 A2 B2 A3 B3 A4 B4 A5 B5]);
        
        summeT = summeT + V16((1i*n*alfak*(besselj(np+1,alfak*R3))^2*(conj(sol.A3)*sol.B3-sol.A3*conj(sol.B3))));
    end
    summeT = V16(summeT);
end


T = V16( (Q(pi)/2)*u0*R3^2*p*real(summeT) );

Getting to this point does not take too long.

Notice that the plotting is not present. fplot() can be very slow. It can be a lot more efficient to

OM = V16(0 : .5 : 5);
Y = subs(T, OMEGA, OM);
plot(OM, Y)

When I tested, I got a straight line, but the T equations is clearly non-linear. That could indicate that it just is not obvious that there is a curve in the limited range that I plotted. Even just taking a single derivative takes a long time, so estimating the curve has to be done numerically.

[Luigi Stragapede]

I have to obtain this RED curve:

In particular the range must be from 0 to 209 (in the graph the range is from 0 to 2000 because it is rpm, it only a different scale).

If I reduce the iterations, I obtain more quickly the graph, but I nedd at least 7 iteretions for k and n in order to have an accurate result.

Now I try your code. You are very kind

[Luigi Stragapede]

Which is the meaning of the first 2 codes?

Q = @(v) sym(v);

V16 = @(v) vpa(v, 16);

[Walter Roberson]

I get a quite different plot

    OT = [0 0.780189980083088;3 0.755430990896907;6 0.731442157697117;9 0.710036424788487;12 0.684280102527316;15 0.653856319703993;18 0.619574167961039;21 0.58233456579523;24 0.543020606219544;27 0.502475470603902;30 0.461477022475297;33 0.420714176591535;36 0.380771555309333;39 0.342123546085244;42 0.305136241153018;45 0.270075032704127;48 0.237115735452024;51 0.206357475581531;54 0.177836008357966;57 0.151536524675646;60 0.127405347444595;63 0.105360190826173;66 0.0852988590170856;69 0.0671064032686978;72 0.050660846343487;75 0.0358376340058246;78 0.0225129942476625;81 0.0105663861101363;84 -0.000117791338377977;87 -0.00964907745802768;90 -0.01813028456402;93 -0.025657178496177;96 -0.0323184152464476;99 -0.0381956615657892;102 -0.0433638428993567;105 -0.0478914749239129;108 -0.0518410455471591;111 -0.0552694227342512;114 -0.0582282702468524;117 -0.0607644586082755;120 -0.0629204626163232;123 -0.0647347397530944;126 -0.0662420860928605;129 -0.0674739679552964;132 -0.0684588287306016;135 -0.0692223711266604;138 -0.0697878156443812;141 -0.070176136444506;144 -0.0704062759806561;147 -0.0704953398798443;150 -0.0704587735839264;153 -0.0703105222465887;156 -0.0700631753275887;159 -0.0697280972516447;162 -0.0693155454127023;165 -0.0688347767117329;168 -0.0682941437222885;171 -0.0677011814857774;174 -0.0670626858498591;177 -0.0663847841797106;180 -0.0656729991938642;183 -0.0649323066041588;186 -0.0641671871731083;189 -0.0633816737415182;192 -0.0625793937242093;195 -0.0617636075219076;198 -0.0609372432523592;201 -0.0601029281631467;204 -0.0592630170521526;207 -0.0584196179887523;210 -0.0575746155993032];
    omega = OT(:,1);
    Tval = OT(:,2);
    plot(omega, Tval);

The values decrease along a curve, becoming negative at roughly OMEGA=83