ZIP with MATLAB scripts and note:
exercise 6.18
exercise 6.18 notes:
eps0 = 8.854187817e-12 % [F/m]=[C^2/N*m^2] permittivity in empty space
mu0 = 4*pi*1e-7 % [H/m]=[T*m/A] permeability in empty space
c0 = 1/sqrt(mu0*eps0) % light speed through empty space [m/s]
etha0 = (mu0/eps0)^.5 % ~ 120*pi [ohm] characteristic impedance of empty space
mur=1;er=36.2
a=(7.99/2)*10^-3 % a: cylinder radius
L0=2.14e-3 % L0: cylinder height
N=5;
M=5;
dx=.001; % dx spatial resolution to find J(x) J'(x) roots
1.- this exercise shows that although sometimes coincide, cut-off frequencies not always have same values as resonant frequencies:
The support function intersections.m needs update.
The earlier version has a limitation that for dx 'too' small, for instance dx=.0001 intersections.m attempts to borrow too much RAM exceeding any reasonable RAM request.
Since I have detected it while solving this exercise 6.18 I mention such upgrade in this note, besides the intersections.m used here, also including intersections_OBSOLETE.m in case interested readers would like to have a look at Mr Schwarz's function upgrade that now doesn't crash because of attempting to borrow too much RAM, now instead borrowing time, leaving to the user the decision to whether wait, use more computing power, or find an alternative way.
[p p1]=roots_bessel(dx,N,M)
p =
Columns 1 through 3
2.404825587623834 5.520078116818740 8.653727924339440
3.831705997281383 7.015586687090941 10.173468147302287
5.135622324716259 8.417244151376378 11.619841177883229
6.380161906569109 9.761023131142313 13.015200727876815
7.588342449349022 11.064709497802784 14.372536680265151
Columns 4 through 5
11.791534449561969 14.930917711006265
13.323691944300002 16.470630057943414
14.795951783895344 17.959819499090457
16.223466167991070 19.409415232696745
17.615966050730911 20.826932958455323
>> p1
p1 =
Columns 1 through 3
3.832206002438044 7.016086681407118 10.173968140669018
1.841683794603466 5.331942785453609 8.536816379612862
3.054736905459390 6.706633213085026 9.969967828342625
4.201688902908669 8.015736610027316 11.346424323015336
5.318053105903168 9.282896290396309 12.682408452626531
Columns 4 through 5
13.324191945277404 16.471130056856072
11.706504916619080 14.864088639390282
13.170870863084952 16.348022321444045
14.586348295876196 17.789247873220077
15.964607046454198 19.196528807921194
saving J(x) J'(x) roots for next time in case dx too small no need to run again function roots_bessel
save p_roots_J_bessel.mat p
save p1_roots_J1_bessel.mat p1
.
If we start with the following expressions
k0=2*pi/lambda0 % wave number without dielectric
k_TE=2*pi*f0*(mu0*e0*mur*er)^.5 % wave number with dielectric
k_TM=2*pi*f0*(mu0*e0*mur*er)^.5
beta_TE_nm=(k_TE^2.-(p1/a).^2).^.5 % is there propagation?
beta_TM_nm=(k_TM^2.-(p/a).^2).^.5
the unsorted cut-off frequencies would be
L=repmat([1:M],5,1)
fc_TE_nmL=c0/(2*pi*(mur*er)^.5)*((p1/a).^2+(L*pi/L0).^2).^.5
fc_TM_nmL=c0/(2*pi*(mur*er)^.5)*((p/a).^2+(L*pi/L0).^2).^.5
cut-off frequencies, not resonating frequencies
L=repmat([1:M],5,1)
fc_TE_nmL=c0/(2*pi*(mur*er)^.5)*((p1/a).^2+(L*pi/L0).^2).^.5
fc_TM_nmL=c0/(2*pi*(mur*er)^.5)*((p/a).^2+(L*pi/L0).^2).^.5
2.- to know at what frequency an er=3.5 dielectric cyl with radius a=3.995 %mm and L0=2.14 %mm tall, without metallic walls we need to following expressions
alpha=((pr/a)^2-(k0_)^2)^.5
beta=(er*(k0_)^2-(pr/a)^2)^.5
tan(L0/2*(er*(k0)^2-(pr/a)^2)^.5)-((pr/a)^2-(k0)^2)^.5/(er*(k0)^2-(pr/a)^2)^.5==0
assuming L0 << lambda (pg294), which may not be the case.
[POZAR] also considers TE01x to be the lowest and only resonating mode to be taken into account to then solve (6.70)
1st try with command solve:
syms k0
pr=p(1,1);
eq1=tan(L0/2*(er*(k0)^2-(pr/a)^2)^.5)-((pr/a)^2- ...
(k0)^2)^.5/(er*(k0)^2-(pr/a)^2)^.5==0
k0_=double(solve(eq1,k0))
alpha=((pr/a)^2-(k0_)^2)^.5
beta=(er*(k0_)^2-(pr/a)^2)^.5
a =
5.830608049524744e+02
beta =
6.695731738706555e+02
pr =
3.832206002438044
eq1 =
tan((107*((181*k0^2)/5 - 7904128337556299/8589934592)^(1/2))/100000) - (7904128337556299/8589934592 - k0^2)^(1/2)/((181*k0^2)/5 - 7904128337556299/8589934592)^(1/2) == 0
Warning: Unable to solve symbolically. Returning a numeric solution
using vpasolve.
> In solve (line 304)
command solve cannot solve symbolically, it reverts to numerical solution 2nd try:
pr=p1(1,1)
eq1=tan(L0/2*(er*(k0)^2-(pr/a)^2)^.5)- ...
((pr/a)^2-(k0)^2)^.5/(er*(k0)^2-(pr/a)^2)^.5==0
k0_=double(solve(eq1,k0))
alpha=((pr/a)^2-(k0_)^2)^.5
beta=(er*(k0_)^2-(pr/a)^2)^.5
k0_ =
-2.080945166134460e+02
alpha =
9.364071317185853e+02
beta =
8.046234049516080e+02
3.- but is the resonant frequency alpha dependent?
I know that for x small y<0, but only for some values of alpha, as trying some different values of
y sometimes starts y>0 but for some other alphas, x small enough then y<0.
Therefore I am going to gradually increase x span until y shows sign change,
and there is going to be the 1st zero (from left) of equation tan(beta*L0/2)-alpha/beta = = 0
dx=.0001;x=[dx:dx:3];alpha_1=.01;y=tan(x*L0/2)-alpha_1./x;
y0=find(y>0)
n1=0
while isempty(y0)
n1=n1+dx;
dx=.0001;x=[dx:dx:3+n1];alpha_1=.01;y=tan(x*L0/2)-alpha_1./x;
y0=find(y>0);
end
y0
x0=3+n1
y0 =
30571
x0 =
3.057100000000001
this resonant frequency is way below the expected 7.11GHz and measured 7.8GHz yet not a single clue supplied about how such measurement has been done.
f0=x0/(2*pi)*c0
y0=[];n0=[];f0=[];
x1=.1
dx=.001;x=[dx:dx:x1];
% let's sweep alpha 10 fold from 1e-5 to 1
figure;
s1_range=10.^[-5:1:3];
for s1=s1_range
alpha_1=s1;
y=tan(x*L0/2)-alpha_1./x;
if y(end)>0 y0_=find(y<0); end
if y(end)<0 y0_=find(y>0); end
n0_=0
while isempty(y0_)
n0_=n0_+dx;
dx=.0001;
x=[dx:dx:3+n0_];
alpha_1=s1;
y=tan(x*L0/2)-alpha_1./x;
if y(end)>0 y0_=find(y<0); end
if y(end)<0 y0_=find(y>0); end
end
y0=[y0 y0_] % update meters
n0=[n0 n0_]
plot(x,10*log10(abs(y)));
hold on
end
grid on
showing 10*log10(abs(y)) rather than just y, y alone with quite sharp corner really close to both axes
x0=x1+n0;
f0=x0/(2*pi)*c0
f0 =
1.0e+09 *
Columns 1 through 3
0.047713451594041 0.147911699941527 0.467591825621601
Columns 4 through 5
1.460031618777663 4.613890769143707
that when applying the dielectric relative permittivity to wavelength:
f0*er^.5
f0
ans =
1.0e+10 *
Columns 1 through 3
0.028707483233682 0.088993198024413 0.281333335690080
Columns 4 through 5
0.878448986950665 2.776013628696985
shows one of the resonant frequencies to be quite close to the measurement, not to the trial-error poking suggested in the solutions manual
f0
f0 =
1.0e+09 *
Columns 1 through 3
0.047713451594041 0.147911699941527 0.467591825621601
Columns 4 through 5
1.460031618777663 4.613890769143707
for the record, if some one claims a resulting figure out of RF measurements, the basic thing to do is to give readers also a basic layout of such measurement, a photo, dimensions of the test board, circuit schematic, test instruments basic details .. in similar way that whoever wrote the solutions manuals quotes here and there that this or that electromagnetic simulator has been used.
% generating tag for command legend
L10=s1_range;
tg1={};
for s1=1:1:numel(f0)
tg1=[tg1 ['\alpha=' num2str(L10(s1))]];
end
legend(tg1)
title('y=tan(x*L0/2)-\alpha/x');
xlabel('x = k0 [1/m]');ylabel('10*log10(|y|)');
4.- tan(b*L0/2)=alpha/b
let's plot the surface that the solutions manual seems to use to find 7.11GHz resonance at alpha=583 and beta=664.
a1=.01;a2=100
b1=20;b2=1700 % beta range
Nmesh=1e2 % same amount of steps
alpha_range=linspace(a1,a2,Nmesh);
beta_range=linspace(b1,b2,Nmesh);
[ALF,BET]=meshgrid(alpha_range,beta_range);
% build 2D plane points
Y=tan(BET*L0/2)-ALF./BET;
% calculate eq (6.70)
figure;hm1=surf(ALF,BET,Y)
grid on
hm1.EdgeColor='none'
xlabel('\alpha');ylabel('\beta');zlabel('tan(\beta*L0/2)- ...
\alpha/\beta');
title(['y=tan(\beta*L0/2)- ...
\alpha/\beta \alpha:[' num2str(a1) ' , ...
' num2str(a2) '] \beta:[' num2str(b1) ' , ...
' num2str(b2) ']']);
min(min(Y))
I have used alpha values way too small, this dielectric material has to have heavy loss
campos(1.0e+02 *[-8.159611968576154 ...
7.108311333195441 ...
0.526126500110708])
one can see that effectively there's a zero with beta somewhere between [400 1200] with alpha constant.
a1=300;a2=700
b1=200;b2=1000 % beta range
Nmesh=1e3 % same amount of steps
alpha_range=linspace(a1,a2,Nmesh);
beta_range=linspace(b1,b2,Nmesh);
[ALF,BET]=meshgrid(alpha_range,beta_range);
% build 2D plane points
Y=tan(BET*L0/2)-ALF./BET;
% calculate eq (6.70)
figure;hm1=surf(ALF,BET,Y)
grid on
hm1.EdgeColor='none'
xlabel('\alpha');ylabel('\beta');
zlabel('tan(\beta*L0/2)-\alpha/\beta');
title(['y=tan(\beta*L0/2)-...
\alpha/\beta \alpha:[' num2str(a1) ' ,...
' num2str(a2) '] \beta:[' num2str(b1) ' ,...
' num2str(b2) ']']);
min(min(Y))
hold on
log_min=[0 0 0];
for k1=1:1:size(Y,1)
slice_Y=abs(Y(k1,:));
[slice_min,slice_loc]=min(slice_Y);
if slice_min<.02
log_min=[log_min ; ...
[alpha_range(k1) ...
beta_range(slice_loc) slice_min]];
plot3(alpha_range(slice_loc), ...
beta_range(k1),Y(k1,slice_loc), ...
'ro','MarkerFaceColor','r');
end
end
log_min(1,:)=[];
size(log_min)
alpha_min=log_min(:,1);
beta_min=log_min(:,2);
So assuming just TE01z (z small) then
k0_from_alpha=(alpha_min.^2-(p1_(1,1)/a)^2).^.5;
% p1_(1,1) is p'(01) 1st (smallest) root of J'(x) Bessel 1st kind 1st order.
k0_from_beta=(1/er*(beta_min.^2+(p1_(1,1)/a)^2)).^.5;
f0_ov_alpha=k0_from_alpha*c0/(2*pi);
f0_ov_beta=k0_from_beta*c0/(2*pi);
figure;
plot(alpha_min,f0_ov_alpha/1e9)
title('f0(\alpha_m_i_n) [GHz]')
grid on;xlabel('\alpha_m_i_n [dB/m]');ylabel('f0 [GHz]')
figure;
plot(beta_min,f0_ov_beta)
title('f0(\beta_min) [GHz]')
grid on;title('f0(\beta_m_i_n) [GHz]')
grid on;xlabel('\beta_m_i_n');ylabel('f0 [GHz]')
5.-
When looking for lowest resonant frequency of a device, the key physical dimensions of the device have to be present somewhere in the equations to solve. So again, the more or less general transcendental equation with alpha beta and L0 cylinder tall is valid point to start understanding dielectric resonators, but in my opinion, the equation to work on is:
tan(L0/2*(er*(k0)^2-(pr/a)^2)^.5)-((pr/a)^2-(k0)^2)^.5/(er*(k0)^2-(pr/a)^2)^.5
As shown in point 1.- command solve, at least for this exercise, doesn't work.
So from example 6.5 and taking advantage that alpha and beta have to be both real let's check, not just p(1,1) = p_01 = 2.40482448 ~ 2.405, let's find J(x) and J'(x) 12th lowest roots:
M0=3;N0=4;
% building matrix with roots and indices on each row:
p_=[0 0 0];p1_=[0 0 0];
for s1=1:1:M0
for s2=1:1:N0
p_=[p_ ; p(s1,s2) s2-1 s1];
p1_=[p1_; p1(s1,s2) s2-1 s1];
end
end
p_(1,:)=[];p1_(1,:)=[];
p_
p1_
sortrows(p_) % TM modes roots
sortrows(p1_) % TE roots
p_=[p_ zeros(M0*N0,1)];
p1_=[p1_ ones(M0*N0,1)];
pr_=[p_;p1_];
an extra column has been added, 0: TM mode, 1: TE mode, and sorting TE TM modes in single matrix.
for this device TE TM modes are interleaved:
pr_ =
Columns 1 through 3
2.404825587623834 0 1.000000000000000
5.520078116818740 1.000000000000000 1.000000000000000
8.653727924339440 2.000000000000000 1.000000000000000
11.791534449561969 3.000000000000000 1.000000000000000
3.831705997281383 0 2.000000000000000
7.015586687090941 1.000000000000000 2.000000000000000
10.173468147302287 2.000000000000000 2.000000000000000
13.323691944300002 3.000000000000000 2.000000000000000
5.135622324716259 0 3.000000000000000
8.417244151376378 1.000000000000000 3.000000000000000
11.619841177883229 2.000000000000000 3.000000000000000
14.795951783895344 3.000000000000000 3.000000000000000
3.832206002438044 0 1.000000000000000
7.016086681407118 1.000000000000000 1.000000000000000
10.173968140669018 2.000000000000000 1.000000000000000
13.324191945277404 3.000000000000000 1.000000000000000
1.841683794603466 0 2.000000000000000
5.331942785453609 1.000000000000000 2.000000000000000
8.536816379612862 2.000000000000000 2.000000000000000
11.706504916619080 3.000000000000000 2.000000000000000
3.054736905459390 0 3.000000000000000
6.706633213085026 1.000000000000000 3.000000000000000
9.969967828342625 2.000000000000000 3.000000000000000
13.170870863084952 3.000000000000000 3.000000000000000
Column 4
0
0
0
0
0
0
0
0
0
0
0
0
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
pr_TE_TM=sortrows(pr_)
pr_TE_TM =
Columns 1 through 3
1.841683794603466 0 2.000000000000000
2.404825587623834 0 1.000000000000000
3.054736905459390 0 3.000000000000000
3.831705997281383 0 2.000000000000000
3.832206002438044 0 1.000000000000000
5.135622324716259 0 3.000000000000000
5.331942785453609 1.000000000000000 2.000000000000000
5.520078116818740 1.000000000000000 1.000000000000000
6.706633213085026 1.000000000000000 3.000000000000000
7.015586687090941 1.000000000000000 2.000000000000000
7.016086681407118 1.000000000000000 1.000000000000000
8.417244151376378 1.000000000000000 3.000000000000000
8.536816379612862 2.000000000000000 2.000000000000000
8.653727924339440 2.000000000000000 1.000000000000000
9.969967828342625 2.000000000000000 3.000000000000000
10.173468147302287 2.000000000000000 2.000000000000000
10.173968140669018 2.000000000000000 1.000000000000000
11.619841177883229 2.000000000000000 3.000000000000000
11.706504916619080 3.000000000000000 2.000000000000000
11.791534449561969 3.000000000000000 1.000000000000000
13.170870863084952 3.000000000000000 3.000000000000000
13.323691944300002 3.000000000000000 2.000000000000000
13.324191945277404 3.000000000000000 1.000000000000000
14.795951783895344 3.000000000000000 3.000000000000000
Column 4
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
1.000000000000000
0
each root in turn, for alpha and beta to be both real, generates lower and upper frequency limits.
BW=[pr_TE_TM(:,1)*c0/(2*pi*a*er^.5) pr_TE_TM(:,1)*c0/(2*pi*a)]
BW =
1.0e+11 *
0.036558202349153 0.219957673555300
0.047736805147761 0.287215342345935
0.060637819717177 0.364836148863998
0.076060984845251 0.457631836355149
0.076070910159111 0.457691553431091
0.101944275497212 0.613362370952755
0.105841319685767 0.636809496860918
0.109575885594000 0.659279048866473
0.133129506163845 0.800992788914944
0.139262363309180 0.837892004501407
0.139272288407853 0.837951720282639
0.167085856871810 1.005296048490417
0.169459416180969 1.019576908875781
0.171780159796996 1.033539996060002
0.197908078656671 1.190742371393672
0.201947646067996 1.215047008735484
0.201957571147821 1.215106724403320
0.230658762536131 1.387791563382668
0.232379074408518 1.398142066769731
0.234066946604651 1.408297392186838
0.261447357867914 1.573035567902801
0.264480920972854 1.591287435890544
0.264490846203755 1.591347152467345
0.293705901549961 1.767123727731302
checking 1st root with Newton-Raphson
fun1=@(f0_) tan(L0/2*(er*(2*pi*f0_/c0).^2-...
(pr_TE_TM(1,1)/a).^2).^.5)-((pr_TE_TM(1,1)/a).^2-...
(2*pi*f0_/c0).^2).^.5./(er*(2*pi*f0_/c0).^2-...
(pr_TE_TM(1,1)/a).^2).^.5
x0_f0_ = [BW(1,1)+1e6];
% newtonraphson initial guess
options1 = optimset('TolX',1e-12);
[f0_newtonraphson, resnorm, f3, exitflag, output, jacob] =... newtonraphson(fun1, x0_f0_, options1);
Niter resnorm stepnorm lambda rcond convergence
-------------------------------------------------------------------
0 42.14 0 1 1 Inf
1 24.42 1.974e+06 1 1 0.5456
2 14.33 5.634e+06 1 1 0.5334
3 8.432 1.596e+07 1 1 0.5301
4 4.788 4.901e+07 1 1 0.5659
5 2.678 1.383e+08 1 1 0.5809
6 1.428 3.429e+08 1 1 0.6289
7 0.5873 7.348e+08 1 1 0.8885
8 0.09362 8.43e+08 1 1 1.836
9 0.01317 2.22e+08 1 1 1.961
10 0.001244 3.057e+07 1 1 2.36
11 1.026e-05 2.609e+06 1 1 4.798
12 5.357e-07 2.284e+04 1 1 2.952
6.- Now let's find the lowest 12 resonant frequencies, from roots in pr_, with newtonraphson
% f0 format:
% [resonant_freq N_modeindex M_modeindex mode]
% mode being 1: TE, 0: TM
options11 = optimset('TolX',1e-12);
f0_nr=[0 0 0 0];
for s4=1:1:N0*M0
x0_f0_ = [BW(s4,1)+1e6];
[f0_newtonraphson, resnorm, f, exitflag, output, jacob] = ...
newtonraphson(fun1, x0_f0_, options11);
f0_nr=[f0_nr; f0_newtonraphson pr_TE_TM(s4,2) pr_TE_TM(s4,4) pr_TE_TM(s4,4)]
end
f0_nr(1,:)=[];
f0_nr(:,1)
f0_nr =
1.0e+10 *
Column 1
0.000000000000000 + 0.000000000000000i
0.598246503527381 + 0.000000000000000i
0.598246412363155 + 0.000000000000000i
0.598246365019343 + 0.000000000000000i
0.598246519010487 + 0.000000000000000i
0.598246363143861 + 0.000000000000000i
0.598246418812999 + 0.000000000000000i
0.598246357760106 + 0.000000000000000i
0.598246394545926 + 0.000000000000000i
2.338546288922569 + 0.000000000005164i
2.347149417425641 - 0.000000000007259i
2.304960174314115 + 0.000000000003941i
2.314852331070381 + 0.000000000004785i
Column 2
0.000000000000000 + 0.000000000000000i
0.000000000000000 + 0.000000000000000i
..
0.000000000100000 + 0.000000000000000i
0.000000000100000 + 0.000000000000000i
Column 3
0.000000000000000 + 0.000000000000000i
0.000000000100000 + 0.000000000000000i
..
0.000000000100000 + 0.000000000000000i
0.000000000000000 + 0.000000000000000i
Column 4
0.000000000000000 + 0.000000000000000i
0.000000000100000 + 0.000000000000000i
..
0.000000000100000 + 0.000000000000000i
0.000000000000000 + 0.000000000000000i
ans =
1.0e+10 *
0.598246503527381 + 0.000000000000000i
0.598246412363155 + 0.000000000000000i
0.598246365019343 + 0.000000000000000i
0.598246519010487 + 0.000000000000000i
0.598246363143861 + 0.000000000000000i
0.598246418812999 + 0.000000000000000i
0.598246357760106 + 0.000000000000000i
0.598246394545926 + 0.000000000000000i
2.338546288922569 + 0.000000000005164i
2.347149417425641 - 0.000000000007259i
2.304960174314115 + 0.000000000003941i
2.314852331070381 + 0.000000000004785i
a bit better with TolX up to .1Hz resolution
options12 = optimset('TolX',.1);
f0_nr=[0 0 0 0];
for s4=1:1:N0*M0
x0_f0_ = [BW(s4,1)+1e6];
[f0_newtonraphson, resnorm, f, exitflag, output, jacob] = ...
newtonraphson(fun1, x0_f0_, options12);
f0_nr=[f0_nr; f0_newtonraphson pr_TE_TM(s4,2) pr_TE_TM(s4,4) pr_TE_TM(s4,4)]
end
f0_nr(1,:)=[];
f0_nr(:,1)
saving NR results to text file:
str_pr_TE_TM=num2str(pr_TE_TM([1:N0*M0],4));
% build string for fprintf
str_pr_TE_TM=strrep(str_pr_TE_TM','1','E')'
str_TEM=strrep(str_pr_TE_TM','0','M')'
str_f0_nr= ...
[repmat('f0 T',N0*M0,1) ...
str_TEM ...
repmat(' ',N0*M0,1) ...
num2str(pr_TE_TM([1:N0*M0],2)) ...
repmat(' ',N0*M0,1) ...
num2str(pr_TE_TM([1:N0*M0],3)) ...
repmat(' = ',N0*M0,1) ...
num2str(round(floor(f0_nr(:,1))/1e9,4)) ...
repmat(' GHz',N0*M0,1)]
file_id=fopen('f0_nr.txt','w'); % writing to text file
for s5=1:1:size(str_f0_nr,1) fprintf(file_id,'%s\n',str_f0_nr(s5,:)); end
fclose(file_id)
f0_nr =
1.0e+10 *
Columns 1 through 3
0 0 0
0.365682023491535 0 0.000000000100000
0.564989820779805 0 0
0.606478197171770 0 0.000000000100000
0.612705763558529 0 0
0.612707491116259 0 0.000000000100000
0.661222918903519 0 0
0.648448958477496 0.000000000100000 0.000000000100000
0.609253255839162 0.000000000100000 0
2.278229515968897 0.000000000100000 0.000000000100000
2.290302142223282 0.000000000100000 0
2.304960174310541 0.000000000100000 0.000000000100000
2.131478925028222 0.000000000100000 0
Column 4
0
0.000000000100000
0
0.000000000100000
0
0.000000000100000
0
0.000000000100000
0
0.000000000100000
0
0.000000000100000
0
f0 TE 0 2 = 3.6568 GHz
f0 TM 0 1 = 5.6499 GHz
f0 TE 0 3 = 6.0648 GHz
f0 TM 0 2 = 6.1271 GHz
f0 TE 0 1 = 6.1271 GHz
f0 TM 0 3 = 6.6122 GHz
f0 TE 1 2 = 6.4845 GHz
f0 TM 1 1 = 6.0925 GHz
f0 TE 1 3 = 22.7823 GHz
f0 TM 1 2 = 22.903 GHz
f0 TE 1 1 = 23.0496 GHz
f0 TM 1 3 = 21.3148 GHz
Newton-Raphson is known to be faster but to require the existence of 1st derivative.
Rewording; functions with completely vertical slopes, always saying this depending on how the MATLAB function has been implemented, without having had a look at the function itself, may not render expected results.
7.-
In pages 139 and 140 [POZAR] mentions the bisection method, the classic entry point to numerical methods for roots finding, a method that intuitively can solve more functions that Newton-Raphson.
options10 = optimset('TolX', .1);
% frequency resolution 0.1Hz ok
target1=0;
f0_bisect=bisection(fun1,BW(1,1),BW(1,2),target1,options10)
f0_bisect =
1.220238084617503e+10
repeating 12 lowest roots, now with bisection, and reducing (improving) tolerance step size from 0.1 down to 1e-4, the resulting resonant frequencies are the same, now with bisection:
options10 = optimset('TolX', 1e-5);
f0_bsc=[0 0 0 0];
for s4=1:1:N0*M0
x0_f0_ = [BW(s4,1)+1e6];
f0_bisect=bisection(fun1,BW(s4,1),BW(s4,2),target1,options10)
f0_bsc=[f0_bsc; f0_bisect pr_TE_TM(s4,2) pr_TE_TM(s4,4) pr_TE_TM(s4,4)]
end
f0_bsc(1,:)=[];
f0_bsc(:,1)
ans =
1.0e+11 *
0.122023808461750
0.287215342345064
0.364836148863445
0.457631836354454
0.457691553430397
0.613362370951825
0.636809496859953
0.659279048865473
0.800992788914336
NaN
NaN
1.005296048489654
but if resolution step below 1e-5, i.e. 1e-5 the loop trapped endless iterations, stopped manually.
8.-
Function fzero, for same inputs as newtonraphson, gets lost and returns NaN
fzero(fun1,x0_f0_, options1)
fzero(fun1,x0_f0_, options1)
=
1.220238084620664e+10
9.- COMMENT 1: Comparing Mathworks' fzero, Sky's bisection, and Mikofski's newtonraphson
It's really what is commonly known about don't buy a book by the cover, or when it comes to rifles and cars, no matter what the leaflet reads, because it's unlikely that you are going to buy exactly the one used in the leaflet review, always try before buy.
Numerical methods literature usually relegates the bisection method to some kind of early slow method, then spending loads of time on more complicated methods,
Bringing tolerance for bisection until loop trapped was a bit of test.
When comparing these 3 numerical methods and playing a bit with tolerances, then not even giving advantage of a coarse tolerance to newtonraphson and fzero by 5 orders of magnitude, this is, forcing bisection to refine down to 1e-5 while allowing newtonraphson and fzero using broader steps, broadest 1e-5, no even then Sky's implementation of the Bisection method, it doesn't fall behind.
options2 = optimset('TolX', 1e-9);
options3 = optimset('TolX', 1e-5);
target = [(-100:.1:100)' (-1000:1:1000)'];
tic;
xfz = zeros(size(target));
for ii = 1:numel(target)
xfz(ii) = fzero(@(x) x.^3-target(ii), [-20 20], options3);
end
fzero_time = toc
tic;
xbis = bisection(@(x) x.^3, -20, 20, target, options3);
bisection_time = toc
tic
[x, resnorm, f, exitflag, output, jacob] = ...
newtonraphson(@(x) x.^3, -20, options3);
nr_time = toc
fzero_time = 0.759273200000000
bisection_time = 0.014720300000000
nr_time = 0.017901500000000
It seems as if the solutions manual has exercise 6.18 and example 6.5 pg297 crossed: Example 6.15 has a clear definition of the attenuation, single value, therefore it is possible to apply equation (6.70), same as (6.65a), and obtain resonant frequency from k0. But in exercise 6.18 no value for the attenuation alpha is mentioned in the question.
Unknown attenuation means we cannot get the root value from in the way
shown in the solutions manual and example 6.5, this would be, straight to p01=2.405 (L0<<lambda assumed) as if alpha had been measured, known, and then the root can be calculated.
However without mentioning a specific TE TM mode, or band, or alpha measurement, the start point should be the transcendental equation without assuming any Bessel root known, for then sweeping root values, building one transcendental equation for each root, and solving each equation.
10.- APPLICATION
To close this solution I would like to mention the cylindrical resonator design by [CUEN]:
Although resonators may be key parts of RF microwave and optical devices, it's just one part that is relatively cheap comparatively speaking to the overall cost of for instance a spectrum analyser. And the RF microwave system is usually in turn transporting, acquiring or/and serving data to other systems, part of a major project.
[CUEN] shows a cylindrical resonating test system to characterise doped dielectric to design and manufacture biomedical products.
to obtain complex relative permittivity measurements the material performance while tested inside the cavity that because of the sandy nature of material, the resonating cavity requires metallic walls therefore limiting the possible TEM modes to TM modes only, L0 = 40mm as shown above.
The purpose of [CUEN] is to measure permittivity, but regarding cylindrical cavity resonant frequencies, with the above lay-out [CUEN] acknowledges the resonant frequencies, simulated and measured, shown in these 2 tables:
REFERENCES:
newtonraphson.m
NewtonRaphson, by Mark Mikofski. v1.6.0.1 May 28th 2019
https://uk.mathworks.com/matlabcentral/fileexchange/43097-newtonraphson
bisection.m
Bisection Method Root Finding, by Sky Sartorius. v1.16 Oct 23rd 2019
https://uk.mathworks.com/matlabcentral/fileexchange/28150-bisection-method-root-finding
[CUEN]
Investigating the Broadband Microwave Absorption of Nanodiamond Impurities
Jerome A. Cuenca, Evan Thomas, Soumen Mandal, Oliver Williams, and Adrian Porch
IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 63, NO. 12, DEC. 2015