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

on 10 Mar 2019

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

on 11 Mar 2019
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.
Torsten

### Torsten (view profile)

on 11 Mar 2019
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).