目录
一、理论基础
差分迭代法是一种数值求解离散微分方程的方法,其基本原理是通过将微分方程转化为差分方程,并利用迭代的方式进行求解。
1.1、差分迭代法的基本原理
假设我们要求解如下形式的离散微分方程:
y(n+1) = f(n, y(n))
其中,y(n) 表示变量 y 在时刻 n 的取值,f(n, y(n)) 是一个关于 n 和 y(n) 的函数。
差分迭代法的基本思想是利用差分方程来逼近微分方程,从而将微分方程的求解问题转化为差分方程的求解问题。具体地,我们可以将微分方程转化为如下形式的差分方程:
y(n+1) = y(n) + h * f(n, y(n))
其中,h 表示时间步长。
然后,我们可以利用迭代的方式来求解差分方程。假设我们已经知道 y(0) 的值,那么可以通过以下步骤来迭代求解 y(n) 的值:
- 将 y(0) 代入差分方程,计算得到 y(1) 的值。
- 将 y(1) 代入差分方程,计算得到 y(2) 的值。
- 以此类推,直到计算得到所需的 y(n) 的值。
1.2、差分迭代法的数学公式
假设我们要求解如下形式的离散微分方程:
y'(t) = f(t, y(t))
其中,y'(t) 表示变量 y 在时刻 t 的导数,f(t, y(t)) 是一个关于 t 和 y(t) 的函数。
为了将微分方程转化为差分方程,我们需要引入时间步长 h,并将 t 离散化为 t = nh,其中 n 表示时刻。然后,我们可以将微分方程转化为如下形式的差分方程:
y((n+1)h) = y(nh) + h * f(nh, y(nh))
其中,y((n+1)h) 表示变量 y 在时刻 (n+1)h 的取值。
接下来,我们可以利用迭代的方式来求解差分方程。假设我们已经知道 y(0) 的值,那么可以通过以下步骤来迭代求解 y((n+1)h) 的值:
- 将 y(0) 代入差分方程,计算得到 y(h) 的值。
- 将 y(h) 代入差分方程,计算得到 y(2h) 的值。
- 以此类推,直到计算得到所需的 y((n+1)h) 的值。
需要注意的是,差分迭代法的精度取决于时间步长 h 的选取。如果 h 过大,会导致差分方程与微分方程的误差较大;如果 h 过小,会导致计算量增加,同时机器精度也可能影响计算结果的准确性。因此,在实际应用中需要根据具体情况选择合适的时间步长 h。
1.3 复杂方程组的差分方程
这类复杂方程组,其解方程通常不能直接使用MATLAB内部的自带函数进行求解,往往需要使用别的算法进行,本课题涉及到的方程组的基本类型如下所示:
“连续微分方程”到“离散微分方程”到“差分方程”,离散微分方程,变成差分方程。建立差分方程时,时间采用一阶显格式,空间采用一阶偏心差分格式。
二、MATLAB仿真程序
clc;
clear;
close all;
warning off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vj = 1.0e14*[2.037671762107052
2.034215151823580
2.030770248941575
2.027336994082841
2.023915328270042
2.020505192923336
2.017106529857023
2.013719281276238
2.010343389773680
2.006978798326360
2.003625450292398
2.000283289407840
1.996952259783514
1.993632305901912
1.990323372614108
1.987025405136702
1.983738349048801
1.980462150289017
1.977196755152515
1.973942110288066
1.970698162695152
1.967464859721083
1.964242149058149
1.961029978740801
1.957828297142857
1.954637052974735
1.951456195280716
1.948285673436231
1.945125437145175
1.941975436437247
1.938835621665319
1.935705943502825
1.932586352941177
1.929476801287208
1.926377240160643
1.923287621491580
1.920207897518015
1.917138020783374
1.914077944134078
1.911027620717132
1.907987003977725
1.904956047656871
1.901934705789056
1.898922932699921];
DVj = 1.0e+011*[3.462483877836962
3.450746652796573
3.439069007521719
3.427450539446898
3.415890849400914
3.404389541572597
3.392946223476911
3.381560505921475
3.370232002973479
3.358960331926962
3.347745113270506
3.336585970655279
3.325482530863469
3.314434423777078
3.303441282347068
3.292502742562887
3.281618443422334
3.270788026901764
3.260011137926652
3.249287424342494
3.238616536886035
3.227998129156823
3.217431857589106
3.206917381424041
3.196454362682216
3.186042466136488
3.175681359285135
3.165370712325313
3.155110198126804
3.144899492206068
3.134738272700597
3.124626220343543
3.114563018438641
3.104548352835411
3.094581911904646
3.084663386514161
3.074792470004828
3.064968858166864
3.055192249216406
3.045462343772321
3.035778844833293
3.026141457755156
3.016549890228479
3.007003852256406];
F_ASE_vj = [0.727687625579220
0.726884730352238
0.726081696237599
0.725278527912882
0.724475230042399
0.723671807275757
0.722868264249905
0.722064605587050
0.721260835896694
0.720456959773481
0.719652981799073
0.718848906541101
0.718044738553185
0.717240482375278
0.716436142533662
0.715631723540323
0.714827229893890
0.714022666078659
0.713218036565208
0.712413345810112
0.711608598256263
0.710803798332059
0.709998950452820
0.709194059019136
0.708389128418137
0.707584163022211
0.706779167191096
0.705974145268722
0.705169101586599
0.704364040461595
0.703558966196364
0.702753883079915
0.701948795386920
0.701143707377818
0.700338623299606
0.699533547384828
0.698728483851725
0.697923436905154
0.697118410735281
0.696313409518148
0.695508437416143
0.694703498577095
0.693898597135282
0.693093737210495];
O12_vj = 1.0e-24*[ 0.11
0.13
0.175
0.215
0.3
0.275
0.288
0.31
0.33
0.34
0.335
0.33
0.335
0.33
0.303
0.28
0.25
0.25
0.248
0.26
0.268
0.28
0.325
0.425
0.57
0.628
0.53
0.426
0.388
0.35
0.31
0.25
0.24
0.22
0.195
0.17
0.15
0.14
0.118
0.1
0.078
0.068
0.06
0.059];
O21_vj = 1.0e-24*[ 0.05
0.06
0.07
0.095
0.11
0.13
0.15
0.16
0.18
0.2
0.21
0.218
0.23
0.231
0.23
0.228
0.217
0.218
0.225
0.246
0.275
0.318
0.378
0.51
0.76
0.89
0.725
0.638
0.618
0.58
0.527
0.485
0.49
0.475
0.42
0.39
0.365
0.346
0.305
0.262
0.221
0.195
0.17
0.169];
%参数初始化
Vp = 3.074794441025641e14;
Vs = 1.955593333333334e14;
DVs = 3.186042466136488e11;
h = 6.626e-34;
Ac = 10.5625e-12;
Fp = 0.8773;
Fs = 0.7078;
NEr = 1.5e26;
NYb = 1.9e27;
O12_vs = 6.5e-25;
O21_vs = 9.0e-25;
O13_vp = 2.58e-25;
O31_vp = 0;
O12Yb_vp= 1.0e-24;
O21Yb_vp= 1.0e-24;
r21 = 7.9e-3; A21 = 1/r21;
r32 = 1.0e-9; A32 = 1/r32;
r43 = 1.0e-9; A43 = 1/r43;
r21Yb = 2.0e-3; A21Yb = 1/r21Yb;
C2 = 5.0e-23;
C3 = 5.0e-23;
C14 = 3.5e-23;
Ktr = 2.0e-22;
Kba = 0;
ap = 0;
as = 0;
M = 44;
%其余参数初始化
L = 0.05; %空间长度
time = 1e-8; %时间长度
Nz = 100; %空间点数
Nt = 100; %时间网点数
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复杂系数的参数表示
W11 = Fp*O13_vp/(Ac*h*Vp);
W12 = Fs*O12_vs/(Ac*h*Vs);
for i = 1:M
W13(i) = F_ASE_vj(i) * O12_vj(i) / (Ac*h*Vj(i));
end
W14 = Fs*O21_vs/(Ac*h*Vs);
for i = 1:M
W15(i) = F_ASE_vj(i) * O21_vj(i) / (Ac*h*Vj(i));
end
W16 = Fp*O31_vp/(Ac*h*Vp);
%方程组1的式子2复杂系数的参数表示
W21 = Fs*O12_vs/(Ac*h*Vs);
for i = 1:M
W22(i) = F_ASE_vj(i) * O21_vj(i) / (Ac*h*Vj(i));
end
W23 = Fs*O21_vs/(Ac*h*Vs);
for i = 1:M
W24(i) = F_ASE_vj(i) * O21_vj(i) / (Ac*h*Vj(i));
end
%方程组1的式子3复杂系数的参数表示
W31 = Fp*O13_vp/(Ac*h*Vp);
W32 = Fp*O31_vp/(Ac*h*Vp);
%方程组1的式子4复杂系数的参数表示
W41 = Fp*O12Yb_vp/(Ac*h*Vp);
W42 = Fp*O21Yb_vp/(Ac*h*Vp);
Ps(1,:) = ones(1,Nt);
Pp_plus(1,:) = ones(1,Nt);
Pp_minus(1,:) = 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)) + W12*Ps(i,j) + sum(W13.*(PASE_plus(:,i,j)+PASE_minus(:,i,j))'))*N1(i,j) +...
(A21 + W14*Ps(i,j) + sum(W15.*(PASE_plus(:,i,j)+PASE_minus(:,i,j))'))*N2(i,j) + ...
C2*N2(i,j)^2 + W16*(Pp_plus(i,j) + Pp_minus(i,j))*N3(i,j) + C3*N3(i,j)^2 - C14*N1(i,j)*N4(i,j)+...
-1*Ktr*N2_Yb(i,j)*N1(i,j)+Kba*N1_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));
%PASE_plus = zeros(M,Nz,Nt);
%PASE_minus = zeros(M,Nz,Nt);
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) = 10*log10(1/G(z,t) + PASE_plus2(z,t)/(G(z,t)*Vs*DVs) );
end
end
%计算得到N1~t,N2~t,N3~t,N4~t,........
figure;
subplot(221);
plot(N1(1,2:end),'b-','LineWidth',2);
xlabel('t');
ylabel('N1(Z)');
title('N1(Z)&t');
grid on;
subplot(222);
plot(N2(1,2:end),'b-','LineWidth',2);
xlabel('t');
ylabel('N2(Z)');
title('N2(Z)&t');
grid on;
subplot(223);
plot(N3(1,2:end),'b-','LineWidth',2);
xlabel('t');
ylabel('N3(Z)');
title('N3(Z)&t');
grid on;
subplot(224);
plot(N4(1,2:end),'b-','LineWidth',2);
xlabel('t');
ylabel('N4(Z)');
title('N4(Z)&t');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
subplot(211);
plot(N1_Yb(1,2:end),'b-','LineWidth',2);
xlabel('t');
ylabel('N1Yb(Z)');
title('N1Yb(Z)&t');
grid on;
subplot(212);
plot(N2_Yb(1,2:end),'b-','LineWidth',2);
xlabel('t');
ylabel('N2Yb(Z)');
title('N2Yb(Z)&t');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
subplot(211);
plot(PASE_plus2(3,1:end-2),'b-','LineWidth',2);
xlabel('t');
ylabel('PASE+(Z)');
title('PASE+(Z)&t');
grid on;
subplot(212);
plot(PASE_minus2(3,1:end-2),'b-','LineWidth',2);
xlabel('t');
ylabel('PASE-(Z)');
title('PASE-(Z)&t');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
plot(Ps(2,1:end-1),'b-','LineWidth',2);
xlabel('t');
ylabel('Ps(Z)');
title('Ps(Z)&t');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
subplot(211);
plot(Pp_plus(20,3:end-1),'b-','LineWidth',2);
xlabel('t');
ylabel('Pp+(Z)');
title('Pp+(Z)&t');
grid on;
subplot(212);
plot(Pp_minus(20,3:end-1),'b-','LineWidth',2);
xlabel('t');
ylabel('Pp-(Z)');
title('Pp-(Z)&t');
grid on;
figure;
subplot(211);
plot(Pp_plus(1:end,6),'b-','LineWidth',2);
xlabel('t');
ylabel('Pp+(Z)');
title('Pp+(Z)&t');
grid on;
subplot(212);
plot(Pp_minus(1:end,6),'b-','LineWidth',2);
xlabel('t');
ylabel('Pp-(Z)');
title('Pp-(Z)&t');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
plot(G(20,3:end-1),'b-','LineWidth',2);
xlabel('t');
ylabel('G(Z)dB');
title('G(Z)&t');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
plot(NF(2,3:end-1),'b-','LineWidth',2);
xlabel('t');
ylabel('NF(Z)dB');
title('NF(Z)&t');
grid on;
三、仿真结果
A16-15