(adsbygoogle = window.adsbygoogle || []).push({}); 1. The problem statement, all variables and given/known data

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?

2. Relevant equations

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

3. 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

**Physics Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

Can you offer guidance or do you also need help?

**Physics Forums | Science Articles, Homework Help, Discussion**