1.程序功能描述
基于离散差分法的复杂微分方程组求解.“连续微分方程”到“离散微分方程”到“差分方程”，离散微分方程，变成差分方程。建立差分方程时，时间采用一阶显格式，空间采用一阶偏心差分格式。
2.测试软件版本以及运行结果展示
MATLAB2022a版本运行
3.核心程序
```% ʼ
L = 0.05; % ռ䳤
time = 1e-8; %ʱ 䳤
Nz = 10; % ռ
Nt = 10; %ʱ
dt = time/Nt;%t ֵ ۼ
dz = L/Nz;%z ֵ ۼ
N1 = zeros(Nz,Nt);
N2 = zeros(Nz,Nt);
N3 = zeros(Nz,Nt);
N4 = zeros(Nz,Nt);
N1_Yb = zeros(Nz,Nt);
N2_Yb = zeros(Nz,Nt);
Ps = zeros(Nz,Nt);
PASE_plus = zeros(M,Nz,Nt);
PASE_minus = zeros(M,Nz,Nt);
Pp_plus = zeros(Nz,Nt);
Pp_minus = zeros(Nz,Nt);
G = zeros(Nz,Nt);
NF = zeros(Nz,Nt);
% 1 ʽ 1 ϵ IJ ʾ
W11 = FpO13_vp/(AchVp);
W12 = FsO12_vs/(AchVs);
for i = 1:M
W13(i) = F_ASE_vj(i) O12_vj(i) / (Ach*Vj(i));
end
W14 = FsO21_vs/(AchVs);
for i = 1:M
W15(i) = F_ASE_vj(i) O21_vj(i) / (AchVj(i));
end
W16 = FpO31_vp/(Ach*Vp);
% 1 ʽ 2 ϵ IJ ʾ
W21 = FsO12_vs/(Ach*Vs);
for i = 1:M
W22(i) = F_ASE_vj(i) O12_vj(i) / (AchVj(i));
end
W23 = FsO21_vs/(AchVs);
for i = 1:M
W24(i) = F_ASE_vj(i) O21_vj(i) / (Ach*Vj(i));
end
% 1 ʽ 3 ϵ IJ ʾ
W31 = FpO13_vp/(AchVp);
W32 = FpO31_vp/(AchVp);
% 1 ʽ 4 ϵ IJ ʾ
W41 = FpO12Yb_vp/(AchVp);
W42 = FpO21Yb_vp/(AchVp);
Ps(1,:) = 0.001ones(1,Nt);
Pp_plus(1,:) = 0.06ones(1,Nt);
Pp_minus(1,:) = 0.04*ones(1,Nt);
for j = 1:Nt-1
for i = 1:Nz-1
% 1ʽ 1
N1(i,j+1) = N1(i,j) + ...
dt( -1(W11(Pp_plus(i,j) + Pp_minus(i,j)) + W12Ps(i,j) + sum(W13.(PASE_plus(:,i,j)+PASE_minus(:,i,j))'))N1(i,j) +...
(A21 + W14Ps(i,j) + sum(W15.(PASE_plus(:,i,j)+PASE_minus(:,i,j))'))N2(i,j) + ...
C2N2(i,j)^2 + W16(Pp_plus(i,j) + Pp_minus(i,j))N3(i,j) + C3N3(i,j)^2 - C14N1(i,j)N4(i,j)+...
-1KtrN2_Yb(i,j)N1(i,j)+KbaN1_Yb(i,j)N3(i,j) );
% 1ʽ 2
N2(i,j+1) = N2(i,j) + ...
dt*( (W21*Ps(i,j)+sum(W22.*(PASE_plus(:,i,j)+PASE_minus(:,i,j))'))*N1(i,j) +...
-1*(A21 + W23*Ps(i,j) + sum( W24.*(PASE_plus(:,i,j)+PASE_minus(:,i,j))' ))*N2(i,j) +...
A32*N3(i,j) - 2*C2*N2(i,j)^2 + 2*C14*N1(i,j)*N4(i,j) );
% 1ʽ 3
N3(i,j+1) = N3(i,j) + ...
dt*( W31*(Pp_plus(i,j) + Pp_minus(i,j))*N1(i,j) - A32*N3(i,j) - W32*(Pp_plus(i,j) + Pp_minus(i,j))*N3(i,j) -...
2*C3*N3(i,j)^2 + A43*N4(i,j) + Ktr*N2_Yb(i,j)*N1(i,j) - Kba*N1_Yb(i,j)*N3(i,j) );
% 1ʽ 4
N1_Yb(i,j+1) = N1_Yb(i,j) + ...
dt*(-1*W41*(Pp_plus(i,j) + Pp_minus(i,j))*N1_Yb(i,j) + W42*(Pp_plus(i,j) + Pp_minus(i,j))*N2_Yb(i,j) +...
A21Yb*N2_Yb(i,j) + Ktr*N2_Yb(i,j)*N1(i,j) - Kba*N1_Yb(i,j)*N3(i,j));
% 1ʽ 5
N4(i,j+1) = NEr - (N1(i,j+1) + N2(i,j+1) + N3(i,j+1));
% 1ʽ 6
N2_Yb(i,j+1) = NYb - N1_Yb(i,j+1);
if N1(i,j+1) > NEr,N1(i,j+1) = NEr;end
if N2(i,j+1) > NEr,N2(i,j+1) = NEr;end
if N3(i,j+1) > NEr,N3(i,j+1) = NEr;end
if N4(i,j+1) > NEr,N4(i,j+1) = NEr;end
if N1_Yb(i,j+1) > NYb,N1_Yb(i,j+1) = NYb;end
if N2_Yb(i,j+1) > NYb,N2_Yb(i,j+1) = NYb;end
if N1(i,j+1) < 0,N1(i,j+1) = 0;end
if N2(i,j+1) < 0,N2(i,j+1) = 0;end
if N3(i,j+1) < 0,N3(i,j+1) = 0;end
if N4(i,j+1) < 0,N4(i,j+1) = 0;end
if N1_Yb(i,j+1) < 0,N1_Yb(i,j+1) = 0;end
if N2_Yb(i,j+1) < 0,N2_Yb(i,j+1) = 0;end
% Ϸ ̼ õ N1 N2 N3 N4 N1Yb N2Yb
% 2
Pp_plus(i+1,j) = Pp_plus(i,j) + dz*(-Fp*(O13_vp*N1(i,j) - O31_vp*N3(i,j) + O12Yb_vp*N1_Yb(i,j) - O21Yb_vp*N2_Yb(i,j))*Pp_plus(i,j) - ap*Pp_plus(i,j));
Pp_minus(i+1,j) = Pp_minus(i,j) + dz*(Fp*(O13_vp*N1(i,j) - O31_vp*N3(i,j) + O12Yb_vp*N1_Yb(i,j) - O21Yb_vp*N2_Yb(i,j))*Pp_minus(i,j) + ap*Pp_plus(i,j));
Ps(i+1,j) = Ps(i,j) + dz*(Fs*( O21_vs*N2(i,j) - O12_vs*N1(i,j) )*Ps(i,j) - as*Ps(i,j));
for ii = 1:M
PASE_plus(ii,i+1,j) = PASE_plus(ii,i,j)+dz*(F_ASE_vj(ii)*( O21_vj(ii)*N2(i,j) - O12_vj(ii)*N1(i,j) ) * PASE_plus(ii,i,j) +...
2*h*Vj(ii)*DVj(ii)*F_ASE_vj(ii)*O21_vj(ii)*N2(i,j)-as*PASE_plus(ii,i,j));
PASE_minus(ii,i+1,j) = PASE_minus(ii,i,j)+dz*(-1*F_ASE_vj(ii)*( O21_vj(ii)*N2(i,j) - O12_vj(ii)*N1(i,j) ) * PASE_minus(ii,i,j) -...
2*h*Vj(ii)*DVj(ii)*F_ASE_vj(ii)*O21_vj(ii)*N2(i,j)+as*PASE_minus(ii,i,j));
end
if Pp_plus(i+1,j) < 0,Pp_plus(i+1,j) = 0;end
if Pp_minus(i+1,j) < 0,Pp_minus(i+1,j) = 0;end
if Ps(i+1,j) < 0,Ps(i+1,j) = 0;end
%ͨ ̬ õ Pp+ Pp- Pase+ Pase- Ps
end
end
for z = 1:Nz
for t = 1:Nt
PASE_plus2(z,t) = sum(PASE_plus(:,z,t));
PASE_minus2(z,t) = sum(PASE_minus(:,z,t));
end
end
for z = 1:Nz
for t = 1:Nt
G(z,t) = 10*log10(Ps(z,t)/Ps(1,1));
end
end
for z = 1:Nz
for t = 1:Nt
NF(z,t) = 10log10(1/G(z,t) + PASE_plus2(z,t)/(G(z,t)Vs*DVs) );
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pp_plus2 = interp1(dz:dz:L,Pp_plus(1:end,Nz),0:dz/10:L,'cubic');
Pp_minus2 = interp1(dz:dz:L,Pp_minus(1:end,Nz),0:dz/10:L,'cubic');
figure;
subplot(211);
plot(0:dz/10:L,Pp_plus2,'g-','LineWidth',3);
xlabel('z');
ylabel('Pp+(Z)');
title('Pp+(Z)&z');
grid on;
subplot(212);
plot(L:-dz/10:0,Pp_minus2,'m--','LineWidth',2);
xlabel('z');
ylabel('Pp-(Z)');
title('Pp-(Z)&z');
grid on;
16_015m
```
4.本算法原理
本课题求解的方程组表达式如下：
基于离散差分法的复杂微分方程组求解.“连续微分方程”到“离散微分方程”到“差分方程”，离散微分方程，变成差分方程。建立差分方程时，时间采用一阶显格式，空间采用一阶偏心差分格式。
