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

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

[Geoff Hayes]

Milad - what are you expecting to see? Or, what do you mean by guitar string won't converge? Have you tried increasing the number of timesteps? (What is the reason for choosing 300?)

[Milad Emadi]

Hi, thank you for your reply!

I haven't really considered what I was expecting, haha! I suppose I expected it to sort of oscillate a bit and then at some point reach a point where it stops flat in the resting position of the string (what I meant with converge), but what I am getting is this very wild, strange distortion after the first 15 iterations or so that seems to keep growing, it confused me a lot so I felt like there's definitely some mistake I've made somewhere along the way, even though I still haven't been able to find what I'm doing wrong.

I haven't tried increasing the timesteps to more than 300, I'll try that! I don't really have a reason for picking 300 - it just seemed to be long enough as the string seemed to see a fair bit of action during that time. I'll report back with what I see!

EDIT: I tried increasing the timestep to 1000 and seeing what happens - the values continute to grow wildly, by the end of 1000 iterations they have an order of magnitude of 10^32. This seems really strange to me, do you have any clue what the cause of these oscillations might be?

[Torsten]

Your boundary condition is only set for the first 5 times ...

[Milad Emadi]

Good point, you're right! The loops that later fill in the matrix are skipping position 1 and the final position in the matrix so they still become zeros I believe, I tried adjusting the matrix definition to match the amount of timesteps but it still didn't do what I was expecting it to, any other ideas?

[Torsten]

I think you will have to tell us first which partial differential equation you are trying to solve. The wave equation

d^2u/dt^2 - c^2*d^2u/dx^2 = 0

?

[Milad Emadi]

Yes exactly that one!

[Torsten]

And where can we find a link to the discretized form you use in your code ?

[Milad Emadi]

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]

From

http://hplgit.github.io/num-methods-for-PDEs/doc/pub/wave/html/._wave001.html#wave:string:fd

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

[Milad Emadi]

Oh what a silly mistake haha, thank you so much, I now get a nice smooth little simulation for the string, thanks!