超越方程的数值求解

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

对于单一解的简单方程，我们可以使用fsolve直接求解，当我们碰到多支解的超越方程，直接使用fsolve，运行结果往往并不理想。我们分享一种求解超越方程的数值方法，并且展示MATLAB代码和运行结果。此外，欢迎交流效率更高的算法。

一、方法

求解超越方程，形如$$f(x)=0，$$（1）用数值方法初筛出一系列接近零点的解，按照大小排序分类，挑选出最接近零点的值；
（2）在此解附近取初筛精度的间距，加密细筛零点解；
（3）细筛出的零点解作为初始值，使用fsolve求解。

二、示例：代码和结果

这里我就不把示例中的超越方程公式敲出来了，所有必要的代码注释都在程序中，自行理解。

%输入参数
L_lambda_p=0.04;
options=optimoptions("fsolve","Algorithm","trust-region");

%定义要求解的方程 f=0
f=@(K_para,Omega) 1-(sqrt(1-1./(Omega.^2-K_para.^2))-(1-1./Omega.^2)).^2./(sqrt(1-1./(Omega.^2-K_para.^2))+(1-1./Omega.^2)).^2.*exp(2j*2*pi.*sqrt(Omega.^2-K_para.^2).*L_lambda_p);

%求解
Omega=0.0002:0.0005:1.2002;Y={};
for K_para=0:0.1:20
    %初筛选
    Dispersion=@(Omega)f(K_para,Omega);
    Omega_scatter=Dispersion(Omega);
    [g_scatter,index]=sort(abs(Omega_scatter));   %选取函数绝对值从小到大排序
    Omega_index=Omega(index(1:8));   %取函数绝对值最小的几个频率取值
    Omega_index=sort(Omega_index);   %对第一次筛选出来的频率值从小到大排序
    c=1;C={};count=1;   %数有多少相邻取值区域
    for i=2:length(Omega_index)
        if Omega_index(i)-Omega_index(i-1)>=0.00051   %第一次筛选出来的频率取值不相邻判据
            C=[C,Omega_index(c:i-1)];   %把相邻的频率值以数组形式放进Cell里
            c=i;
            count=count+1;
        end
    end
    C=[C,Omega_index(c:length(Omega_index))];

    % 细筛选
    Y_index2=[];   %第二次筛选值
    for i=1:count
        D=C{i};   %从Cell中从小到大取出连续数组
        Omega_scatter2=Dispersion(D);
        [g_scatter2,index2]=sort(abs(Omega_scatter2));
        Omega_index2=D(index2(1));
        Y_index2=[Y_index2,Omega_index2];
    end

    %fsolve精确求解
    Omega_K_para=fsolve(Dispersion,Y_index2,options);   %对第二次筛选值求解【数组】
    Y=[Y,Omega_K_para];
end

%处理求解结果
for j=1:length(Y)
    Y{j}=real(Y{j});
end
Y_TM_p=Y;
% save("TM_004.mat","Y_TM_p")

结果：

% plot
K_para=0:0.1:20;
for i=1:length(Y_TM_p)
    y=transpose(Y_TM_p{i});
    scatter(K_para(i),y,"red","filled","SizeData",14);
    hold on
end
xlabel(gca,"$ck_\parallel/\omega_p$","Interpreter","latex","FontName","Times New Roman","FontSize",14)
ylabel(gca,"$\omega/\omega_p$","Interpreter","latex","FontName","Times New Roman","FontSize",14)
set(gca,"FontName","Times New Roman","FontSize",14,"Box","on","TickLength",[0.025,0.025],"XMinorTick","on","YMinorTick","on")
axis([-0.5,20.5,-0.05,1.05])

总结

我们提供一种数值求解超越方程的方法。