Is there a way of removing these for loops to speed up my code?
9 views (last 30 days)
Show older comments
Matthew Hunt
on 16 Oct 2019
Answered: Fabio Freschi
on 18 Oct 2019
I have the analytical solution to the following PDE:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242960/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242961/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242962/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242963/image.png)
in the form of a Green's function solution:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242964/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242965/image.png)
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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242966/image.png)
2 Comments
Fabio Freschi
on 17 Oct 2019
Can you also provide a suitable set of f,mu_n,r,t,Gamma,D,u_0 to run and profile the code?
Accepted Answer
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](https://www.mathworks.com/matlabcentral/answers/uploaded_files/243294/Screen%20Shot%202019-10-18%20at%2008.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](https://www.mathworks.com/matlabcentral/answers/uploaded_files/243295/Screen%20Shot%202019-10-18%20at%2008.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!
0 Comments
More Answers (0)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!