function string9 % wave equation for driven damped string, trapezoidal method close all; % global vect % vect = img2dat; strlength = 0.32; % length of string dx = 0.01; % size of x step n = round(strlength/dx); % number of x iterations tfin = 50; % length of time dt = 0.01; % size of t step, must be <= 0.02 m = round(tfin/dt); % number of t iterations % coefficients and other constant variables z = 0.00063; % linear mass density k = 52; % tension xf = 1/3*strlength; % position of finger on string xb = 1/2*strlength; % position of bow on string % damping due to finger fingvec = zeros(n-1,1); % initialize empty vector for finger damping fingloc = round(xf/dx); % index for location of damping on string fingvec(fingloc) = 1; % zero at all places except fingloc % bow action bowvec = zeros(n-1,1); % initialize empty vector for bow action bowloc = round(xb/dx); % index for location of bow on string bowvec(bowloc) = 1; % zero at all places except bowloc % initial conditions, t=0 x = (dx:dx:strlength-dx)'; % cut off ends of string for matching dimensions V = zeros(2*(n-1),m); V(1:n-1,1) = f(x); Q = ones(n-1,3); Q(:,2) = -2*Q(:,2); A = sparse(zeros(2*(n-1),2*(n-1))); A(1:n-1,n:2*(n-1)) = speye(n-1); A(n:2*(n-1),1:n-1) = spdiags(Q,[-1 0 1],n-1,n-1)*(k/z*dt^2/dx^2); A0 = A; A0(n:2*(n-1),n:2*(n-1)) = -spdiags(fingvec,0,n-1,n-1)*c(0)*dt/(dx*z); B = sparse(zeros(2*(n-1),1)); B0 = B; B0(n:2*(n-1),1) = b(0,dt)*bowvec*dt^2/(dx*z); for i=2:m t = (i-1)*dt; A(n-1+fingloc,n-1+fingloc) = -c(t)*dt/(dx*z); B(n-1+bowloc,1) = b(t,dt)*dt^2/(dx*z); D = speye(2*(n-1)) - dt/2*A; E = speye(2*(n-1)) + dt/2*A0; V(:,i) = D\(E*V(:,i-1) + dt/2*(B+B0)); % reset "previous" values for next iteration A0 = A; B0 = B; end V = [zeros(1,m); V(1:n-1,:); zeros(2,m); V(n:2*(n-1),:); zeros(1,m)]; % add rows of zeros for string's ends figure(1); rotate3d; % allow 3D plot to rotate pinc = 100; position = (0:dx:strlength)'*ones(1,m/pinc); time = ones(n+1,1)*(0:pinc:m-1)*dt; % length m/20 for a less dense plot vplot = V(1:n+1,1:pinc:end); hu = plot3(position,time,vplot,'Color',1/256*[139 131 120]); hold on; fingtime = dt:dt:tfin; fingpos = xf*ones(1,length(fingtime)); fingplot = V(fingloc,1:length(fingtime)); hf = plot3(fingpos,fingtime,fingplot,'Color',1/256*[238 203 173], ... 'LineWidth',3); bowtime = dt:dt:tfin; bowpos = xb*ones(1,length(bowtime)); bowplot = V(bowloc,1:length(bowtime)); hb = plot3(bowpos,bowtime,bowplot,'Color',1/256*[139 69 19],'LineWidth',2); grid on; legend([hu(1) hf(1) hb(1)],'String''s Path','Finger Damping','Bow Action'); title('String''s Path Over Time'); xlabel('position on string (m)'); ylabel('time (sec)'); zlabel('vertical displacement of string (m)'); axis([0 strlength 0 tfin]); % % sound % figure(2); % vsample = V(round((n+1)/2),1:end); % %vol=1; % increment to increase volume % %soundvec=reshape(vol*U',q*w,1); % reshape U matrix into a vector % plot(vsample); % soundsc(vsample,44000); % play string's sound % % movie % figure(3); % % find minimum and maximum vertical height of string to set axes % U = V(1:n+1,:)'; % [a1,b1]=size(U); % maxi=U(1,1); % mini=maxi; % for i=1:a1 % for j=1:b1 % if U(i,j)>maxi % maxi=U(i,j); % end % if U(i,j)tlen % out=0; % else % out=vect(round(t./tlen*veclen)); % end % end % % function vect=img2dat % % B=imread('wavdat.bmp'); % B=sparse(1-B); % A=imrotate(B,180); % [rows,cols]=size(A); % vect=zeros(cols,1); % for i=1:cols % tot=sum(A(:,i)); % if tot~=0 % vect(i)=mean(find(A(:,i))); % end % % end % vect=fliplr(vect); % vect=vect/max(vect); % % figure(4); % plot(vect); % % end