besselj: NU and Z must be the same size or one must be a scalar.

Plotting the following equation i.e. jx vs. omega (w) (find below) numerically but say NU and z of the besselj must have same size or one must be scalar

This is what I have done so far

n      = 5*10e12;
vf     = 1.0*10e6;
mo     = 9.109*10e-31;
m      = 0.44*mo;
tau    = 10e-12;
e      = 1.6022e-19;
E0     = 10; % Unit in KV/m
s      = 3.5*10e3;
hBar   = 6.52*10e-16;
Ed     = 50;
Eo     = 8.85*10-12;
B      = [0.1, 0.2, 0.3]; % Unit in Tesla
W      = linspace(0,10, 101); % However many you want.
k      = W./s;
a0     = (2*pi*(hBar^2)*Eo)/m/e^2;
s0     = ((e^2)*n*tau)/m;
jo     = 1./e/n/vf;

[L]   = meshgrid(-10000:200:10000);  % 1:0.1:3 span for 'q'
legendStrings = cell(length(B), 1);
for k1 = 1:length(B)
	thisN = B(k1);
	for i = 1:length(W)
        Wc = ((e.*thisN)./m);
        g      = 1 + (1i.*Wc.*tau); gstar = 1 - (1i.*Wc.*tau);
        BF     = (k.*vf)./Wc;
        R = (L.*(besselj(L,BF).^2))./(1-1i.*(W(i)-(L.*Wc)).*tau);
        Rx(i) = (Wc./k).*sum(R(:));
        Sx = ((L.*besselj(L,BF)).^2)./(1-1i.*(W(i)-(L.*Wc)).*tau);
        Sxx(i) = ((2*s0)./(BF.^2)).*sum(Sx(:));
        gkw(i) = 1 + (1i.*(1./Eo(Ed+1)).*(Sxx(i)./(s - Rx(i))));
        J = besselj(L,BF)./(1-1i.*(W(i)-(L.*Wc)).*tau).*(L+((k.*a0.*Sxx(i))./(Wc.*tau)./(Eo.*(s-Rx(i))))).*((g.*(L+1).*besselj(L+1,BF))+(gstar.*(L-1).*besselj(L-1,BF)));
        Jx(i) = jo.*((abs((s0*E0)./gkw(i))).^2).*(1./(BF.^2)./(1+(Wc.*tau).^2)).*sum(J(:));
  end

	legendStrings{k1} = sprintf('B = %.3f', thisN);
	plot(W, real(Jx), '-', 'LineWidth', 2, 'MarkerSize', 15);
	hold on;
	drawnow;
end

grid on;
fontSize = 15;
xlabel('\omega (ns^{-1})', 'FontSize', fontSize)
ylabel('\it j_{x} (\mu A/m)', 'FontSize', fontSize)
title('\it j_{x} (\mu A/m) vs. \omega (ns^{-1}) ', 'FontSize', fontSize)
legend(legendStrings, 'Location', 'Best')

Answer 1

[Walter Roberson]

Yes, your calculations are not being careful about whether you are operating on grids or vectors.

You can repmat(BF) to be a grid like L is, but you have more problems.

Your k is at least a vector but you have

        Rx(i) = (Wc./k).*sum(R(:));

With k being non-scalar, the right hand side is non-scalar.

Why are you saving to Rx(i), Sxx(i), gkw(i)? You are not using any of those after the loop and you never index at other elements within the loop.

Why are you calculating Wc, g, gstar, BF each time through the for i loop? They do not depend upon the value of i. Nothing depends upon the value of i until the denominator of R .

[Samuel Suakye]

Problem half - solved, after corrections and running I had this showing

Index exceeds array bounds.

Error in ben.m (line 37)
        gkw = 1+(1i.*(1./Eo(Ed+1)).*(Sxx./(s - Rx)));

[Walter Roberson]

MATLAB has no implied multiplication. A(B) is either function A called with parameter B, or else array A indexed at location B. Never in MATLAB is it ever A multiplied by B.

[Samuel Suakye]

thank you very much but the graph is not showing ....

n      = 5*10e12;
vf     = 1.0*10e6;
mo     = 9.109*10e-31;
m      = 0.44*mo;
tau    = 10e-12;
e      = 1; %1.6022e-19;
E0     = 10; % Unit in KV/m
s      = 3.5*10e3;
hBar   = 1; %6.52*10e-16;
Ed     = 50;
Eo     = 8.85*10-12;
B      = [0.1, 0.2, 0.3]; % Unit in Tesla
W      = linspace(0,10, 101); % However many you want.
k      = W./s;
a0     = (2*pi*(hBar^2)*Eo)/m/e^2;
s0     = ((e^2)*n*tau)/m;
jo     = 1./(e*n*vf);

