How do I use muller method for solving multivariable equations?

I have two equations of 2 variables. I tried using 'solve' but it keeps on calculating for hours together with no results. I would like to use Muller method as I have used it before and I can define start points and number of iterations. I can also check the residual value. Can anyone please suggest, how can I use Muller for solving multivariable equations?

[Alan Stevens]

Easier to do this if we know what your equations are.

[Shreya Menon]

Sorry... The equations are:

1.) (epir*(diff(besselj(1,V))/(V*k0a*besselj(1,V)))-(diff(besselk(1,W))/(W*k0a*besselk(1,W))))*(mewr*(diff(besselj(1,V))/(V*k0a*besselj(1,V)))-(diff(besselk(1,W))/(W*k0a*besselk(1,W))))=((V^2+W^2)*(V^2+mewr*epir*W^2))/(V^4*W^4*k0a^4)

2.) ((2*(b+L*tand(alpha))/lambda0)^2)*(pi^2)*(mewr*epir-1)==(V^2+W^2)*k0a^2

V and W is to be found. These are equations for propagation characteristics of solid dielectric rod antenna. Kindly help in this regard.

Answer 1

[Alan Stevens]

I guess there are a few options.

1. If you have the Opimisation toolbox, use fsolve.
2. In your second equation replace V^2 + W^2 by, say, Rsq and solve for Rsq. Then express V as a function of W, knowing Rsq. Then use fzero to find W. The code structure might look something like the following (I'm unable to test it because I don't have your constants).

Rsq = ((2*(b+L*tand(alpha))/lambda0)^2)*(pi^2)*(mewr*epir-1)/k0a^2;
Vfn = @(W) sqrt(Rsq - W.^2);
W0 = ....;  % Insert your initial guess
W = fzero(@Wfn, W0);
V = Vfn(W);
    
function WW = Wfn(W)

V = Vfn(W);

WW = (epir*(diff(besselj(1,V))/(V*k0a*besselj(1,V))) ...
    -(diff(besselk(1,W))/(W*k0a*besselk(1,W))))*(mewr*(diff(besselj(1,V))/(V*k0a*besselj(1,V))) ...
    -(diff(besselk(1,W))/(W*k0a*besselk(1,W))))-((V^2+W^2)*(V^2+mewr*epir*W^2))/(V^4*W^4*k0a^4);

end

3. An alternative to using fzero with option 2 is to program the Muller method yourself. However, I suspect fzero is the better option.

[Shreya Menon]

Thank you for the reply and suggestions. This is the value of all constants:

epi0=8.854e-12;

mew0=4*pi*(10^(-7));

z0=sqrt(mew0/epi0);

c0=3.8e08;

alpha=15;

b=15.5677e-03;

L=40e-03;

epir=3.4;

epi1=epir*epi0;

n1=sqrt(epi1/epi0);

mewr=1;

mew1=mewr*mew0;

cntr=1;

for fcntr= 0.1e09:0.10e09:12.5e09

k0a=((2*pi*fcntr)/(3e08));

k1a=((epi1/epi0)^(1/2))*(2*pi*fcntr/c0);

lambda0=c0/fcntr;

omega=2*pi*fcntr;

rsq=((2*(b+L*tand(alpha))/lambda0)^2)*(pi^2)*(mewr*epir-1);

Vfn = @(W) sqrt(Rsq - W.^2 );

W0 = 1000;

W = fzero(@Wfn, W0 );

V = Vfn(W);

beta(cntr)=sqrt(omega^2*mew1*epi1/k0a^2-(Vb(cntr)/(b+L*tand(alpha)))^2);

cntr=cntr+1;

end

function WW = Wfn(W)

fcntr=1e09;

epi0=8.854e-12;

mew0=4*pi*(10^(-7));

z0=sqrt(mew0/epi0);

c0=3.8e08;

alpha=15;

b=15.5677e-03;

L=40e-03;

epir=3.4;

epi1=epir*epi0;

n1=sqrt(epi1/epi0);

mewr=1;

mew1=mewr*mew0;

lambda0=c0/fcntr;

omega=2*pi*fcntr;

k0a=((2*pi*fcntr)/(3e08));

Rsq=((2*(b+L*tand(alpha))/lambda0)^2)*(pi^2)*(mewr*epir-1);

Vfn = @(W) sqrt(Rsq - W.^2 );

V = Vfn(W);

WW = (epir*(diff(besselj(1,V))/(V*k0a*besselj(1,V ))) ...

-(diff(besselk(1,W))/(W*k0a*besselk(1,W))))*(mewr*(diff(besselj(1,V))/(V*k0a*besselj(1,V ))) ...

-(diff(besselk(1,W))/(W*k0a*besselk(1,W))))-((V^2+W^2)*(V^2+mewr*epir*W^2))/(V^4*W^4*k0a^4);

end

I tried what you suggested for 1GHz but getting this error:

"Operands to the || and && operators must be convertible to logical scalar values.

Error in fzero (line 308)

elseif ~isfinite(fx) || ~isreal(fx)

Error in propa_HE11_Chen (line 37)

W = fzero(@Wfn, W0 );"

I can't understand what the error is about.

[Alan Stevens]

You hadn't quite got Rsq right. It should be

Rsq=((2*(b+L*tand(alpha))/lambda0)^2)*(pi^2)*(mewr*epir-1)/k0a^2;

However, this doesn't help a lot! A reason for the error message could be that with the values for the constants you use it turns out that Rsq is a lot smaller than W^2 (about 0.001 compared with 1000^2). This results in a complex value for sqrt(Rsq - W^2).

[Shreya Menon]

Sir, V and W can be complex according to theory. So is there any other solver I can use to generate accurate answers? I am trying to match the theoretical results to 3D FEM simulation tool (HFSS) results. I know the value of beta should be around 1.84+0i where 1.84 is the normalized phase constant and 0i represents normalized attenuation constant.

[Alan Stevens]

Ah, fzero only deals with real numbers I'm afraid. I guess you need to look at the Optimisation toolbox.

[Shreya Menon]

Thank you sir for all the help you provided.