一维传输矩阵的推导以及其在超材料里的应用（代码解读）

秋学好问

已于 2024-11-08 18:26:27 修改

阅读量546

点赞数 6

CC 4.0 BY-SA版权

文章标签：矩阵算法线性代数

于 2023-12-03 14:32:30 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_44220665/article/details/134755570

第一部分一维传输矩阵的推导

第一层到最后一层的总的传输矩阵，可以通过计算求得。

剩下的关键就是要确认材料的厚度、入射角度（假设为平行光垂直入射，角度取0）、以及每一层材料关于波长的折射率。

第二部分数据生成环节

先在FDTD软件中导出SiO2、Si3N4材料的折射率数据，再以.mat文件导入到matlab中

以TE模式举例：

2.1、加载数据

clc;
clear;
close all
SiO2 =load('SiO2.mat');
Si3N4 = load('Si3N4.mat');

2.2、定义、初始化

拟生成50000组数据

times = 50000;%指定模拟的迭代次数
data = cell(times,2);%初始化一个细胞数组，用于储存模拟结果

layer = 10;%定义结构的层数
d = ones(1,layer);%初始化一个数组，用于储存各层的厚度
S_temp = cell(1,layer);
P_temp = cell(1,layer);

2.3、计算（各层叠加后）总的传输矩阵

%计算每个波长的反射和透射系数
for t=1:times
theta = 0;
S = ones(2,2,length(SiO2.lambda));
P = ones(2,2,length(SiO2.lambda));
   
     %第2i-1层
     for j = 1:2:layer
        d(1,j) = fix(unifrnd(1,150))*1e-9;
 %调用function transmission_matrix_TE
      S_temp{1,j} = transmission_matrix_TE(SiO2.n(:,2),1,d(1,j),theta,SiO2.lambda);
     
     %第2i层
       d(1,j+1) = fix(unifrnd(1,150))*1e-9;

     S_temp{1,j+1} = transmission_matrix_TE(Si3N4.n(:,2),1,d(1,j+1),theta,Si3N4.lambda);
     end
     
    n0=ones(length(SiO2.lambda),1);   %相对真空折射率
    mum0=ones(length(SiO2.lambda),1) ;
    mum4=ones(length(SiO2.lambda),1) ; 
    eta0_s=n0./sqrt(mum0).*cos(theta);
    eta4_s=n0./sqrt(mum4).*cos(theta);
    
for i=1:length(SiO2.lambda)
    %总
   
    for n = 1:layer
          %S(:,:,i)=(S1(:,:,i)*S2(:,:,i))^10;将各层矩阵相乘，计算总的传输矩阵
    temp_S = S_temp{1,n};
    S(:,:,i) = (temp_S(:,:,i) * S(:,:,i));
    end
    r_TE(i)=((S(1,1,i)+S(1,2,i)*eta4_s(i))*eta0_s(i)-(S(2,1,i)+S(2,2,i)*eta4_s(i)))/((S(1,1,i)+S(1,2,i)*eta4_s(i))*eta0_s(i)+(S(2,1,i)+S(2,2,i)*eta4_s(i)));
    rc_TE(i)=conj(r_TE(i));
    Ref_TE(i)=r_TE(i)*rc_TE(i);
    t_TE(i)=(2*eta0_s(i))/((S(1,1,i)+S(1,2,i)*eta4_s(i))*eta0_s(i)+(S(2,1,i)+S(2,2,i)*eta4_s(i)));
    tc_TE(i)=conj(t_TE(i));
    T_TE(i)=(eta4_s(i)/eta0_s(i))*t_TE(i)*tc_TE(i);
end

2.4、上述 tramission_matrix_TE 的调用函数

function outputArg1 = transmission_matrix_TE(inputArg1,inputArg2,inputArg3,inputArg4,inputArg5)
index = inputArg1;
mum= inputArg2;
d=inputArg3;
theta=inputArg4;
lambda=inputArg5;

eta=index./sqrt(mum).*cos(theta);    
delte=2*pi.*index.*d./lambda.*cos(theta);
M=ones(2,2,length(lambda));
M(1,1,:)=cos(delte) ;   %传输矩阵A
M(1,2,:)=-1i*sin(delte)./eta;%1i：-1开根号    虚数   %传输矩阵B
M(2,1,:)=-eta.*1i.*sin(delte); %传输矩阵C
M(2,2,:)= cos(delte); %传输矩阵D

outputArg1=M;  %Mk 
end

2.5、计算并生成关于每个波长的反射系数、每层的厚度，保存到data数组中，以备后续神经网络训练

data{t,1} = Ref_TE;
data{t,2} = d;
end

第三部分：LSTM预测模型

3.1 导入数据 “load matlab” 直接调用上部分的data数据

clc;
clear;
close all
%% 导入数据
load  matlab

%% 训练数据
train_num=0.7;%训练数据占比
N=length(data);
for i=1:N*train_num
%   YTrain(i,1)=data{i,5};
%   YTrain(i,2)=data{i,6};
for j = 1:10
    temp = data{i,2};
  XTrain(i,j) = temp(:,j);
end
  for k=1:150
      train = data{i,1};
      YTrain(i,k)=train(:,k);
  end
  %YTrain{i,:}=train;
%   XTrain(i,:,:)=train';
end

3.2 测试数据分配

