Matlab, wave equation, and a flick of a string

[joshtanga]

Homework Statement

You know how if you flick one end of a garden hose you can watch the wave travel down, and then back to you? I wrote this MATLAB code to solve the associated PDE via the Fourier method and the resulting animation looks good for a few time steps, but then the solution curves are no longer smooth. It seems the Fourier coefficients are decaying too rapidly, but I can't tell where my code goes wrong.

Have a look, run the code, watch the animation and tell me why the curve suddenly changes character and looks so piecewise linear?

Homework Equations

To follow the code, you should be up to speed with the Fourier method for solving PDE's

The Attempt at a Solution

clear U
plotme = 1;
N = 200; % length of partial sum
c = 2000; % wave speed
l = 100; % length of guitar string is 1 meter
e = .01 % parameter of flick speed
t = 0:.001:.2; % parameterization of time
x = 1:l; % length of string

% The non-homogeneous Dirichlet condition on one endpoint
% This is a flick
f = e^(-1)*t -3*e^(-2)*t.^2 + 3*e^(-3)*t.^3 - e^(-4)*t.^4;
ix = find(t>e);
f(ix) = 0;

% The shifting function and the new forceing function, to make the Dirichlet condition homogeneous.
for i = 1:length(x)
p(:,i) = x(i) .* f;
g(:,i) = -x(i)/l*(-6*e^(-2)+18*e^(-3)*t - 12*e^(-4)*t.^2);
end
g(ix,:) = 1/100000*g(ix,:);


% The resulting forcing function on the pde
% gxt = -X/100*(-6*e^(-2)+18*e^(-3)*T - 12*e^(-4)*T^2)

% compute Fourier coeffs one by one
for n = 1:N

% Initial Velocity and a handy constant
th = c*n*pi/l;
du_0 = -x / (l*e);
dn = 2*(-1)^n / (n*pi*e); %2/l * trapz(x,du_0 .* sin(n*pi*x/l));
d_n(n) = dn; % verified

% This loop approximate the an(t) function (Fourier coef of shifted solution)
% at discrete values of t
for z = 1:length(t)

% FOURIER COEFFS of NEW FORCING FUNCTION
cnt = 2/l*trapz(x,g(z,:).*sin(n*pi*x/l));
CNT(z) = cnt;

% VARS FOR ALL TIME, AND A SPECIFIED TIME
tt = t(z);
s = linspace(0,tt,100);

% The gxt function is acutally a piecewise function. Dealt with here
ix = find(s>e);

% Computation of a_n(t)
%an = dn/th*(sin(th*tt) + trapz(s,sin(th*(tt-s)) .* cnt));
an = dn/th*(sin(th*tt) + trapz(s,sin(th*(tt-s)) .* cnt));

% Silly thing with 0
if an == 0
U(z,:,n) = zeros(size(x)) .* sin(n*pi*x/l); % RHS is a vector
else
U(z,:,n) = an * sin(n*pi*x/l); % RHS is a vector
end

% Save for later to look at
a_n(n,z) = an;
end

% display
My_n = n
end

% Computing the partial sum (summing U along the 2nd dimension)
u = p+sum(U,3); %partial sum, shifted

% Force the solution to satisfy boundary conditions
%u(:,1) = 0;
%u(:,l) = 0;

% Display the animation
if plotme == 1

% Good Animation Plotting Essentials
h = [];
A = [0 100 -10 10]; % set axis = [xmin xmax ymin ymax]
fig = figure; hold on; title('Solution to the Wave Equation'); axis(A);

% Create an empty .avi file
% A = avifile(AVI_FileName);

for i = 1:size(u,1)

delete(h);
h = plot(x,u(i,:),'.-'); % Index i ranges over time, thus: snapshots
drawnow;
% Display the current elaspsed time
disp(['Time = ',num2str((i))]);
pause(.1);

% Grab the current frame, so one can later watch it with "movie(F)"
% Fr(tt) = getframe;

% Extract the frame into the .avi file we created
% A = addframe(A,Fr(tt));
end

% Close the figure and the .avi file
%close(fig);
%A = close(A);

end % plotme

 

[Valkarie]

The solution curves suddenly become piecewise linear because the Fourier coefficients are decaying too rapidly. This is due to the fact that the Fourier method approximates solutions of the PDE by expanding them in a series of sine and cosine functions. As the number of terms in the series increases, the coefficients will decrease in magnitude, leading to less accurate approximations. In this case, the coefficients are decaying too rapidly, leading to the piecewise linear appearance of the solution curves.