4.4 Program 4.1, Matlab Code
61
%dx∗ = 8nm/(2N) = 0.04nm.
%dx = dx ∗ /L = 0.04nm/8nm = 5x10 − 3 dimensionless dx.
x0=20.; N=100; tmax=100; dt=0.111; dx=5.e-3; D=9.e-4; Force=0.; Vo=10.;
sigma=4.;
ka = (1/2) ∗ D ∗ (dt/dx 2 );
ka1 = V o ∗ dx 2 ;
ka2=Vo*dx/2;
k12dt=6.e-3;
k21dt=6.e-3;
k12=5.376e-2;
k21=5.376e-2;
s2 = 1/(2 ∗ (sigma 2 ));
s0=1./(sigma*sqrt(2*pi));
T=floor(tmax/dt)+1; % T é o num. de pontos temporais
x=((1:2*N)’)*dx*8;
f = s0 ∗ exp(−((x − x0). 2 ) ∗ s2);
p1=f;
F=zeros(2*N,9);
E=zeros(2*N);E1=zeros(2*N); E2=zeros(2*N); M=zeros(2*N);M1=zeros(2*N);
M1=zeros(2*N);K12=zeros(2*N); K21=zeros(2*N); L1=zeros(2*N); L2=zeros(2*N);
KR21=zeros(2*N);KL12=zeros(2*N);
%STATE 1
for j=1:N
d21=(8*pi)*(-sin(2*pi*j*dx)+2*sin(4*pi*j*dx)-. . .
3*sin(6*pi*j*dx));
a1=1-ka*(-2 + ka1*d21);
b1=2-a1;
K12(j,j)=k12*dt;
E1(j,j)=a1;
L1(j,j)=b1;
end
%STATE 2
for j=N+1:2*N
d22= 0.;
a2=1-ka*(-2 + ka1*d22);
b2=2-a2;
K21(j,j)=k21*dt;
E2(j,j)=a2;
L2(j,j)=b2;
end
%STATE 1, DIAGONAL NEIGBOURGH
for j=2:N
d11=4.*(cos(2.*pi*j*dx)-cos(4.*pi*j*dx)+cos(6.*pi*j*dx))-Force;
c1=ka*(1.+ka2*d11);
d1=ka*(1.-ka2*d11);