how to plot a function where there is inside 1x1 sym?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Luigi Stragapede
le 12 Mai 2020
Commenté : Walter Roberson
le 12 Mai 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 commentaires
Walter Roberson
le 12 Mai 2020
Are you using https://www.mathworks.com/matlabcentral/fileexchange/48403-bessel-zero-solver for besselzero() ?
Réponse acceptée
Walter Roberson
le 12 Mai 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 commentaires
Walter Roberson
le 12 Mai 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
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Conversion Between Symbolic and Numeric dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!