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

12 次查看（过去 30 天）
Matthew Hunt
Matthew Hunt 2019-10-16
回答： Fabio Freschi 2019-10-18
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 个评论
Fabio Freschi
Fabio Freschi 2019-10-17
Can you also provide a suitable set of f,mu_n,r,t,Gamma,D,u_0 to run and profile the code?
Matthew Hunt
Matthew Hunt 2019-10-17
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Fabio Freschi
Fabio Freschi 2019-10-18
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!

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by