Finite Difference Method for a guitar string won't converge and I can't see why
5 views (last 30 days)
Show older comments
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
Answers (0)
See Also
Categories
Find more on DSP Algorithm Acceleration in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!