solution of equation code of transcedental equation

1 Ansicht (letzte 30 Tage)
shiv gaur
shiv gaur am 12 Dez. 2021
Kommentiert: shiv gaur am 3 Jan. 2022
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
  6 Kommentare
shiv gaur
shiv gaur am 12 Dez. 2021
plot between t2 vs real(p)
shiv gaur
shiv gaur am 12 Dez. 2021
plot is the problem and loop and store the value if any one have idea plot the graph

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Torsten
Torsten am 12 Dez. 2021
Bearbeitet: Torsten am 12 Dez. 2021
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).
  5 Kommentare
shiv gaur
shiv gaur am 3 Jan. 2022
pl plot any one
shiv gaur
shiv gaur am 3 Jan. 2022
graph is like that

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by