transcedental equation solve by muller method

I have written program but having some problems for result and I have to plot the variation of t2 with 0.1e-6 to 1e-6 with real value of roots

function y=f(x)
 lamda=2*pi/0.6328e-6   
k0=lamda;

n1=1.521;

n2=4.1-1i*0.211;
ns=1.512;

nc=1;
k1=sqrt(n1.^2*k0.^2-x.^2);
k2=sqrt(n2.^2*k0.^2-x.^2);
t1=1.5e-6;

m11= cos(t1*k1)*cos(t2*k2)+(k2/k1)*sin(t1*k1)*sin(t2*k2);
m12= (cos(t2*k2)*sin(t1*k1)*1i)/k1 + (cos(t1*k1)*sin(t2*k2)*1i)/k2;
m21= k1*cos(t2*k2)*sin(t1*k1)*1i + k2*cos(t1*k1)*sin(t2*k2)*1i;
m22= cos(t1*k1)*cos(t2*k2) - (k1/k2)*sin(t1.*k1)*sin(t2*k2);
gs=x.^2-ns.^2*k0.^2;
gc= x.^2-nc.^2*k0.^2;
 
   y=  @(x) 1i*(gs*m11+gc*m22)-m21+gc*gs*m12 
 

    
    
    
end
  t2=0.081

p0 = 0;
p1 = 1;
p2 = 2;
TOL = 10^-5;
N0 = 100;


h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1) - f(p0))/h1;
DELTA2 = (f(p2) - f(p1))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;

while i <= N0

      b = DELTA2 + h2*d;
      D = (b^2 - 4*f(p2)*d)^(1/2);

      if abs(b-D) < abs(b+D)
          E = b + D;
      else
          E = b - D;
      end

      h = -2*f(p2)/E;
      p = p2 + h;

      if abs(h) < TOL
          p
         
         disp(p) 
          break
      end

      p0 = p1;
      p1 = p2;
      p2 = p;
      h1 = p1 - p0;
     h2 = p2 - p1;
     DELTA1 = (f(p1) - f(p0))/h1;
      DELTA2 = (f(p2) - f(p1))/h2;
     d = (DELTA2 - DELTA1)/(h2 + h1);
      i = i + 1;
     
   
end

Answer 1

[Walter Roberson]

Your k2 is complex by design.

The product of t2 and k2 turns out to be in the tens of thousands for your starting x value, 1.

cos() and sin() of complex numbers in the tens of thousands gives you infinities and nan values.

So you cannot solve this in double precision. You would need to move to symbolic computing to have a chance.

Also, if you need to initialize t2 outside of f, then you need to pass it in to f.

And you should not have the @() on the assignment to y inside f.

[shiv gaur]

real part of this is to be 1.55 approx

Answer 2

[shiv gaur]

function kps3

p0 = 0.5;

p1 = 1;

p2 = 1.5;

TOL = 10^-8;

N0 = 100; format long

h1 = p1 - p0;

h2 = p2 - p1;

DELTA1 = (f(p1) - f(p0))/h1;

DELTA2 = (f(p2) - f(p1))/h2;

d = (DELTA2 - DELTA1)/(h2 + h1);

i=3;

while i <= N0

b = DELTA2 + h2*d;

D = (b^2 - 4*f(p2)*d)^(1/2);

if abs(b-D) < abs(b+D)

E = b + D;

else

E = b - D;

end

h = -2*f(p2)/E;

p = p2 + h;

if abs(h) < TOL

disp(p)

break

end

p0 = p1;

p1 = p2;

p2 = p;

h1 = p1 - p0;

h2 = p2 - p1;

DELTA1 = (f(p1) - f(p0))/h1;

DELTA2 = (f(p2) - f(p1))/h2;

d = (DELTA2 - DELTA1)/(h2 + h1);

i=i+1

end

if i > N0

formatSpec = string('The method failed after N0 iterations,N0= %d \n');

// fprintf(formatSpec,N0);

end

function y=f(x)

t2=1e-9

k0=(2*pi/0.6328)*1e6;

n1=1.521;

n2=2.66;

ns=1.512;

nc=0.15-1i*3.2;

k1=k0*sqrt(n1.^2-x.^2);

k2=k0*sqrt(n2.^2-x.^2);

t1=1.5e-6;

m11= cos(t1*k1)*cos(t2*k2)-(k2/k1)*sin(t1*k1)*sin(t2*k2);

m12=(1/k2)*(cos(t1*k1)*sin(t2*k2)*1i) +(1/k1)*(cos(t2*k2)*sin(t1*k1)*1i);

m21= (k1)*cos(t2*k2)*sin(t1*k1)*1i +(k2)*cos(t1*k1)*sin(t2*k2)*1i;

m22=cos(t1*k1)*cos(t2*k2)-(k1/k2)*sin(t1*k1)*sin(t2*k2);

gs=(x.^2-ns.^2)*k0.^2;

gc= (x.^2-nc.^2)*k0.^2;

y= 1i*(gs*m11+gc*m22)-m21+gc*gs*m12 ;

end

end