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

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?

[Fabio Freschi]

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]

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.

Answer 1

[Fabio Freschi]

Profiling your code it comes up that the most demanding routine is the Green function calculation, with 0.374s:

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:

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!