Reduction of repeated lines in code

Need help making the code more readable by reducing the repeated lines(dfunc lines) in section 1b:

%% 1a velocity 0-1kHz
Lx=0.286; % length of plate(m)
Lz=0.198; % breadth of plate(m)
h=3.25e-3;% thickness of plate(m)(also d)
E=1.4e9 ;%Youngs modulus
M=265*h;%Mass per unit area( rowS surface density)
mu=0.3;%Poisson?s ratio
eta=4e-2;%Damping Factor
x0=0.16;%exc=(0.16,0.083), point of excitation(m)
z0=0.083;%exc=(0.16,0.083), point of excitation(m)
F=1;% Point load assumed to be 1 newton @ (x0,z0)
f=linspace(1000,1,10); % frequency array
w=2*pi.*f;% angular frequency array
t=1; %  time 
m=1:1000;% mode numbers width
n=10; % mode numbers length
func=zeros(n,length(m));% initial values as place holders
for i=1:n 
func(i,:)=(i*pi./Lx).^2 +(m*pi/Lz).^2; % modes for m=1:1000 & n=1:10
end
D=E*h^3/(12*(1-mu^2));% flexural density D eq 86 Notes
wnm=sqrt(D/M)*func; % for loop for wnm is not required as it ia a scalar multiple of func
 
% for all frequencies common vfunc is defined. vfunc is multi-dimensional
% array with dimensions as 10,10,1000. 
vfunc=zeros(length(w),n,length(m));% 3D array initial values as place holders
for k=1:length(w)
for i=1:n
vfunc(k,i,:)= 1j*w(k)*((sin(i*pi*x0/Lx).*sin(m.*pi*z0./Lz)).^2./(-w(k)'.^2 ...
                        +wnm(i,:).^2.*(1+1j*eta))).*exp(1j.*w(k).*t);          % velocity function  
end
end
v=(4/(Lx*Lz))*(F/M)*vfunc; %diffrentiated Eq 98 Notes
vSum=squeeze(sum(v,2)); % squeeze-removes the singleton dimension
RlvSum=sum(real(vSum),2);% vector of all modes at chosen freq array
figure('NumberTitle','off','Name','Velocity 0-1kHz')
plot(f,RlvSum);
xlabel('Frequency Hz')
ylabel('Velocity m/s')
 
 
 