[L]   = meshgrid(-10000:200:10000);  % 1:0.1:3 span for 'q'
legendStrings = cell(length(B), 1);
for k1 = 1:length(B)
	thisN = B(k1);
    Wc = ((e.*thisN)./m);
    g  = 1 + (1i.*Wc.*tau); gstar = 1 - (1i.*Wc.*tau);
    BF = repmat(((k.*vf)./Wc),101,1);
	for i = 1:length(W)
        R = (L.*(besselj(L,BF).^2))./(1-1i.*(W(i)-(L.*Wc)).*tau);
        Rx = ((Wc./k).*sum(R(:)));
        Sx = ((L.*besselj(L,BF)).^2)./(1-1i.*(W(i)-(L.*Wc)).*tau);
        Sxx = ((2*s0)./(BF.^2)).*sum(Sx(:));
        gkw = 1+(1i.*(1./(Eo*(Ed+1))).*(Sxx./(s - Rx)));
        J = besselj(L,BF)./(1-1i.*(W(i)-(L.*Wc)).*tau).*(L+((k.*a0.*Sxx)./(Wc.*tau)./(Eo.*(s-Rx)))).*((g.*(L+1).*besselj(L+1,BF))+(gstar.*(L-1).*besselj(L-1,BF)));
        Jx = jo.*((abs((s0*E0)./gkw)).^2).*(1./(BF.^2)./(1+(Wc.*tau).^2)).*sum(J(:));
  end

	legendStrings{k1} = sprintf('B = %.3f', thisN);
	plot(W, real(Jx), '-', 'LineWidth', 2, 'MarkerSize', 15);
	hold on;
	drawnow;
end

grid on;
fontSize = 15;
xlabel('\omega (ns^{-1})', 'FontSize', fontSize)
ylabel('\it j_{x} (\mu A/m)', 'FontSize', fontSize)
title('\it j_{x} (\mu A/m) vs. \omega (ns^{-1}) ', 'FontSize', fontSize)
legend(legendStrings, 'Location', 'Best')

[Samuel Suakye]

Please, is there any other way to plot the equations numerically because of the summations?

[Samuel Suakye]

also is there a way to use vectorization for numerical summution instead of the grid?

[Walter Roberson]

In

        Sx = ((L.*besselj(L,BF)).^2)./(1-1i.*(W(i)-(L.*Wc)).*tau);

The (L.*besselj(L,BF)) is coming out all zero. The besselj(L,BF) part is coming out all 0 except where L is 0; with BF ranging from 1e-37 to 1e-26 then the besselj(0, BF) is coming out as 1. But the L.* part is multiplying that and L is 0, so those 1's are being multiplied by 0. So one way or another (L.*besselj(L,BF)) is all zero, so Sx is all zero, so Sxx becomes all zero, and you divide by that zero...

[Samuel Suakye]

please is there any other way of solving that

[Walter Roberson]

        Jx = jo.*((abs((s0*E0)./gkw)).^2).*(1./(BF.^2)./(1+(Wc.*tau).^2)).*sum(J(:));
        

That line is in a for loop. You are not using Jx inside the for loop, so every iteration of the loop you are going to overwrite all of Jx. Aftr the loop you expect to plot Jx, but it is going to be a 2D matrix because it is built using BF which is a 2D matrix. There are cases where it is valid to use plot() with a 2D matrix: the effect is to draw one line per column, which in this case would be 101 columns. Are you sure that is what you want to do??

Please recheck your code.

[Walter Roberson]

The code below returns non-zero values for Jx.

Please pay attention to the for i loop, which I have altered to

for i = length(W):length(W)  %was 1: but all the variables are completely overwritten so no point running the rest

I did that because you completely overwrite all the output variables for each different i value, so the final result of the loop is the same as if you only did the last iteration, so to save quite a bit of wasted calculation, I chop out everything except that last iteration.

