sandbox/easystab/stab2014/vibrating_string_time_convergence
Convergence study of the time step
%{ Done by Alexandre Bilczewski and Fadil Boodoo. In this work, we do a convergence study of the time step of the code vibrating_string.m The scheme used for the march in time is Crank-Nicolson.
clear all; clf % parameters N=101; % number of gridpoints L=2pi;
% domain length c=1; % wave velocity % dt=0.5; % time step tmax=20; %
final time %Loop in time step for i=1:30 dt=i/10 ; % differentiation
matrices scale=-2/L; [x,DM] = chebdif(N,2);
dx=DM(:,:,1)scale;
dxx=DM(:,:,2)scale^2;
x=(x-1)/scale; Z=zeros(N,N); I=eye(N); II=eye(2N); % system
matrices E=[I,Z; Z,I]; A=[Z,c^2*dxx; I, Z]; % locations in the state
v=1:N; f=N+1:2N; % boundary conditions loc=[f(1),f(N)]; E(loc,:)=0;
A(loc,:)=II(loc,:); % march in time matrix
M=(E-Adt/2)(E+Adt/2); % initial condition q=[zeros(N,1);
sin(pi*x/L)]; % marching loop tvec=dt:dt:tmax; Nt=length(tvec);
e=zeros(Nt,1); for ind=1:Nt
q=Mq; % one step forward e(ind)=q(f((N+1)/2)); % store center point
end % time evolution of central point subplot(1,2,1); Ttheo=2L/c;
tt=linspace(0,tmax,500); etheo=cos(2pitt/Ttheo); % calculation
of the error etheoNt=cos(2pi*tvec/Ttheo); eetheo = etheoNt-e.’ ;
dtstock(i)=dt; error(i)=max(eetheo); end % plot the error xx2=dtstock;
plot(dtstock,error,‘b.-’); title(‘error in function of the time step’);
subplot(1,2,2); % comparison of the numerical error and the theorical
error xx=1./dtstock; loglog(xx,error,‘b.-’,xx,xx.^-2,‘r’);
legend(‘numerical error’,‘theorical error order 2’); title(‘Convergence
of the theorical error and the numerical error in loglog’) grid on
xlabel(‘dt’); ylabel(‘error’);
