MATLAB Answers

Is there a way of removing these for loops to speed up my code?

38 views (last 30 days)
I have the analytical solution to the following PDE:
with boundary condition and and initial condition
in the form of a Green's function solution:
Where is an infinite series summing over n. I have written code with lots of for loops and I was wondering if there was a way I could speed up my code? My code for the Greens function is:
function G = Greens_fn(D,mu_n,r,t,r_bar)
%This is the Green's function for spherical diffusion equation
%Solution can be found in A.D. Polyanin. The Handbook of Linear Partial Differential Equations
%for Engineers and Scientists, Chapman & Hall, CRC 17
G=zeros(length(t),length(r_bar));
N_t=length(t);
N=length(mu_n); %Number of terms used in Green's function series
for i=1:N
for j=1:N_t
G_i(j,:)=(2*r_bar/r).*((1+mu_n(i)^-2).*sin(mu_n(i)*r).*sin(mu_n(i)*r_bar).*exp(-D*mu_n(i)^2*t(j))); %Ther series part of the Greens function
G(j,:)=G(j,:)+G_i(j,:); %Adding uo the terms of the series
end
end
for i=1:N_t
G(i,:)=G(i,:)+3*r_bar.^2; %Adding the final part to give the Green's function
end
end
My code for the solution is:
function sol = u(f,mu_n,r,t,Gamma,D,u_0)
sol=zeros(length(t),1);
sol(1)=u_0;
%This is the solution of the spherical diffusion equation using Greens
%Function.
if (length(u_0)==1)
N_r=100; %Steps used in the r_bar integration
else
N_r=length(u_0);
end
N_t=length(f); %Steps used in the tau integration.
r_bar=linspace(0,1,N_r);
G_1=Greens_fn(D,mu_n,r,t,r_bar);
for i=2:N_t
G_2=Greens_fn(D,mu_n,1,t(i)-t(1:i),1);
sol(i)=trapz(r_bar,u_0.*G_1(i,:))+(D*Gamma)*trapz(t(1:i),f(1:i).*G_2);
end
end
I should point out that I'm only interested in . Is there a way I can speed this up?
  2 Comments
Matthew Hunt
Matthew Hunt on 17 Oct 2019
So some values which will help you to get an idea are:
N_t=300;
t=linspace(0,1,N_t);
r=1;
f=5*ones(N_t,1);
D=3.403;
Gamma=0.0183;
u_0=0.98;
N=20;
old=0;error=10^-8;
mu_n=zeros(N,1); %array which stores the solutions
for i=1:N
err=10;
while err>error %This solves the iteration x_n+1=atan(x_n)+n*pi
y_n=atan(old)+i*pi;
err=abs(tan(y_n)-y_n);
old=y_n;
end
mu_n(i)=y_n;
end
That should be some realistic values for you to get working with.

Sign in to comment.

Accepted Answer

Fabio Freschi
Fabio Freschi on 18 Oct 2019
Profiling your code it comes up that the most demanding routine is the Green function calculation, with 0.374s:
Screen Shot 2019-10-18 at 08.33.25.png
I have rewritten your function, eliminating the innermost for loop and vectorizing the code. I made explicit use of bsxfun even if, starting from Matlab2016b, it is possible to use implicit expansion. I prefer to clearly have under control what is happening in the code.
I also removed the last loop to update the final part of the greens function. The result is the following. I renamed the function as Greens_fn2 to avoid confusion (Remember to change the call in your u function).
function G = Greens_fn2(D,mu_n,r,t,r_bar)
%This is the Green's function for spherical diffusion equation
%Solution can be found in A.D. Polyanin. The Handbook of Linear Partial Differential Equations
%for Engineers and Scientists, Chapman & Hall, CRC 17
G=zeros(length(t),length(r_bar));
N=length(mu_n); %Number of terms used in Green's function series
for i=1:N
G = G+bsxfun(@times,bsxfun(@times,(1+mu_n(i)^-2).*sin(mu_n(i)*r).*sin(mu_n(i)*r_bar),exp(-D*mu_n(i)^2*t(:))),(2*r_bar/r));
end
G = bsxfun(@plus,G,3*r_bar.^2);
end
A new profile, gives the new performance time:
Screen Shot 2019-10-18 at 08.33.44.png
that is 0.094s, with a speedup of about 4. The results are identical.
Possible additional measures to speedup the code:
  • I am not sure but maybe you can also expand with respect to N_t, removing also the remaining loop.
  • Use a mex-file for the calculation of the Green's function
Befor any other attempt, I would analyze how this code scales with respect to N_r and N_t (I think you are interested in more complex cases). Please let us know!

More Answers (0)

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!