%Code 1D equation String propagation FD method.
%DIEGO RODRIGUEZ #STUDENT 2014295001
clear all
clc
%define the medium
d1=1;
d2=4;
c1=3;
c2=1.5;
%discretized space
x=40; %input('Insert how much space do you want it? (x)(m): ')
dx=0.1; %input('Insert how much divide the space (dx)DEF (0.1): ');
nx=x/dx+1;
%discretized time
dt=0.009; %input('Insert how much divide the time (dt)DEF (0.01): ');
t=15; %input('Insert how much time do you want it? (t)(seg): ')
nt=t/dt+1;
%source conditions
f=1.6; %input('Insert how much frecuency do you want it? (Hz): ')
start=6.5; %input('Insert what space do you want to star it the signal? (m): ')
dts=0.0035;
T=1/f;
ts=0:dts:T
tao=0.07;
Wn=20;
A=exp(-(tao^2*Wn^2));
w1=10;
S=A*cos(w1*ts);
ns=start/dx+1;
U=zeros (nt,nx);
den=zeros (1,nx);
vel=zeros (1,nx);
den(1:round(nx/2)-1)=d1;
den(round(nx/2):nx)=d2;
vel(1:round(nx/2)-1)=c1;
vel(round(nx/2):nx)=c2;
Z=zeros (1,nx);
%Centered Method Second order
for k=2:nt-1
if k<=length(S);
U(k,ns)=U(k,ns)+S(k);
end
for i=2:nx-1
%Centered Method Second order
U(k+1,i)=(((den(i)*vel(i)^2*dt^2)/(den(i)*dx^2))*(U(k,i+1)-2*U(k,i)+U(k,i-1)))+2*U(k,i)-U(k-1,i);
end
end
%Centered Method with order 2M
%Coeficients order 2M product "a" is the matrix with the coeficients
for m=1:nx
Z=1;
for q=1:nx
if(q~=m)
Z=Z*((q^2-(vel(q)*dt/dx).^2)/(q^2-m^2));
end
end
a(m)=(((-1)^(m+1))/m^2)*Z;
end
a
%Solution wave propagation with 2M order
%for k=2:nt-1
% if k<=length(S);
% U(k,round(ns/100))=U(k,round(ns/100))+S(k);
% end
% for l=1:nx-1
% for qq=1:round(nx/2)-1
% for r=1:round(nt/2)-1
% for jh=1:nx
% U(round(nt/2)+r,round(nx/2))=((dt^2*vel(jh).^2/dx^2)*(a(1,nx/nx)*U(round(nt/2),round(nx/2))+sum(a(1,nx/nx+l)*(U(round(nt/2),round(nx/2)-qq)-U(round(nt/2),round(nx/2)+qq)))))+2*U(round(nt/2),round(nx/2))-U(round(nt/2)-r,round(nx/2));
% end
% end
% end
% end
%end
for k=1:nt
plot(U(k,:));
title('1D Wave Propargation')
xlabel('Space X')
ylabel('Amplitude')
getframe;
end
%Numerical dispersion
n=100;
nd=0:1:n;
phaseVel=ones(1,n);
phaseVel2=ones(1,n);
for z=1
for j=1:n
phaseVel(z,j)=(c1*cos(pi/j))/((1-(c1*dt/dx*sin(pi/j).^2)).^(1/2));
phaseVel2(z,j)=(c2*cos(pi/j))/((1-(c2*dt/dx*sin(pi/j).^2)).^(1/2));
end
end
for j=1:n
plot (real(phaseVel(j,:)),'b')
hold all
plot (real(phaseVel2(j,:)),'r')
xlabel('Number of Wavelength')
ylabel('Phase Velocity')
title('Numerical Dispersion')
text(30,3,'\leftarrow Vel 1','HorizontalAlignment','left')
text(30,1.5,'\leftarrow Vel 2','HorizontalAlignment','left')
end