code to be optimised

Question

sriram amarnath 2019-6-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/469208-code-to-be-optimised

Code under %% 4000Hz takes too much time to run. Can you please help me optimise it for quicker results!

 %%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 (to be optimized)
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