# Matlab, wave equation, and a flick of a string

## 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