Effacer les filtres
Effacer les filtres

Solving transcendental equation numerically

26 vues (au cours des 30 derniers jours)
Jose Iglesias
Jose Iglesias le 29 Mar 2022
Modifié(e) : Jose Iglesias le 30 Mar 2022
I am trying to plot a transcendental equation numerically in order to plot the propagation constants B+ and B-. I revamped my code after modifying thd equation where I can enter in the H0 values, but when I enter in the H0 values I get an error in line 19. My code and error are found below. Can anyone help me alleviate my mistake which is a matrix dimension issue? Thank you in advance. The H0 values you would input are [0 200 400 600 800 1000 1200 1400 1600]. The plot should look like the one below.
Here is the transcendental equation being used to create the Matlab code:
Here is the Matlab code with the given parameters to plot the forward and reverse propagation constants B+ and B-:
clear all;
H0=input('Enter value of H0:')
f0=2.8e6.*H0; %larmor precession frequency
fm=2.8e6*1700; %saturation magnetization frequency
f=10e9; %operating frequency
w=2*pi*f; %2pif
c=0;
a=1e-2;
t=a/2;
d=t;
k0=209.4395; %1/m
b=0*k0:0.000025*k0:4*k0;
e0=8.854e-12; %F/m permittivity of free space
u0=4*pi*1e-7; %H/m permeability of free space
e=e0*u0;
u=u0.*(1+ ((f0.*fm)./((f0.^2)-f^2))); %permeability tensor element
k=u0.*((f*fm)./((f0.^2)-f^2)); %permeability tensor element
ue=(u.^2-k.^2)./(u); %effective permeability
kf=sqrt(w.^2*ue.*e-b.^2); %kf is cutoff wavenumber for ferrite
ka=sqrt(k0.^2-b.^2); %ka is cutoff wavenumber for air
yplus=((u.*ue./u0).*(ka./k).*cot(ka.*d))+((u.*kf./k).*cot(kf.*t))+b;
yminus=((u.*ue./u0).*(ka./-k).*cot(ka.*d))+((u.*kf./-k).*cot(kf.*t))+b;
yp=real(yplus);
ym=real(yminus);
plot(b/k0,yp,'-',b/k0,ym,'--');
grid on;
axis([0 4 -1 1]);
Matrix dimensions must agree.
Error in ME_8p9 (line 19)
kf=sqrt(w.^2*ue.*e-b.^2); %kf is cutoff wavenumber for ferrite
  3 commentaires
Jan
Jan le 29 Mar 2022
@Jose Iglesias: Whenever you mention an error in the forum, post an exact copy of the message. Paraphrasing the message is less useful.
Torsten
Torsten le 29 Mar 2022
Modifié(e) : Torsten le 29 Mar 2022
Define H0 first as a symbolic variable and substitute values for H0 after you got a solution for your system of equations depending on B and H0.
The variable d is undefined.

Connectez-vous pour commenter.

Réponse acceptée

Jan
Jan le 29 Mar 2022
You can check this using the debugger. Type this in the command window:
dbstop if error
Now run the code again. When Matlab stops at the error, evaluate the parts of the currently processed line:
eqn3 = ((ka.*cot(ka.*c))*((kf*cot(kf.*t))/(u0.*ue))+((k.*B)/(u0.*u.*ue)));
size((ka .* cot(ka .* c)))
size((kf * cot(kf.*t)) / (u0 .* ue))
size((k.*B) / (u0.*u.*ue))
Maybe you want to replace the * by a the elementwise .*

Plus de réponses (1)

Torsten
Torsten le 29 Mar 2022
Modifié(e) : Torsten le 30 Mar 2022
Your code returns NaN or +/- Inf when trying to evaluate equation (9.79) for the constants you gave.
Before you do not solve this problem, it does not make sense to continue.
Ho = 750;
B=10:0.1:100;
for i=1:numel(B)
f(i) = fun(B(i),Ho);
end
plot(B,f)
function res = fun(B,Ho)
c=0;
d=5e-3;
a=1e-2;
t=a/2;
%Ho=0:1500;
fo=2.8e6.*Ho;%larmor precession frequency
fm=2.8e6*1700;%saturation magnetization frequency
u0=4*pi*1e-7; %H/m permeability of free space
f=10e9;%operating frequency
u=u0*(1+ ((fo.*fm)./((fo.^2)-f^2))); %permeability tensor element
k=u0*((f*fm)./((fo.^2)-f^2)); %permeability tensor element
ue=(u.^2-k.^2)./(u);%effective permeability
ko=209.4395;%1/m
ka=sqrt(ko.^2-B.^2);%ka is cutoff wavenumber for air
w=6.283e10;%2pif
eo=8.85e-12; %F/m permittivity of free space
er=13;
e=eo*er;
kf=sqrt(w.^2*ue*e-B^2); %kf is cutoff wavenumber for ferrite
eqn1=(kf./ue).^2;
eqn2=((k.*B)./(u.*ue)).^2;
eqn3=(ka*cot(ka*c).*((kf./(u0.*ue)).*cot(kf.*t)+((k.*B)./(u0.*u.*ue))));
eqn4=(ka/u0).^2;
eqn5=((cot(ka.*c).*cot(ka.*d)));
eqn6=(ka*cot(ka*d).*((kf./(u0.*ue)).*cot(kf.*t)-((k.*B)./(u0.*u.*ue))));
eqn7=(eqn1+eqn2-eqn3-eqn4.*eqn5)-eqn6;
res = eqn7;
end

Catégories

En savoir plus sur Mathematics dans Help Center et File Exchange

Produits


Version

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by