Solving transcendental equation numerically

10 visualizzazioni (ultimi 30 giorni)
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 Commenti
Jan
Jan il 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 il 29 Mar 2022
Modificato: Torsten il 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.

Accedi per commentare.

Risposta accettata

Jan
Jan il 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 .*

Più risposte (1)

Torsten
Torsten il 29 Mar 2022
Modificato: Torsten il 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

Categorie

Scopri di più su Mathematics in Help Center e File Exchange

Prodotti


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by