solution of equation code of transcedental equation

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

[shiv gaur]

this is the solution of the equation my question is that plot by varying the value of t2 using for loop from 1e-3 micron to 1 micron this this equation is the the solution giving complex root we have to plot t2 vs real(p) and t2 vs imag (p)

pl plot by varying t2

[shiv gaur]

at value t2=1e-9 the roots of equation is 1.512 - 2.999D-16i

[shiv gaur]

this solution gives all the transcedental equation containing complex roots also

[Torsten]

And what is your problem ?

You don't know how to build a loop from t2=1e-9 to t2=1e-6 and store the corresponding results for p ?

[shiv gaur]

plot between t2 vs real(p)

[shiv gaur]

plot is the problem and loop and store the value if any one have idea plot the graph

Answer 1

[Torsten]

function kps3
T2 = 1e-9:1e-9:1e-6;
for j=1:numel(T2)
    
    t2 = T2(j);
    p0 = 0.5;
    p1 = 1;
    p2 = 1.5;
    TOL = 10^-8;
    N0 = 100; format long
    h1 = p1 - p0;
    h2 = p2 - p1;
    DELTA1 = (f(p1,t2) - f(p0,t2))/h1;
    DELTA2 = (f(p2,t2) - f(p1,t2))/h2;
    d = (DELTA2 - DELTA1)/(h2 + h1);
    i=3;
    while i <= N0 
        
        b = DELTA2 + h2*d;
        D = (b^2 - 4*f(p2,t2)*d)^(1/2);
        
        if abs(b-D) < abs(b+D)
            E = b + D;
        else
            E = b - D;
        end
        
        h = -2*f(p2,t2)/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,t2) - f(p0,t2))/h1;
        DELTA2 = (f(p2,t2) - f(p1,t2))/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
    P(j)=real(p);
end
plot(T2,P)
end
function  y=f(x,t2)

%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

Maybe you can use p_old = P(j-1) as starting point for the solution of your equation with t2_new = T2(j).

[shiv gaur]

sir you are requested to match the corrected fig

[shiv gaur]

function kps5
K0 = 1550e-6:1e-6:1600e-6;
for j=1:numel(K0)
    
    k0 = K0(j);
    p0 = 0.5;
    p1 = 1;
    p2 = 1.5;
    TOL = 10^-8;
    N0 = 100; format long
    h1 = p1 - p0;
    h2 = p2 - p1;
    DELTA1 = (f(p1,k0) - f(p0,k0))/h1;
    DELTA2 = (f(p2,k0) - f(p1,k0))/h2;
    d = (DELTA2 - DELTA1)/(h2 + h1);
    i=3;
    while i <= N0 
        
        b = DELTA2 + h2*d;
        D = (b^2 - 4*f(p2,k0)*d)^(1/2);
        
        if abs(b-D) < abs(b+D)
            E = b + D;
        else
            E = b - D;
        end
        
        h = -2*f(p2,k0)/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,k0) - f(p0,k0))/h1;
        DELTA2 = (f(p2,k0) - f(p1,k0))/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
    P(j)=real(p);
    R(j)=real(r)
end
plot(K0,R)
end
function  y=f(x,k0)
n0=1.3707;
n1=1.3;
n2=1.59;
n3=1.45
n4=3.46;
na=1.36;
t1=2.10e-6;
t2=0.198e-6;
t3=1e-6;
t4=0.002e-6;
k1=k0*sqrt(n1.^2-x.^2);
k2=k0*sqrt(n2.^2-x.^2);
k3=k0*sqrt(n3.^2-x.^2);
k4=k0*sqrt(n4.^2-x.^2);
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)*%i +(k2)*cos(t1*k1)*sin(t2*k2)*%i;
m22=cos(t1*k1)*cos(t2*k2)-(k1/k2)*sin(t1*k1)*sin(t2*k2);
m34= cos(t3*k3)*cos(t4*k4)-(k4/k3)*sin(t3*k3)*sin(t4*k4);
m32=(1/k4)*(cos(t3*k3)*sin(t4*k4)*1i) +(1/k3)*(cos(t4*k4)*sin(t3*k3)*1i);
m23= (k3)*cos(t4*k4)*sin(t3*k3)*1i +(k4)*cos(t3*k3)*sin(t4*k4)*1i;
m33=cos(t3*k3)*cos(t4*k4)-(k3/k4)*sin(t3*k3)*sin(t4*k4);
M11=m11*m34+m12*m23;
M12=m11*m32+m12*m33;
M21=m21*m34+m22*m23;
M22=m21*m32+m22*m33;
g0=sqrt(x.^2-n0.^2)*k0;
ga= sqrt(x.^2-na.^2)*k0;
y= 1i*(g0*M11/n0^2+ga*M22/na^2)-M21+(g0*ga/n0^2*na^2)*M12 ;
r=(na^2*g0*M11-n0^2*ga*M22+g0*ga*M12-na^2*n0^2*M21)/(na^2*g0*M11+n0^2*ga*M22+g0*ga*M12+na^2*n0^2*M21);
end

[shiv gaur]

same type of problem pl help to plot between k0 vs r that is given in equation

[shiv gaur]

pl plot any one

[shiv gaur]

graph is like that