Dismiss Notice
Join Physics Forums Today!
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

  1. Oct 28, 2011 #1
    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);
    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)

    cnt = 2/l*trapz(x,g(z,:).*sin(n*pi*x/l));
    CNT(z) = cnt;

    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
    U(z,:,n) = an * sin(n*pi*x/l); % RHS is a vector

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

    % display
    My_n = n

    % 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)

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

    % 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));

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

    end % plotme
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?