how to plot a function where there is inside 1x1 sym?
1 view (last 30 days)
Show older comments
Luigi Stragapede
on 12 May 2020
Commented: Walter Roberson
on 12 May 2020
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
2 Comments
Walter Roberson
on 12 May 2020
Are you using https://www.mathworks.com/matlabcentral/fileexchange/48403-bessel-zero-solver for besselzero() ?
Accepted Answer
Walter Roberson
on 12 May 2020
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
7 Comments
Walter Roberson
on 12 May 2020
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
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!