%% 测试数据
test_num=1-train_num;%测试数据占比
for i=1:N*test_num
%   YTest(i,1)=data{i,5};
%   YTest(i,2)=data{i,6};
for j = 1:10
    temp = data{7000+i,2};
  XTest(i,j) = temp(:,j);
end
for k=1:150
      test = data{7000+i,1};
      YTest(i,k)=test(:,k);
  end
%   XTest(i,:,:)=test';
end

3.3 数据归一化

XTrain=XTrain'*10^6;%训练输入
XTest=XTest'*10^6;%测试输入
YTrain=YTrain';%训练输出
YTest=YTest';%测试输出

3.4 参数设置

神经网络的学习率为 0.001。令收敛阈值设置为：RMSE 小于 0.5，并损失值小于 0.1

rng(0)
numFeatures = 10;%输入节点数
numResponses = 150;%输出节点数
miniBatchSize = 64; %batchsize
numHiddenUnits = 300;%隐含节点
maxEpochs=10000;%训练次数
learning_rate=0.001;%学习率
layers = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    lstmLayer(200)
    lstmLayer(200)
    reluLayer
    fullyConnectedLayer(numResponses)
    regressionLayer];
options = trainingOptions('adam', ...
    'ExecutionEnvironment','cpu',...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'InitialLearnRate',learning_rate, ...
    'GradientThreshold',1, ...
    'Shuffle','once', ...
    'Verbose',true,...
    'Plots','training-progress',...
    'ValidationData',{XTest YTest});
net = trainNetwork(XTrain,YTrain,layers,options);%建立神经网络模型
Pred = predict(net,XTrain,'ExecutionEnvironment','cpu');
Pred=double(Pred);
YPred = predict(net,XTest,'ExecutionEnvironment','cpu');
YPred=double(YPred);

3.5 绘图分析对比测试集、训练集中的原数据与预测数据的差异


figure(1)
plot(lambda,YTrain(:,12345))%1-35000
hold on
plot(lambda,Pred(:,12345))
grid on 
title('训练集LSTM')
legend('真实值','预测值')
xlabel('波长/米');ylabel('反射率')
figure(2)
plot(lambda,YTest(:,7500))
hold on
plot(lambda,YPred(:,7500))
grid on
title('测试集LSTM')
legend('真实值','预测值')
xlabel('波长/米');ylabel('反射率')
figure(3)
plot(lambda,YTest(:,1000))
hold on
plot(lambda,YPred(:,1000))
grid on
title('测试集LSTM')
legend('真实值','预测值')
xlabel('波长/米');ylabel('反射率')

3.6 效果对比

损失图：

测试集预测值与真实值差异：

训练集预测值与真实值差异：

第四部分反向神经网络模型建立

4.1 导入数据---分配测试集、训练集---数据归一化--建立神经网络

结构上大致与正向神经网络一致，反向的含义将输入输出端口互换

%% LSTM时间序列预测
clc;
clear;
close all
%% 导入数据
load  matlab
% plot(data{1,4})
%% 训练数据
train_num=0.7;%训练数据占比
N=length(data);
for i=1:N*train_num

for j = 1:10
    temp = data{i,2};
  YTrain(i,j) = temp(:,j);
end
  for k=1:150
      train = data{i,1};
      XTrain(i,k)=train(:,k);
  end
  %YTrain{i,:}=train;
%   XTrain(i,:,:)=train';
end
%% 测试数据
test_num=1-train_num;
for i=1:N*test_num
%   YTest(i,1)=data{i,5};
%   YTest(i,2)=data{i,6};
for j = 1:10
    temp = data{7000+i,2};
  YTest(i,j) = temp(:,j);
end
for k=1:150
      test = data{7000+i,1};
      XTest(i,k)=test(:,k);
  end
  XTest(i,:,:)=test';
end
%% 数据归一化

XTrain=XTrain';%训练输入
XTest=XTest';%测试输入
YTrain=YTrain'*10^7;%训练输出
YTest=YTest'*10^7;%测试输出

输入点变为150 即反射谱对应的150个点，输出点为10 代表十层材料各自的厚度

%% 参数设置
rng(0)
numFeatures = 150;%输入节点数
numResponses = 10;%输出节点数
miniBatchSize = 64; %batchsize
numHiddenUnits = 400;
maxEpochs=10000;
learning_rate=0.001;
layers = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    lstmLayer(100)
    reluLayer
    fullyConnectedLayer(numResponses)
    regressionLayer];
options = trainingOptions('adam', ...
    'ExecutionEnvironment','gpu',...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'InitialLearnRate',learning_rate, ...
    'GradientThreshold',1, ...
    'Shuffle','never', ...
    'Verbose',true,...
    'Plots','training-progress',...
    'ValidationData',{XTest YTest});
net = trainNetwork(XTrain,YTrain,layers,options);
Pred = predict(net,XTrain,'ExecutionEnvironment','cpu');
Pred=double(Pred);
YPred = predict(net,XTest,'ExecutionEnvironment','cpu');
YPred=double(YPred);

4.2 绘图分析

% % 绘图分析
figure(1)
plot(YTrain(:,35000)/10)%1-35000
hold on
plot(Pred(:,35000)/10)
grid on
title('训练集LSTM')
legend('真实值','预测值')
xlabel('层数');ylabel('厚度/微米')
figure(2)
plot(YTest(:,150)/10)%1-15000
hold on
plot(YPred(:,150)/10)
grid on
title('测试集LSTM')
legend('真实值','预测值')
xlabel('层数');ylabel('厚度/微米')