The Jx that is output is a 101 x 101 complex symbolic array. If you take double() of it, you will get all zero. The values that get stored in Jx vary in magnitude by roughly 10^300 depending on i, and roughly 10^40 depending on BF; expect values on the order of 1E-10800 to 1E-12000. The magnitude difference is large enough within any one array that you will not be able to do any useful contouring. You can, though, surf(double(log10(real(Jx(:,2:end))) ) to get something -- only a 10^45 range difference over the plot.

I would suggest to you that if you are expecting meaning values on the order of 1e-11050 then you need much higher precision in your input constants. Expecting 2 to 4 digits of precision of inputs to be able to reliably calculate down to below 1E-10000 is fundamentally Garbage In Garbage Out.

Q = @(v) sym(v);
n      = Q(5)*10e12;
vf     = Q(1.0)*10e6;
mo     = Q(9.109)*10e-31;
m      = Q(0.44)*mo;
tau    = Q(10e-12);
e      = Q(1); %1.6022e-19;
E0     = Q(10); % Unit in KV/m
s      = Q(3.5)*10e3;
hBar   = Q(1); %6.52*10e-16;
Ed     = Q(50);
Eo     = Q(8.85)*10-12;
B      = Q([0.1, 0.2, 0.3]); % Unit in Tesla
N = 101;
W      = linspace(Q(1e-10),Q(10), N); % However many you want.
k      = W./s;
a0     = (2*Q(pi)*(hBar^2)*Eo)/m/e^2;
s0     = ((e^2)*n*tau)/m;
jo     = 1./(e*n*vf);

[L]   = meshgrid(-Q(10000):Q(200):Q(10000));  % 1:0.1:3 span for 'q'
legendStrings = cell(length(B), 1);

digits(50);
Besselj = @(nu,z) vpa(besselj(sym(nu), sym(z)));

for k1 = 1:length(B)
	thisN = B(k1);
    Wc = ((e.*thisN)./m);
    g  = 1 + (1i.*Wc.*tau); gstar = 1 - (1i.*Wc.*tau);
    BF = repmat(((k.*vf)./Wc),N,1);
	for i = length(W):length(W)  %was 1: but all the variables are completely overwritten so no point running the rest
        R = (L.*(Besselj(L,BF).^2))./(1-1i.*(W(i)-(L.*Wc)).*tau);
        Rx = ((Wc./k).*sum(R(:)));
        Sx = ((L.*Besselj(L,BF)).^2)./(1-1i.*(W(i)-(L.*Wc)).*tau);
        Sxx = ((2*s0)./(BF.^2)).*sum(Sx(:));
        gkw = 1+(1i.*(1./(Eo*(Ed+1))).*(Sxx./(s - Rx)));
        J = Besselj(L,BF)./(1-1i.*(W(i)-(L.*Wc)).*tau).*(L+((k.*a0.*Sxx)./(Wc.*tau)./(Eo.*(s-Rx)))).*((g.*(L+1).*Besselj(L+1,BF))+(gstar.*(L-1).*Besselj(L-1,BF)));
        Jx = jo.*((abs((s0*E0)./gkw)).^2).*(1./(BF.^2)./(1+(Wc.*tau).^2)).*sum(J(:));
  end

	legendStrings{k1} = sprintf('B = %.3f', thisN);
	plot(W, real(Jx), '-', 'LineWidth', 2, 'MarkerSize', 15, 'DisplayName', legendStrings{k1});
	hold on;
	drawnow;
end

grid on;
fontSize = 15;
xlabel('\omega (ns^{-1})', 'FontSize', fontSize)
ylabel('\it j_{x} (\mu A/m)', 'FontSize', fontSize)
title('\it j_{x} (\mu A/m) vs. \omega (ns^{-1}) ', 'FontSize', fontSize)
legend(legendStrings, 'Location', 'Best')

[Samuel Suakye]

Thank you very much but I tried it on similar equations and still the graph is not showing

hFig = figure
clc;
% alignComments
Q = @(v) sym(v);
n      = Q(5)*1e12;
t      = Q(2.7);    % in eV
delta  = Q(0.26);    % in eV
mo     = Q(9.109)*1e-31;
m      = Q(0.44)*mo;
e      = Q(1); %1.6022e-19;
E0     = Q(10); % Unit in KV/m
s      = Q(3.5)*1e3;
hBar   = Q(1); %6.52e-16;
Ed     = Q(50);
Eo     = Q(8.85)*1e-12;
b      = Q(0.142)*1e-9;
a      = (3*b)/2/hBar;
p      = (2*Q(pi)*hBar)/3/b;
mu     = Q(1)*1e4;
cps    = cos(Q(pi)/sqrt(3));
cap    = cos(Q(pi));
sps    = sin(Q(pi)/sqrt(3));
sap    = sin(Q(pi));
Eqs    = sqrt(delta^2 + t^2.*(1 + 4*cps.*(cap + cps)));
vx      = -(2*a*(t^2)*cps*sap)/Eqs;
tau    = (mu*p)/e/vx;
tau1   = tau/2;
Sg     = ((e^2)*n*vx*tau)/p;
bg     = (2*Q(pi)*hBar*Eo*vx)/(e^2)/p;
U      = cps*cap - ((1/sqrt(3))*sap*sps);
B      = Q([0.1, 0.2, 0.3]); % Unit in Tesla
N = 101;
W      = linspace(Q(0),Q(10),N);%Q(0):Q(0.01):Q(10);  % However many you want.
k      = W./s;
jo     = -1./(2*e*n*vx);
SE     = (abs(Sg*E0))^2;
BB     = vx*p/4/m/U/(a^2)/(t^2);
R      = (-vx*Q(pi)*hBar)/3/b/(a^2)/(t^2)/U/m;
S      = -((vx^3)*(e^2)*tau)/3/b/(a^2)/(t^2)/U/hBar;

%L   = linspace(-10000,10000,201);  % 1:0.1:3 span for 'q'
[r] = meshgrid(-Q(10000):Q(200):Q(10000));

digits(50);
Besselj = @(nu,z) vpa(besselj(sym(nu), sym(z)));

legendStrings = cell(length(B), 1);
for k1 = 1:length(B)
	thisN = B(k1);
    Wc = ((e.*thisN.*vx)./p);
    g  = 1 + (1i.*Wc.*tau1); gstar = 1 - (1i.*Wc.*tau1);
    Bg = repmat(((k.*vx)./Wc),N,1);
	for i = length(W):length(W)
        Rx = R.*(Wc./k).*sum(r.*(Besselj(Bg,r).^2)./(1-(1i.*(W(i)-(r.*Wc)).*tau)));
        Sxx = S.*(Bg.^2).*sum(((r.*Besselj(Bg,r)).^2)./(1-(1i.*(W(i)-(r.*Wc)).*tau)));
        gkw = 1+(1i.*(1./(Eo*(Ed+1))).*(Sxx./(s-Rx)));
        CC = ((vx^2).*tau.*Wc.*(p^2))./(8*pi*(hBar^2)*(a^2)*(t^2)*U*n);
        DD = ((1./Bg)./(1 + (Wc.*tau1).^2)).^2;
        Jx = jo.*SE.*CC.*DD.*sum((Besselj(Bg,r)./(1-1i.*(W(i)-(r.*Wc)).*tau)).*(r-((k.*bg.*Sxx)./(Wc.*tau)./(Eo.*(s-Rx).*gkw))).*BB.*((-1i.*(g.^2).*(r+1).*Besselj(Bg,r+1))+(1i.*(gstar.^2).*(r-1).*Besselj(Bg,r-1))));
  end

	legendStrings{k1} = sprintf('B = %.3f', thisN);
	plot(W, real(Jx), '-', 'LineWidth', 2, 'MarkerSize', 15);
	hold on;
	drawnow;
end

grid on;
fontSize = 15;
xlabel('\omega (ns^{-1})', 'FontSize', fontSize)
ylabel('\it j_{x} (\mu A/m)', 'FontSize', fontSize)
title('\it j_{x} (\mu A/m) vs. \omega (ns^{-1}) ', 'FontSize', fontSize)
%legend(legendStrings, 'Location', 'Best')
legend('B = 0.1T','B = 0.2T','B = 0.3T','Location','Best')
% Maximize the figure window.
hFig.WindowState = 'maximized';

[Walter Roberson]

I have not compared line by line; visually the code looks quite similar. I would expect that it is going to produce similar magnitude of results, Jx values on the order of 1e-10000. plot() of values converts them to floating-point and the floating-point limit is about 1e-324, many orders of magnitude higher. So all of the Jx values get converted to 0.

You could try plotting log10(Jx)

[Samuel Suakye]

Same occurring problem

[Walter Roberson]

sap    = sin(Q(pi));

Exactly 0.

vx      = -(2*a*(t^2)*cps*sap)/Eqs;

Multiplies by the exact 0, so the result is exact 0.

tau    = (mu*p)/e/vx;

Divide by 0.

jo     = -1./(2*e*n*vx);

Divide by 0.

Jx = jo.*SE.*CC.*DD.*sum((Besselj(Bg,r)./(1-1i.*(W(i)-(r.*Wc)).*tau)).*(r-((k.*bg.*Sxx)./(Wc.*tau)./(Eo.*(s-Rx).*gkw))).*BB.*((-1i.*(g.^2).*(r+1).*Besselj(Bg,r+1))+(1i.*(gstar.^2).*(r-1).*Besselj(Bg,r-1))));

Multiplies by the result of division by 0, so without looking further you can predict that Jx will either be +inf or -inf or NaN. It turns out that it becomes NaN.