- #1

- 1

- 0

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