请问FindRoot外面套一个For循环的问题
lamda = 1.55 10^-6;k0 = 2*Pi/lamda;
n1 = 1.4677;(*纤芯折射率*)
n2 = 1.4628;(*包层折射率*)
n3 = 0.469 + 9.32*I;(*银折射率*)
a1 = 4.1 10^-6;(*纤芯半径*)
a2 = 62.5 10^-6;(*包层半径*)
d = 40 10^-9;(*金属厚度*)
a3 = a2 + d;
mu = Pi*4 10^-7;(*真空磁导率*)
epsi0 = 8.85 10^-12;(*介电常数*)
n4 = 1.330;
neffcl = neffclre + neffclim*I;
betacl = k0*neffcl;
omega = 2*Pi*299792458/lamda;
epsi1 = n1^2*epsi0;
epsi2 = n2^2*epsi0;
epsi3 = n3^2*epsi0;
epsi4 = n4^2*epsi0;
u1 = k0*Sqrt;
u2 = k0*Sqrt;
u3 = k0*Sqrt;
w4 = k0*Sqrt;
Iua111 = BesselI;
Iua121 = BesselI;
Iua122 = BesselI;
Iua132 = BesselI;
Iua133 = BesselI;
IIua111 = (BesselI + BesselI )/2;
IIua121 = (BesselI + BesselI )/2;
IIua122 = (BesselI + BesselI)/2;
IIua132 = (BesselI + BesselI)/2;
IIua133 = (BesselI + BesselI)/2;
Kua121 = BesselK ;
Kua122 = BesselK ;
Kua132 = BesselK ;
Kua133 = BesselK ;
Kwa143 = BesselK ;
KKua121 = -(BesselK + BesselK )/2;
KKua122 = -(BesselK + BesselK )/2;
KKua132 = -(BesselK + BesselK )/2;
KKua133 = -(BesselK + BesselK )/2;
KKwa143 = -(BesselK + BesselK )/2;
H1 = (betacl*Kwa143*
Kua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3)*(u3/u2*IIua132*
Kua122 - u3^2/u2^2*Iua132*KKua122) - (betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3)*(u3/u2*KKua132*
Kua122 - u3^2/u2^2*Kua132*KKua122) + (betacl*Iua132*
Kua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(w4/u3*KKwa143*
Kua133 - w4^2*epsi3/u3^2/epsi4*Kwa143*KKua133) - (betacl*
Kua132*Kua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(w4/u3*KKwa143*
Iua133 - w4^2*epsi3/u3^2/epsi4*Kwa143*IIua133);
H2 = (betacl*Kwa143*
Kua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3)*(u3/u2*IIua132*
Iua122 - u3^2/u2^2*Iua132*IIua122) - (betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3)*(u3/u2*KKua132*
Iua122 - u3^2/u2^2*Kua132*IIua122) + (betacl*Iua132*
Iua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(w4/u3*KKwa143*
Kua133 - w4^2*epsi3/u3^2/epsi4*Kwa143*KKua133) - (betacl*
Kua132*Iua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(w4/u3*KKwa143*
Iua133 - w4^2*epsi3/u3^2/epsi4*Kwa143*IIua133);
H3 = (betacl*Iua132*Kua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(betacl*
Kwa143*Kua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3) - (betacl*
Kua132*Kua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3) + (u3/u2*IIua132*
Kua122 -
u3^2*epsi2/u2^2/epsi3*Iua132*KKua122)*(w4/u3*KKwa143*Kua133 -
w4^2*epsi3/u3^2/epsi4*Kwa143*KKua133) - (u3/u2*KKua132*Kua122 -
u3^2*epsi2/u2^2/epsi3*Kua132*KKua122)*(w4/u3*KKwa143*Iua133 -
w4^2*epsi3/u3^2/epsi4*Kwa143*IIua133);
H4 = (betacl*Iua132*Iua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(betacl*
Kwa143*Kua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3) - (betacl*
Kua132*Iua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/epsi4/u3/a3) + (u3/u2*IIua132*
Iua122 -
u3^2*epsi2/u2^2/epsi3*Iua132*IIua122)*(w4/u3*KKwa143*Kua133 -
w4^2*epsi3/u3^2/epsi4*Kwa143*KKua133) - (u3/u2*KKua132*Iua122 -
u3^2*epsi2/u2^2/epsi3*Kua132*IIua122)*(w4/u3*KKwa143*Iua133 -
w4^2*epsi3/u3^2/epsi4*Kwa143*IIua133);
M1 = (betacl*Iua132*Kua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(betacl*
Kwa143*Kua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3) - (betacl*Kua132*
Kua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3) + (u3/u2*IIua132*Kua122 -
u3^2/u2^2*Iua132*KKua122)*(w4/u3*KKwa143*Kua133 -
w4^2/u3^2*Kwa143*KKua133) - (u3/u2*KKua132*Kua122 -
u3^2/u2^2*Kua132*KKua122)*(w4/u3*KKwa143*Iua133 -
w4^2/u3^2*Kwa143*IIua133);
M2 = (betacl*Iua132*Iua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(betacl*
Kwa143*Kua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3) - (betacl*Kua132*
Iua122*(u3^2/u2^2 - 1)/omega/epsi3/u2/a2)*(betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3) + (u3/u2*IIua132*Iua122 -
u3^2/u2^2*Iua132*IIua122)*(w4/u3*KKwa143*Kua133 -
w4^2/u3^2*Kwa143*KKua133) - (u3/u2*KKua132*Iua122 -
u3^2/u2^2*Kua132*IIua122)*(w4/u3*KKwa143*Iua133 -
w4^2/u3^2*Kwa143*IIua133);
M3 = (betacl*Kwa143*
Kua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3)*(u3/u2*IIua132*Kua122 -
u3^2*epsi2/u2^2/epsi3*Iua132*KKua122) - (betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3)*(u3/u2*KKua132*Kua122 -
u3^2*epsi2/u2^2/epsi3*Kua132*KKua122) + (betacl*Iua132*
Kua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(w4/u3*KKwa143*Kua133 -
w4^2/u3^2*Kwa143*KKua133) - (betacl*Kua132*
Kua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(w4/u3*KKwa143*Iua133 -
w4^2/u3^2*Kwa143*IIua133);
M4 = (betacl*Kwa143*
Kua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3)*(u3/u2*IIua132*Iua122 -
u3^2*epsi2/u2^2/epsi3*Iua132*IIua122) - (betacl*Kwa143*
Iua133*(w4^2/u3^2 - 1)/omega/mu/u3/a3)*(u3/u2*KKua132*Iua122 -
u3^2*epsi2/u2^2/epsi3*Kua132*IIua122) + (betacl*Iua132*
Kua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(w4/u3*KKwa143*Kua133 -
w4^2/u3^2*Kwa143*KKua133) - (betacl*Kua132*
Iua122*(u3^2/u2^2 - 1)/omega/mu/u2/a2)*(w4/u3*KKwa143*Iua133 -
w4^2/u3^2*Kwa143*IIua133);
R1 = u2^2/u1^2*Iua121*IIua111 - u2/u1*IIua121*Iua111;
T1 = u2^2/u1^2*Kua121*IIua111 - u2/u1*KKua121*Iua111;
U1 = betacl*Iua121*Iua111*(u2^2/u1^2 - 1)/omega/epsi2/u1/a1;
V1 = betacl*Kua121*Iua111*(u2^2/u1^2 - 1)/omega/epsi2/u1/a1;
R2 = u2^2/u1^2*epsi1/epsi2*Iua121*IIua111 - u2/u1*IIua121*Iua111;
T2 = u2^2/u1^2*epsi1/epsi2*Kua121*IIua111 - u2/u1*KKua121*Iua111;
U2 = betacl*Iua121*Iua111*(u2^2/u1^2 - 1)/omega/mu/u1/a1;
V2 = betacl*Kua121*Iua111*(u2^2/u1^2 - 1)/omega/mu/u1/a1;
xicl1 = (-R1*H1 + T1*H2 + U1*H3 - V1*H4)/(R1*M1 - T1*M2 - U1*M3 +
V1*M4);
xicl2 = (-R2*H3 + T2*H4 + U2*H1 - V2*H2)/(R2*M3 - T2*M4 - U2*M1 +
V2*M2);
x = xicl1 - xicl2;
x1 = Re;
x2 = Im;
FindRoot[{x1,x2},{{neffclre,1.333},{neffclim,0.00001}}];
]
代码如上，结果是{neffclre -> 1.33017, neffclim -> 0.0000172055}
但我把FindRoot[{x1,x2},{{neffclre,1.333},{neffclim,0.00001}}];
换成
For[i = 1, i < 133, i++, neffclbase = 1.330 + 0.001*i;
FindRoot[{x1, x2}, {{neffclre, neffclbase}, {neffclim, 0.00001}}];
]
就会出现
FindRoot::lstol: 线搜索把步长降低到由 AccuracyGoal 和 PrecisionGoal 指定的容差范围内，但是无法找到 merit 函数的充足的降低. 您可能需要多于 MachinePrecision 位工作精度以满足这些容差.
请问是怎么回事?
页:
[1]