3.7 Programs
43
d2=-Vo*(2*pi)*(-sin(2*pi*j*dx) + sin(4*pi*j*dx));%dimensionless derivative
a=1-k*(d2*(dx 2 )-2.);
b=2-a;
E(j,j)=a;
M(j,j)=b;
end
for j=2:L;
d1=-Vo*(cos(2.*pi*j*dx) - 0.5*cos(4.*pi*j*dx))-F;%dimensionless derivative
c=k*(1.+ 0.5*dx*d1);
d=2.*k - c;
E(j,j-1)=-d;
E(j-1,j)=-c;
M(j,j-1)=d;
M(j-1,j)=c;
end
IE=inv(E);
for tp=1:8
for it=1:pt
p=IE*(M*p);
P(:,tp)=p; %distribution at times tp
end
end
figure(1); plot(x,P(:,1),x,P(:,2),x,P(:,3),x,P(:,4), . . .
x,P(:,5),x,P(:,6),x,P(:,7),x,P(:,8))
F1=P(:,1)
F2=P(:,2)
F3=P(:,3)
F4=P(:,4)
F5=P(:,5)
F6=P(:,6)
F7=P(:,7)
F8=P(:,8)
F9=p1’
save F1.dat F1 -ascii;
save F2.dat F2 -ascii;
save F3.dat F3 -ascii;
save F4.dat F4 -ascii;
save F5.dat F5 -ascii;
save F6.dat F6 -ascii;
save F7.dat F7 -ascii;
save F8.dat F8 -ascii;
save F9.dat F9 -ascii;
save x.dat x -ascii;