带有求和的绘图函数会产生错误的结果

Question

我有一个方程式需要绘制，但绘制不正确。

等式如下：

情节应该是这样的：

但是我的代码：

clear; clc; close all;

eta      = 376.7303134617706554679; % 120pi
ka       = 4;
N        = 24;
coeff    = (2)/(pi*eta*ka);
Jz       = 0;

theta = [0;0.0351015938948580;0.0702031877897160;0.105304781684574;0.140406375579432;0.175507969474290;0.210609563369148;0.245711157264006;0.280812751158864;0.315914345053722;0.351015938948580;0.386117532843438;0.421219126738296;0.456320720633154;0.491422314528012;0.526523908422870;0.561625502317728;0.596727096212586;0.631828690107444;0.666930284002302;0.702031877897160;0.737133471792019;0.772235065686877;0.807336659581734;0.842438253476592;0.877539847371451;0.912641441266309;0.947743035161167;0.982844629056025;1.01794622295088;1.05304781684574;1.08814941074060;1.12325100463546;1.15835259853031;1.19345419242517;1.22855578632003;1.26365738021489;1.29875897410975;1.33386056800460;1.36896216189946;1.40406375579432;1.43916534968918;1.47426694358404;1.50936853747890;1.54447013137375;1.57957172526861;1.61467331916347;1.64977491305833;1.68487650695319;1.71997810084804;1.75507969474290;1.79018128863776;1.82528288253262;1.86038447642748;1.89548607032233;1.93058766421719;1.96568925811205;2.00079085200691;2.03589244590177;2.07099403979662;2.10609563369148;2.14119722758634;2.17629882148120;2.21140041537606;2.24650200927091;2.28160360316577;2.31670519706063;2.35180679095549;2.38690838485035;2.42200997874520;2.45711157264006;2.49221316653492;2.52731476042978;2.56241635432464;2.59751794821949;2.63261954211435;2.66772113600921;2.70282272990407;2.73792432379893;2.77302591769378;2.80812751158864;2.84322910548350;2.87833069937836;2.91343229327322;2.94853388716807;2.98363548106293;3.01873707495779;3.05383866885265;3.08894026274751;3.12404185664236;-3.12404185664236;-3.08894026274751;-3.05383866885265;-3.01873707495779;-2.98363548106293;-2.94853388716807;-2.91343229327322;-2.87833069937836;-2.84322910548350;-2.80812751158864;-2.77302591769378;-2.73792432379893;-2.70282272990407;-2.66772113600921;-2.63261954211435;-2.59751794821949;-2.56241635432464;-2.52731476042978;-2.49221316653492;-2.45711157264006;-2.42200997874520;-2.38690838485035;-2.35180679095549;-2.31670519706063;-2.28160360316577;-2.24650200927091;-2.21140041537605;-2.17629882148120;-2.14119722758634;-2.10609563369148;-2.07099403979662;-2.03589244590177;-2.00079085200691;-1.96568925811205;-1.93058766421719;-1.89548607032233;-1.86038447642748;-1.82528288253262;-1.79018128863776;-1.75507969474290;-1.71997810084804;-1.68487650695319;-1.64977491305833;-1.61467331916347;-1.57957172526861;-1.54447013137375;-1.50936853747890;-1.47426694358404;-1.43916534968918;-1.40406375579432;-1.36896216189946;-1.33386056800461;-1.29875897410975;-1.26365738021489;-1.22855578632003;-1.19345419242517;-1.15835259853032;-1.12325100463546;-1.08814941074060;-1.05304781684574;-1.01794622295088;-0.982844629056025;-0.947743035161167;-0.912641441266309;-0.877539847371451;-0.842438253476592;-0.807336659581735;-0.772235065686877;-0.737133471792019;-0.702031877897161;-0.666930284002303;-0.631828690107445;-0.596727096212586;-0.561625502317728;-0.526523908422871;-0.491422314528013;-0.456320720633154;-0.421219126738296;-0.386117532843439;-0.351015938948581;-0.315914345053722;-0.280812751158864;-0.245711157264007;-0.210609563369149;-0.175507969474290;-0.140406375579432;-0.105304781684575;-0.0702031877897167;-0.0351015938948580;-2.44929359829471e-16];

for n = 0:N
    if n == 0
        kappa  = 1;
    else
        kappa  = 2;
    end

    num    = (-1.^(n)).*(1i.^(n)).*(cos(n.*theta)).*(kappa);
    Hankel = besselh(n,2,ka);
    Jz     = Jz + ((num./Hankel));
end

Jz = Jz.*coeff;
x = linspace(0,2*pi,length(theta));
plot(x,abs(Jz));

生成以下不正确的图：

请注意，theta 的值是围绕圆柱体的离散角度。该方程是二维 TMz 极化圆柱体电流密度的解析解。

Answer 1

我认为您的结果实际上正确，这是一个简单的绘图问题或您如何指定 theta。由于这是一个周期函数，我们再画几个周期：

function q52693512

eta      = 376.7303134617706554679; % 120pi
ka       = 4;
N        = 24;
coeff    = (2)/(pi*eta*ka);
Jz       = 0;

theta = linspace(-3*pi, 3*pi, 180);
for n = 0:N
    kappa  = 1 + (n>0);
    num    = (-1.^(n)).*(1i.^(n)).*(cos(n.*theta)).*(kappa);
    Hankel = besselh(n,2,ka);
    Jz     = Jz + ((num./Hankel));
end

Jz = Jz.*coeff;
figure(); plot(theta, abs(Jz));

您可能已经能够看到所需的结果在那里但是相对于我们的结果偏移了半个周期 .如果我们再看一下中心会更清楚（如果忽略水平轴值，它正是您想要的形状）。

尝试寻找 φ 等于 theta ± π/2（或类似值）的理由。

带有求和的绘图函数会产生错误的结果

Plotting function with a summation produces a wrong result

matlab

plot

integral

trigonometry

bessel-functions