%% 1b Displacement plot 
xr=linspace(0,Lx,1000);% 1000 x co-ordinate points on surface
zr=linspace(0,Lz,1000); % 1000  z co-ordinate points on surface
wp=2*pi*500; % pressure plot @ 500 Hz
[X,Z]=meshgrid(xr,zr);% plane to be plotted on 
dfunc1=(sin(pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(1,:).*(1+1j*eta));% mode n=1 & m=1:1000
dfunc2=(sin(2.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((2.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(2,:).*(1+1j*eta));% mode n=2 & m=1:1000
dfunc3=(sin(3.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((3.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(3,:).*(1+1j*eta));% mode n=3 & m=1:1000
dfunc4=(sin(4.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((4.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(4,:).*(1+1j*eta));% mode n=4 & m=1:1000
dfunc5=(sin(5.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((5.*pi.*X))/Lx).*sin(m.*pi.*Z/Lz).*exp(1j*wp*t)./(-wp^2+wnm(5,:).*(1+1j*eta));% mode n=5 & m=1:1000
dfunc6=(sin(6.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((6.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(6,:).*(1+1j*eta));% mode n=6 & m=1:1000
dfunc7=(sin(7.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((7.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(7,:).*(1+1j*eta));% mode n=7 & m=1:1000
dfunc8=(sin(8.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((8.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(8,:).*(1+1j*eta));% mode n=8 & m=1:1000
dfunc9=(sin(9.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((9.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(9,:).*(1+1j*eta));% mode n=9 & m=1:1000
dfunc10=(sin(10.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((10.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(10,:).*(1+1j*eta));% mode n=10 & m=1:1000
dt1=dfunc1+dfunc2+dfunc3+dfunc4+dfunc5+dfunc6+dfunc7+dfunc8+dfunc9+dfunc10;% total all modes
ydt=(4/(Lx*Lz))*(F/M).*real(dt1);% displacement due to all modes
figure('NumberTitle','off','Name','Displacment @ 500 Hz')
surf(X,Z,ydt,'EdgeColor','interp');
view(-90,90)
xlabel('Length')
zlabel('Displacement')
ylabel('Width')
colorbar

Answer 1

[KSSV]

This should do:

dfunc = zeros(1000,1000,10) ; 
for n = 1:10
    dfunc(:,:,i) = (sin(2.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((n.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(n,:).*(1+1j*eta));
end

The whole code:

clc; clear all ;
%% 1a velocity 0-1kHz
Lx=0.286; % length of plate(m)
Lz=0.198; % breadth of plate(m)
h=3.25e-3;% thickness of plate(m)(also d)
E=1.4e9 ;%Youngs modulus
M=265*h;%Mass per unit area( rowS surface density)
mu=0.3;%Poisson?s ratio
eta=4e-2;%Damping Factor
x0=0.16;%exc=(0.16,0.083), point of excitation(m)
z0=0.083;%exc=(0.16,0.083), point of excitation(m)
F=1;% Point load assumed to be 1 newton @ (x0,z0)
f=linspace(1000,1,10); % frequency array
w=2*pi.*f;% angular frequency array
t=1; %  time 
m=1:1000;% mode numbers width
n=10; % mode numbers length
func=zeros(n,length(m));% initial values as place holders
for i=1:n 
func(i,:)=(i*pi./Lx).^2 +(m*pi/Lz).^2; % modes for m=1:1000 & n=1:10
end
D=E*h^3/(12*(1-mu^2));% flexural density D eq 86 Notes
wnm=sqrt(D/M)*func; % for loop for wnm is not required as it ia a scalar multiple of func
 
% for all frequencies common vfunc is defined. vfunc is multi-dimensional
% array with dimensions as 10,10,1000. 
vfunc=zeros(length(w),n,length(m));% 3D array initial values as place holders
for k=1:length(w)
for i=1:n
vfunc(k,i,:)= 1j*w(k)*((sin(i*pi*x0/Lx).*sin(m.*pi*z0./Lz)).^2./(-w(k)'.^2 ...
                        +wnm(i,:).^2.*(1+1j*eta))).*exp(1j.*w(k).*t);          % velocity function  
end
end
v=(4/(Lx*Lz))*(F/M)*vfunc; %diffrentiated Eq 98 Notes
vSum=squeeze(sum(v,2)); % squeeze-removes the singleton dimension
RlvSum=sum(real(vSum),2);% vector of all modes at chosen freq array
figure('NumberTitle','off','Name','Velocity 0-1kHz')
plot(f,RlvSum);
xlabel('Frequency Hz')
ylabel('Velocity m/s')
 
 
 
%% 1b Displacement plot 
xr=linspace(0,Lx,1000);% 1000 x co-ordinate points on surface
zr=linspace(0,Lz,1000); % 1000  z co-ordinate points on surface
wp=2*pi*500; % pressure plot @ 500 Hz
[X,Z]=meshgrid(xr,zr);% plane to be plotted on 
dfunc = zeros(1000,1000,10) ; 
for n = 1:10
    dfunc(:,:,i) = (sin(2.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((n.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(n,:).*(1+1j*eta));
end
dt1 = sum(dfunc,3) ; 
ydt=(4/(Lx*Lz))*(F/M).*real(dt1);% displacement due to all modes
figure('NumberTitle','off','Name','Displacment @ 500 Hz')
surf(X,Z,ydt,'EdgeColor','interp');
view(-90,90)
xlabel('Length')
zlabel('Displacement')
ylabel('Width')
colorbar

[sriram amarnath]

Thank you for the answer can you also help me in refining the code under %% 4000 Hz it takes too much time to run it. Can you please optimize it

1a velocity 0-1kHz

clc; clear all;

Lx=0.286; % length of plate(m)

Lz=0.198; % breadth of plate(m)

h=3.25e-3;% thickness of plate(m)(also d)

E=1.4e9 ;%Youngs modulus

M=265*h;%Mass per unit area(rowS surface density)

mu=0.3;%Poisson?s ratio

eta=4e-2;%Damping Factor

x0=0.16;%exc=(0.16,0.083), point of excitation(m)

z0=0.083;%exc=(0.16,0.083), point of excitation(m)

F=1;% Point load assumed to be 1 newton @ (x0,z0)

f=linspace(1000,1,10); % frequency array

w=2*pi.*f;% angular frequency array

t=1; % time

m=1:1000;% mode numbers

fv=10;

func=zeros(fv,length(m));

for i=1:fv

func(i,:)=(i*pi./Lx).^2 +(m*pi/Lz).^2; %

end

% func1=(1*pi./Lx).^2 +(m.*pi./Lz).^2;% summation function n=1 & m= 1:1000

%%% I have used matrix approach to minimize your code. Your func1, func2,

%%% func3, func4 etc., are row vectors of func matrix in line 21.

D=E*h^3/(12*(1-mu^2));% flexural density D eq 86 Notes

wnm=sqrt(D/M)*func; % for loop for wnm is not required as it ia a scalar multiple of func

% for all frequencies common vfunc is defined. vfunc is multi-dimensional

% array with dimensions as 10,10,1000.

vfunc=zeros(length(w),fv,length(m));

for k=1:length(w)

for i=1:fv

vfunc(k,i,:)= 1j*w(k)*((sin(i*pi*x0/Lx).*sin(m.*pi*z0./Lz)).^2./(-w(k)'.^2 ...

+wnm(i,:).^2.*(1+1j*eta))).*exp(1j.*w(k).*t); % velocity function

end

end

%vfunc2= 1j.*w(1,1)'.*((sin(2*pi*x0/Lx).*sin(m.*pi*z0./Lz)).^2./(-w(1,1)'.^2+wnm2.^2.*(1+1j*eta))).*exp(1j.*w(1,1).*t');% velocity function

v=(4/(Lx*Lz))*(F/M)*vfunc; %diffrentiated Eq 98 Notes

%v2=(4/(Lx*Lz))*(F/M).*vfunc2;% diffrentiated Eq 98 Notes

%v1000= sum(v,1);% total mode contributions at 1kHz

vSum=squeeze(sum(v,2)); % squeeze-removes the signleton dimension

%sv1=sum(real(v1000)); % real part of all modes

RlvSum=sum(real(vSum),2);

% sv=[sv1 sv2 sv3 sv4 sv5 sv6 sv7 sv8 sv9 sv10];% vector of all modes

figure('NumberTitle','off','Name','Velocity 0-1kHz')

plot(f,RlvSum);

xlabel('Frequency Hz')

ylabel('Velocity m/s')

% I didn't go further, you can use a similar approach.

%% 1b Displacement plot

xr=linspace(0,Lx,1000);% 1000 x co-ordinate points on surface

zr=linspace(0,Lz,1000); % 1000 z co-ordinate points on surface

wp=2*pi*500; % pressure plot @ 500 Hz

[X,Z]=meshgrid(xr,zr);% plane to be plotted on

dfunc = zeros(1000,1000,10) ;

for n = 1:10

dfunc(:,:,i) = (sin(2.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((n.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(n,:).*(1+1j*eta));

end

dt1 = sum(dfunc,3) ;

ydt=(4/(Lx*Lz))*(F/M).*real(dt1);% displacement due to all modes

figure('NumberTitle','off','Name','Displacment @ 500 Hz')

surf(X,Z,ydt,'EdgeColor','interp');

view(-90,90)

xlabel('Length')

zlabel('Displacement')

ylabel('Width')

colorbar

%% Shape function

combo=[2 2; 2 3; 5 3]; % (n=2,m=2)using eq 90 & 91 notes % (n=2,m=3)using eq 90 & 91 notes % (n=5,m=3)using eq 90 & 91 notes

psi=zeros(3,length(xr),length(zr));

freqs=[112 445 850];

%Vsh=zeros(length(m),3);

for i=1:length(combo(:,1))

psi(i,:,:)=sin((combo(i,1)*pi.*X)/Lx).*sin((combo(i,2)*pi.*Z)./Lz);

figure

surf(X,Z,squeeze(psi(i,:,:)),'EdgeColor','interp');

xlabel('Length')

zlabel('Displacement')

ylabel('Width')

colorbar

title(['Shape Function',num2str(freqs(i)), 'Hz'])

Vsh(:,i)=getframe(gcf);% save figure to a video frame

end

close all

%% 1c Between Modes

movie(Vsh,2) % video of 3 resonance modes

Msh=VideoWriter('Resonance modes','Uncompressed AVI');

Msh.FrameRate=1; % frame rate

open(Msh);

writeVideo(Msh,Vsh);% writing video

close(Msh);

%% 4000 Hz

wp4k=2*pi*4000; % driving frequency of 4000Hz

dfunc2=zeros(length(xr),length(zr));

tic

for n=1:fv % index k is equiv to n in act. eqn.

for k=1:length(m) % k equiv to m in eqn. note--> m(k)=k

dfunc2=sin(n*pi*x0/Lx).*sin(k*pi*z0/Lz)...

*sin((n*pi*X)/Lx).*sin(k*pi*Z/Lz)...

.*exp(1j*wp*t)./(-wp4k^2+wnm(n,k).^2+2*1j*eta)+dfunc2; % note--> size of wnm is 10x1000 => its indices should be (k,n)

end

end

toc

ydt4k=(4/(Lx*Lz))*(F/M).*real(dfunc2);% displacement of all modes

figure('NumberTitle','off','Name','Driving Frequency 4000 Hz')

surf(X,Z,ydt4k,'EdgeColor','interp');

xlabel('Length')

zlabel('Displacement')

ylabel('Width')

colorbar