solution of equation code of transcedental equation

1 次查看（过去 30 天）
shiv gaur
shiv gaur 2021-12-12
评论： shiv gaur 2022-1-3
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 个评论
显示 4更早的评论
shiv gaur
shiv gaur 2021-12-12
plot is the problem and loop and store the value if any one have idea plot the graph

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2021-12-12
编辑：Torsten 2021-12-12
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 个评论
显示 3更早的评论

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by