MATLAB Answers

0

Finite Difference Method for a guitar string won't converge and I can't see why

Asked by Milad Emadi on 10 Mar 2019
Latest activity Commented on by Milad Emadi on 11 Mar 2019
Hi, I'm trying to simulate a guitar string as it is pulled up a bit and then released. I use the wave equation to derive the finite difference scheme (central difference in both time and space) and the equations I get are identical to what I've been able to find from googling around, so I assume I've done that part correctly. I also get the first "timestep" by using the initial velocity condition, this also seems to be alright. However, when I try to iterate the rest of the timesteps in the for-loop, it won't converge. It seems to go well at the first few iterations but it just goes completely haywire towards the end. I have ensured that the courant number is below 1, so it should be stable. I've looked through this quite carefully in my opinion, but I just can't see at all what I'm doing wrong. Thanks in advance to anyone who decides to help me out :)
clear all;
clc;
%Some parameters for the guitar string
T = 440;
%String diameter
D = 0.0005;
%String density
rho = 7800;
%String length
L = 0.63;
%Wave velocity on string
c = sqrt(T/(3.14*(D/2)^2*rho));
%Timestep
n = 0.00000118;
%Amount of space points
j = 100;
%Courant number, making sure s < 1
s = (c*n/(L/(j-1)))^2;
%Initial extension of guitar string vertically
uhat = 0.01;
%Initial extension of guitar string horisontally
xhat = 0.45;
%Boundary conditions
u = zeros(j, 5);
u(1,:) = 0;
u(j,:) = 0;
pos = 1;
%Create the initial position in matrix, in u(x,1)
for i = L/(j-1):L/(j-1):L-L/(j-1)
if i <= xhat
u(pos+1,1) = uhat*(i/xhat);
pos = pos+1;
end
if i > xhat && i < L
u(pos+1,1) = uhat*((i-L)/(xhat-L));
pos = pos+1;
end
end
xaxis = 0:L/(j-1):L;
%First iteration is slightly different
%Second position in approximation using the velocity initial condition
for pos = 2:1:j-1
u(pos, 2) = (s*(u(pos+1, 1)-u(pos-1, 1))+2*u(pos, 1)*(1-s))/2;
end
%Compute the rest of the positions while plotting them for 300 timesteps (I think this is where it goes wrong somewhere, but I have no idea why or where)
for t = 2:1:300
for pos = 2:1:j-1
u(pos, t+1) = s*(u(pos+1, t) - u(pos-1, t)) + 2*u(pos,t)*(1-s)-u(pos, t-1);
end
plot(xaxis,u(:,t-1));
axis([0 .7 -.015 .015]);
xlabel('Position on string')
ylabel('Extension of string vertically')
grid on;
drawnow;
end

  10 Comments

I approximate the second order derivatives with the explicit central difference formula:
Then I solve for and get the following recursion formula:
Where I say and then proceed to iterate it over the first timestep using the initial velocity condition to produce the initial condition , which allows me to do the first timestep. The rest of the timesteps are then looped with the recursion formula.
From
use formula (16) for the first timestep and formula (14) for the next timesteps.
Note that you have a sign error in your update formula for u(pos,t+1).
Oh what a silly mistake haha, thank you so much, I now get a nice smooth little simulation for the string, thanks!

Sign in to comment.

0 Answers