Error using plot: Data must be numeric, datetime, duration or an array convertible to double.

clc;
b = 0.142e-9; gammao = 3.0; m = 101;
hbar = 1; e = -1; 
K = 8.617e-16; T = 287.5;
a = ((3*b)/(2*hbar)); Pz = ((2*pi*hbar)/(3*b));
beta2 = 1; beta1 = linspace(0,10, 30); % However many you want.
Wcnzz = sqrt(3);
jo = ((8*e*Wcnzz*gammao)/(3*hbar*m*b));
%%
syms q s 
B1 = q.*beta1; B2 = q.*beta2;
Eqszz = (a./(2*pi)).*((1+(4.*cos(a.*Pz).*cos((pi.*s)./m))+(4.*(cos((pi.*s)./m)).^2)).^0.5);
Fqszz = ((a.^2).*m)./(((2.*(pi.^2).*s.*sqrt(3)).*(a./(2*pi)).*((1+(4.*cos(a.*Pz).*cos((pi.*s)./m))+(4.*(cos((pi.*s)./m)).^2)).^0.5))./(K.*T));
J1 = besselj(0,B1); J2 = besselj(0,B2);
J = q.*Fqszz.*Eqszz.*J1.*J2;
jz = symsum((symsum(J,s,1,m)),q,1,inf);
j = jz./jo;
plot(beta1, j, 'r-', 'LineWidth', 2 );
drawnow;
grid on;
fontSize = 20;
xlabel('T(K)', 'FontSize', fontSize)
ylabel('j_z/j_o', 'FontSize', fontSize)
hold on
%%
b = 0.142e-9; gammao = 3.0; m = 101;
hbar = 1; e = -1; 
K = 8.617e-16; T = 287.5;
a = ((3*b)/(2*hbar)); Pz = ((2*pi*hbar)/(3*b));
beta2 = 1; beta1 = linspace(0,10, 30); % However many you want.
Wcnac = 1; t = sqrt(3); n = 1e-9;
jo = ((8*e*Wcnac*gammao)/(3*hbar*m*b));
%%
syms q s 
B1 = q.*beta1; B2 = q.*beta2;
Eqsac = ((1+(4.*cos((pi.*s.*t)./n).*cos((a.*Pz)./t))+(4.*(cos((a.*Pz)./t)).^2)).^0.5);
Fqszz = ((a.^2).*n)./(((2.*(pi.^2).*s.*t).*((1+(4.*cos((a.*Pz)./t).*cos((a.*Pz)./t))+(4.*(cos((a.*Pz)./t)).^2)).^0.5))./(K.*T));
J1 = besselj(0,B1); J2 = besselj(0,B2);
J = q.*Fqszz.*Eqszz.*J1.*J2;
jz = symsum((symsum(J,s,1,m)),q,1,inf);
j = jz./jo;
plot(beta1, j, 'b-', 'LineWidth', 2);
drawnow;
grid on;
fontSize = 20;
xlabel('T(K)', 'FontSize', fontSize)
ylabel('j_z/j_o', 'FontSize', fontSize)
title('j_x/j_o vs.  \beta_1', 'FontSize', fontSize)
legend('zigzig CNs','armchair CNs','Location','Best');
% Maximize the figure window.
hFig.WindowState = 'maximized';

[madhan ravi]

You are using symbolic variables and use plot(...) :/ ?

[Samuel Suakye]

yeah

Answer 1

[Walter Roberson]

BesselJ oscillates infinitely often as it goes to infinity. For any given finite numeric q, there are an infinite number of zero crossings beyond q. It is like sin(x) in that way. So you cannot predict the sign of BesselJ(0,q) as q goes to infinity without a specific evaluation.

Can we estimate the magnitude of BesselJ(0,q) for large q? Because if we know the order of magnitude of BesselJ(0,q) then we can perhaps estimate the magnitude of q*BesselJ(0,q)

If you plot for some sample large q of the form 10^N, such as 10^100, you can see that the magnitude of BesselJ(0,q) ranges to approximately +/- 1/sqrt(q) -- not exactly by any means, but that order of magnitude.

So for q*Besselj(0,q), we have approximately +/- q*1/sqrt(q) . And as q approaches infinity, that approaches infinity.

Therefore, q*BesselJ(0,q) diverges as q approaches infinity.

At infinity, BesselJ(0,infinity) is 0. But infinity*0 is undefined, so there is no faint hope of somehow definining a result.

With q*besselj(0,q) not converging and oscillating infinitely often between negative and positive, the symsum() cannot be evaluated.

[Samuel Suakye]

so what do you think it is the best to do?

[Walter Roberson]

I think it is best to recheck that your implementation matches the equations. If it does, then the infinite sum cannot be found.

[Samuel Suakye]

Please is there any way apart from the symsum because it is very slow

[Walter Roberson]

Sure, you can just skip the calculation, knowing that it will produce an undefined result.

You can look at it yourself, your J2 is besselj(0,q) and your J multiplies it by q, and your first J1 is 1, so you have a q*besselj(0,q) being summed from 1 to infinity, which I demonstrated earlier does not converge.

More interesting is to start beta1 larger than 0. When you do that, then the first entry for J1 will be a besselj function, so you would have q times two besselJ functions. I already estimated that besselJ varies like 1/sqrt(q), so if you multiply two together you might get a variation like 1/(sqrt(q)^2) times q, which potentially gets you something like q/q -> 1 . And if you plot the function, it does oscillate but the peaks look to be fairly consistent. There would be room for optimism that just maybe the oscillations with consistent peak added up to a finite value.

If you try this, then one of the first steps is to try something like

syms M N q real
limit(q*besselj(0,M*q)*besselj(0,N*q),q,inf)

Unfortunately MATLAB cannot calculate that. But Maple can, and it says

-2/(sqrt(abs(M))*sqrt(abs(N))*Pi) .. 2/(sqrt(abs(M))*sqrt(abs(N))*Pi)

For positive M and N that simplifies to +/- 2/pi*sqrt(1/(M*N)) .

Unfortunately, when your infinite limit is not defined, your integral cannot be defined either.

Furthermore, the problem is structural: for non-zero beta1, none of your values can work, as they will all involve terms of the form symsum(q*besselj(0,M*q)*besselj(0,N*q),q,1,inf)

Fiddling with the other lines cannot help -- not unless you manage to introduce another besselj() or other factor that is less than linear in q.

[Samuel Suakye]

Thank you