基于MATLAB的光纤光传播特性仿真

%% 光纤几何参数
a = 9e-6;       % 核心半径 (单模光纤典型值)
lambda = 1.55e-6; % 波长 (1550 nm)
n_core = 1.46;  % 核心折射率
n_cladding = 1.45; % 包层折射率
​
%% 数值参数
L = 1e-3;       % 光纤长度 (1 mm)
dz = 1e-6;      % 空间步长 (1 μm)
dz_opt = 1e-3;  % 自适应步长 (1 mm)

% 阶跃型折射率分布
n = @(r) n_core^2 - (n_core^2 - n_cladding^2)*(r<=a);
​
% 渐变型折射率分布 (平方律分布)
alpha = 0.2;    % 折射率分布参数
n_grad = @(r) sqrt(n_core^2 - (n_core^2 - n_cladding^2)*exp(-alpha*r.^2));

function dUdz = NLSE(z, U, beta2, gamma, alpha)
    % 参数说明：
    % beta2: 群速度色散系数
    % gamma: 非线性系数
    % alpha: 吸收系数
    
    N = length(U);
    dUdz = zeros(N,1);
    
    % 分段处理
    for i = 1:N-1
        % 空间导数 (4阶有限差分)
        d2Udz2 = (U(i+2) - 4*U(i+1) + 6*U(i) - 4*U(i-1) + U(i-2))/(dz^2);
        
        % 非线性项
        nonlinear = -gamma*abs(U(i)).^2*U(i);
        
        % 色散项
        dispersion = 1i*beta2/2*d2Udz2;
        
        % 合成方程
        dUdz(i) = 1i*(dispersion + nonlinear) - 1i*alpha*U(i);
    end
end

%% 初始条件
U0 = sech(lambda*z).*exp(1i*2*pi*1e9*z); % 初始高斯脉冲
​
%% 参数设置
beta2 = -20e-27;  % 色散系数 (ps²/m)
gamma = 1.2e-3;   % 非线性系数 (W⁻¹m⁻¹)
alpha = 0.01;     % 吸收系数 (dB/m)
​
%% 自适应步长求解
z = 0:dz:1e-3;
U = zeros(size(z));
U(1) = U0;
​
for i = 1:length(z)-1
    % 自适应步长控制
    if i == 1
        dz_step = dz_opt;
    else
        dz_step = dz;
    end
    
    % 求解ODE
    [~, U(i+1)] = ode45(@(dz, U) NLSE(dz, U, beta2, gamma, alpha), dz_step, U(i));
end

%% 群速度色散仿真
beta2 = linspace(-30, 30, 100)*1e-27; % 色散参数范围
pulse_width = zeros(size(beta2));
​
for i = 1:length(beta2)
    [~, U] = simulate_pulse(beta2(i));
    pulse_width(i) = FWHM(U);
end
​
plot(beta2 * 1e27, pulse_width*1e-12);
xlabel('β_2 (ps²/m)'); ylabel('脉冲展宽 (ps)');
title('色散对脉冲展宽的影响');

%% 光孤子仿真
L = 1e-2; % 光纤长度 (10 mm)
z = 0:dz:1e-2;
​
% 初始条件 (高斯脉冲)
U0 = (1/sqrt(2))*sech(lambda*z).*exp(1i*2*pi*1e9*z);
​
% 参数设置
beta2 = -20e-27;  % 色散系数
gamma = 1.2e-3;   % 非线性系数
alpha = 0.01;     % 吸收系数
​
% 求解
[U] = simulate_pulse(beta2, gamma, alpha, L);
​
% 绘制时域波形
figure;
plot(z*1e3, abs(U).^2);
xlabel('传播距离 (mm)'); ylabel('光强 (a.u.)');
title('光孤子自相位调制效应');

%% LP模式场分布
l = 1; m = 1; % LP模式阶数
V = 2.405;    % 归一化频率
​
% Hermite-Gaussian模式
x = linspace(-5*a,5*a,200);
y = linspace(-5*a,5*a,200);
[X,Y] = meshgrid(x,y);
R = sqrt(X.^2 + Y.^2);
Phi = atan2(Y,X);
​
E = (2/(pi*a^2))^(1/2) * (R/a).^l .* exp(-R.^2/(2*a^2)) .* cos(l*Phi);
​
% 可视化
figure;
imagesc(x*1e6, y*1e6, abs(E).^2);
axis image;
title('LP_{11}模场分布');

%% 模式有效折射率计算
beta = zeros(1,100);
lambda_scan = linspace(1.3e-6, 1.6e-6, 100);
​
for i = 1:length(lambda_scan)
    beta(i) = sqrt(n_core^2*(lambda_scan(i)/lambda)^2 - (m*pi/a)^2);
end
​
plot(lambda_scan*1e9, beta*1e-9);
xlabel('波长 (nm)'); ylabel('传播常数 (rad/m)');
title('模式有效折射率色散曲线');

%% GPU并行计算
if gpuDeviceCount > 0
    U_gpu = gpuArray(U);
    beta2_gpu = gpuArray(beta2);
    gamma_gpu = gpuArray(gamma);
    
    [~, U_result] = ode45(@(dz, U) NLSE(dz, U, beta2_gpu, gamma_gpu, alpha), dz_opt, U_gpu);
    U = gather(U_result);
end

%% 三维传播仿真
[X,Y,Z] = ndgrid(linspace(-5*a,5*a,100), linspace(-5*a,5*a,100), linspace(0,L,50));
E3D = repmat(E, [1,1,length(Z)]);
​
for i = 1:length(Z)-1
    E3D(:,:,i+1) = propagate_mode(E3D(:,:,i), dz, beta2, gamma);
end
​
slice(X,Y,Z,abs(E3D).^2, [0.5,0.75](@ref));
shading interp;
title('三维光场传播仿真');