transcedental equation solve by muller method

3 次查看（过去 30 天）
shiv gaur
shiv gaur 2021-1-29
回答： shiv gaur 2021-12-12
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

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

回答（2 个）

Walter Roberson
Walter Roberson 2021-1-29
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
shiv gaur 2021-12-12
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

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by