program is not showing graph

1 次查看（过去 30 天）
shiv gaur
shiv gaur 2022-1-11
回答： Walter Roberson 2022-1-11
function metal3
theta1=linspace(10,70,20);
%T3=1e-9:1e-9:1e-6;
for j=1:numel(theta1)
theta=theta1(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,theta) - f(p0,theta))/h1;
DELTA2 = (f(p2,theta) - f(p1,theta))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;
while i <= N0
b = DELTA2 + h2*d;
D = (b^2 - 4*f(p2,theta)*d)^(1/2);
if abs(b-D) < abs(b+D)
E = b + D;
else
E = b - D;
end
h = -2*f(p2,theta)/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,theta) - f(p0,theta))/h1;
DELTA2 = (f(p2,theta) - f(p1,theta))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i = i + 1;
end
if i > N0
formatSpec = string('The method failed after N0 iterations,N0= %d \n');
end
P(j)=real(p);
end
plot(theta ,r)
end
function y=f(x,theta)
k0=(2*pi/0.6328)*1e6;
n0=1.3707;
n1=1.30;
n2=1.59;
n3=0.065+4*1i;
n4=3.48;
na=1.36;
t1=2.10e-6;
t2=0.198e-6;
t3=1e-6;
t4=0.002e-6;
x=n0*k0*sin(theta);
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)*1i +(k2)*cos(t1*k1)*sin(t2*k2)*1i;
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=(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);
r=y.^2;
end

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

回答（1 个）

Walter Roberson
Walter Roberson 2022-1-11
You have
function y=f(x,theta)
You are passing in a variable named x
k0=(2*pi/0.6328)*1e6;
That's a value in the millions
n0=1.3707;
n1=1.30;
That's a value near 1
n2=1.59;
n3=0.065+4*1i;
n4=3.48;
na=1.36;
t1=2.10e-6;
t2=0.198e-6;
t3=1e-6;
t4=0.002e-6;
x=n0*k0*sin(theta);
You have not used the x that you passed in. Instead you have overwritten x with a value that is determined in part by k0, which is in the millions, so x is going to end up in the millions.
k1=k0*sqrt(n1.^2-x.^2);
Remember n1 is near 1, and you are subtracting a million squared, and taking the square root. So you are going to end up with a value that is roughly complex 1 million.
m11= cos(t1*k1)*cos(t2*k2)-(k2/k1)*sin(t1*k1)*sin(t2*k2);
cos of complex 1 million is infinite . Remember that cos of an imaginary number has to do with exponentials .

类别

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

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by