Damped plucked string wave equation Fourier series assignment/ matrix dimensions problem

Hi there.

I'm having a little trouble implementing the solution to the damped 1-dimensional wave equation for a damped, plucked string. I would like to plot the position of a large number of points along the string (to represent the string itself), at discrete time instants.

My commented code is:

MATLAB code

L_s = 0.33; %   length of string (m)

m = 0.000125;    %   mass of string (kg)

T = 68; %   stretched tension of string (N)

R = 100;    %   linear resistance for string (Ns/m)

rho_L = m/L_s;  %   mass per unit length of string (kg/m)

c = sqrt(T/rho_L);  %   wave speed in string (m/s)

t = 1.3;    %   current time value (s)

d_0 = 0.2;    %   initial displacement value (m)

%   range variables

x = linspace(0,L_s,1000);   %   position points on string (m)

n = linspace(1,length(x),length(x));    %   integers for Fourier series

Omega_n = sqrt((c.*n.*pi./L_s).^2 - R^2);    %   argument coefficient for Fourier sine series

alpha_n = 8*d_0./(n.^2*pi^2).*sin(n.*pi./2);   %   first Fourier series coefficient

beta_n = R.*alpha_n./Omega_n; %   second Fourier series coefficient

w = zeros(1,length(n),length(x)); %   displacement function for string
for i = 1:length(n)
    for j = 1:length(x)
        w(i,j) = w + sin(n(i).*pi.*x(j)./L_s).*exp(-R./t).*(alpha_n(i).*cos(Omega_n(i).*t) ...
        + beta_n(i).*sin(Omega_n(i).*t));
    end
end

hold on

figure(1)
plot(x,w)
xlabel('Position on string (m)')
ylabel('Displacement (m)')
xlim([0 0.33])

The error message is: "Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in ass_assignment2 (line 39) w(i,j) = w + sin(n(i).*pi.*x(j)./L_s).*exp(-R./t).*(alpha_n(i).*cos(Omega_n(i).*t) ..."

Could anyone offer some advice? I understand the error message is telling me that the RHS comprises more dimensions than the LHS, but I cannot see why this would be.

Any help is much appreciated. Any general comments on more efficient and clearer coding techniques are also most welcome.

Regards,

Mike

Answer 1

[Jonathan Epperl]

In here

w(i,j) = w + sin(n(i).*pi.*x(j)./L_s).*exp(-R./t).*alpha_n(i).*cos(Omega_n(i).*t) ...
+ beta_n(i).*sin(Omega_n(i).*t));

the term w is a vector, as you defined it, the rest is all scalars, so what you will get is a vector with size size(w) with the scalars added to each element.

I'm also not sure why you initialize w as a 1 x length(n) x length(x) 3-D array, wouldn't a matrix be sufficient?

[Michael]

Yes, good point about the initialisation - thankyou. But alpha_n, Omega_n and beta_n are all 1x1000 vectors. If I redefine w as the 1000x1000 zeros matrix I get

"Subscripted assignment dimension mismatch."

What are I am trying to achieve is: for every position x along the 1000 points between 0 and L_s, calculate the value of each of the vectors x (for the value stored at the i-th position, alpha_n, beta_n and Omega_n (for their values stored at the j-th positions), combine these to make w (for each i-th and j-th position) as per the formula, then return to the start of the loop, and add this matrix to the next matrix formed by the next set of i-th and j-th positions.

[Michael]

I have found that I get closer to the right result with this slight modification:

%   base data

L_s = 0.33; %   length of string (m)
m = 0.000125;    %   mass of string (kg)
T = 68; %   stretched tension of string (N)
R = 100;    %   linear resistance for string (Ns/m)

rho_L = m/L_s;  %   mass per unit length of string (kg/m)

c = sqrt(T/rho_L);  %   wave speed in string (m/s)

t = 1;    %   current time value (s)
d_0 = 0.2;    %   initial displacement value (m)

%   range variables

x = linspace(0,L_s,1000);   %   position points on string (m)
n = linspace(1,length(x),length(x));    %   integers for Fourier series

Omega_n = sqrt((c.*n.*pi./L_s).^2 - R^2);    %   argument coefficient for Fourier sine series

alpha_n = 8*d_0./(n.^2*pi^2).*sin(n.*pi./2);   %   first Fourier series coefficient

beta_n = R.*alpha_n./Omega_n; %   second Fourier series coefficient

w = zeros(length(x),length(n)); %   displacement function for string
for i = 1:length(x)
    for j = 1:length(n)
        w(i,j) = w(i,j) + sin(n(j).*pi.*x(i)./L_s).*exp(-R./t).*(alpha_n(j).*cos(Omega_n(j).*t) ...
        + beta_n(j).*sin(Omega_n(j).*t));
    end
end

hold on
figure(1)
plot(x,w)
xlabel('Position on string (m)')
ylabel('Displacement (m)')
xlim([min(x) max(x)])

But this seems to be plotting all possible modes of displacement at maximum amplitude, rather than the particular string position at a time t given a certain initial displacement, which is what I'm after.

[Jonathan Epperl]

Hi, so if I understand correctly, then after you're done the value w(i,j) is the j-th Harmonic at x(i), and then to see your string, you'd have to change your call to plot:

plot(x,sum(w,2),'r')

Why the values are so tiny I don't know, but your R is so huge that exp(-R/t) ends up being tiny, maybe that's an error?

Two style-things:

• Your two nested for-loops can be replaced by

w2 = sin(pi/L_s* x' *n)* diag( alpha_n.*cos(Omega_n*t)+beta_n.*sin(Omega_n*t) ) *exp(-R/t);

Waaay faster, on my machine:

Loops:
Elapsed time is 0.381516 seconds.
Vectorized:
Elapsed time is 0.048705 seconds.

Also I don't understand what the addition of w(i,j) is supposed to do, since you initialize w to all zeros before?

• Instead of this

n = linspace(1,length(x),length(x));    %   integers for Fourier series

I recommend

n = 1:numel(x);

I'm not sure whether that's faster, but definitely more compact.

[Michael]

Thanks alot for the tips - much appreciated.

I managed to use the

sum(w,2)

function to combine all the string modes and give the position of the string at time t, as you suggest. The

exp(-R/t)

was indeed an error, it should be

exp(-Rt)

The

w(i,j) + ...

bit of the loop was unnecessary.

Sorted!

I'll look into use of the

diag

function to vectorise - it's new to me.

Thanks a lot for your help.

[Jonathan Epperl]

Glad I could help. Vectorization is probably the most speed-gaining technique in Matlab, you should definitely get familiar with it. Functions that I find often helpful:

ones(), diag(), blkdiag(), repmat(), meshgrid(), ndgrid()