单模光纤的特征方程以及MATLAB求解

在之前的文章中，我们写出了单模光纤的特征方程以及对应的导波模式

这里我们在MATLAB中求解特征方程并表示出几个归一化参数的变化曲线

用到的公式是单模光纤导波模式的特征方程：

clear
close all
tic
Vmax = 10;
N =100;

for j =1:N
    V(j) = j/N*Vmax;
    Vtemp = V(j);

    Utemp = NaN;
    i = 0;

    while (isnan(Utemp) && i<N+1)
        init =Vtemp*(i+1)/N-1e-16;
        try
            Utemp = fzero(@(Utemp) ...
                besselj(0,Utemp)/(Utemp*besselj(1,Utemp)) - ...
                besselk(0,sqrt(Vtemp^2-Utemp^2))/ ...
                (sqrt(Vtemp^2-Utemp^2)*besselk(1,sqrt(Vtemp^2-Utemp^2))),init);
        catch
        end
        i = i+1;
    end
    U(j) = abs(Utemp);
end

W = sqrt(V.^2-U.^2);
Ymax = ceil(max([U,W]));

figure
subplot(1,2,1)
plot(V,U,'r');
axis equal
axis([0 Vmax 0 Ymax])
xlabel('v')
ylabel('u')
title('Lp_{01}  v-u')

subplot(1,2,2)
plot(V,W,'r');
axis equal
axis([0 Vmax 0 Ymax])
xlabel('v')
ylabel('w')
title('Lp_{01}  v-w')

toc

fzero函数能求解fun(x) = 0 的点 x。

使用方式：

x = fzero(fun,x0)

其中x0表示初始值

我们从0-10之间选取了N（100）个点，对于每一个点都有一个Vtemp与之对应

对于每一个确定的Vtemp，特征方程都有一个Utemp与之对应

结果如下：

 try catch语句在java中常见于抛出错误

在预见到会产生错误的前提下，把可能引发异常的语句封装在 try 语句块中，用以捕获可能发生的异常。

catch 语句会依据所拋出异常对象的类型进行捕获，并处理。

处理之后，程序会跳过 try 语句块中剩余的语句，转到 catch 语句块后面的第一条语句开始执行。

如果 try 语句块中没有异常发生，那么 try 块正常结束，后面的 catch 语句块被跳过，程序将从 catch 语句块后的第一条语句开始执行。

ceil函数是取整函数，将输入参数的每个元素四舍五入到大于或等于该元素的最接近整数。这里用来规划坐标轴的范围

代码的输出结果就是单模光纤UVW三个元素之间的关系，我们可以在命令行窗口看到具体的数值结果：

[V' U' W']

截取前一段来观察：

在V很小时，U和W没有数值解，也就是不存在合适的U和W形成导波模式

在曲线的起始端，对应的数值分别为Vc和Uc

随着V不断增加，U和W都增大，但是，W的增加近似于线性（在V较大的情况下），U的增加趋近于一个有限值

我们可以用多项式拟合得到近似的直线，并和通过特征方程求解的结果作对比：

NN = 15:24
x = V(NN);
y = W(NN);
p = polyfit(x,y,1);
f = polyval(p,x);
maxerr = max(y-f)
figure
plot(x,y,'o',x,f,'-')
xlabel('v')
ylabel('w')

polyfit函数用于拟合函数，返回次数为 n 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（基于最小二乘指标）。p 中的系数按降幂排列，p 的长度为 n+1

这里我们设置阶数为1，也就是直线拟合

结果如下：

polyval函数计算多项式 p 在 x 的每个点处的值。参数 p 是长度为 n+1 的向量，其元素是 n 次多项式的系数（降幂排序）

通过这两个函数组合，我们得到了拟合直线的表达式以及拟合直线的在对应点的结果

在命令行窗口读出最大的误差，可以看到最大的误差是数量级。

参考《高等光学仿真——光波导、激光》