Improve code for the numerical aproximation of the 1D fractional heat equation

I need a way to improve the time it takes to calculate the matrix. Right now it takes around 19 seconds for the small interval [-100,100] and h = .125. and I need it to work for h = e-6 and and the interval [a,b] = [-5000,5000]

tic
%Bound 
a = -100;
b = 100;
T = 1;

%Alpha value
alpha = 1;


h = 0.125;
tau = h^alpha/100;

%Boundary conditions
f1 = @(x) 1./(1+x.^2);

%Exact solution
y = @(x,t) (t+1)./((t+1).^2 + x.^2);

%Draw of the aproximation and the error
DFNLE(alpha,a,b,T,h,tau,f1,y)

With the main function being

function DFNLE(alpha,a,b,T,h,tau,f1,y)

n = (b-a)/h - 1;
n = floor(n);
x = a:h:b;
m = T / tau;
m = nearest(m);
t = 0:tau:T;

w = zeros(n+2,1);
v1 = zeros(n,1);
v2 = zeros(n,1);

%Restrain
theta = h^(1/3);

%Constant C
C = 2^alpha * gamma(1/2+alpha/2) / (sqrt(pi)*abs(gamma(-alpha/2)));

%Known values of u
u = zeros(n+2,m+1);
u(:,1) = f1(x);


%Loop to create the weights
 for k = 1:n+2
     if k*h > theta
        w(k) = C*h^(-alpha)/k^(alpha+1);
     end
 end
 

%CFL condition
if tau > 1/sum(w)
   fprintf('Not convergent tau = %s', tau);
   return
end

% Main loop
for j = 1:m
    for i = 1:n+2
        for k = 1:n+2
 
            if i+k <= n+2 
               v1(k) = (u(i+k,j)-u(i,j))*w(k);
            else
               v1(k) = -u(i,j)*w(k);
            
            end
                              
            if  i-k >= 1
               v2(k) = (u(i-k,j)-u(i,j))*w(k);   
            else 
               v2(k) = -u(i,j)*w(k);
            end
                           
        end  
        u(i,j+1) = u(i,j)+tau*(sum(v1)+sum(v2));  
    end
end

final = T;
[X,T] = meshgrid(x,t);

u = sparse(u);

%Error
err = max(max(abs(u' - y(X,T))));
fprintf('Error, epsilon = %d', err)


% Drawing
% surf(X,T,u');
% axis([-10 10 0 final 0 1])
% shading interp       
% title(['Explicit method for alpha  = ', num2str(alpha)])
% xlabel('x')
% ylabel('t')
% zlabel('function u') 

end

Answer 1

[Shadaab Siddiqie]

From my understanding you want to improve the performance of above code. MATLAB is optimized to work on matrices and vectorizing the code improves MATLAB performance compared to using nested FOR loops.

The following two examples use the functions ARRAYFUN and CELLFUN to compute the same final matrix without using FOR loops. When the size of the matrix is small, all these approaches are similar in terms of the time needed to execute them. However, when the matrix size increases, these approaches are much faster.

%% 1st approach: using arrays and ARRAYFUN
% Step 1: Compute the cumulative variance vector (last column of the target matrix Y).
tic
f = @(k) var(Z(1:k)); 
w = arrayfun(f, (1:numel(Z)).');
% Step 2: Assemble the matrix (which is built from zeros
% and subvectors of the cumulative variance vector).
Y3 = triu( w*ones(1, numel(w)) );
toc
%% 2nd approach: using cell arrays and CELLFUN
A = Z';
n = length(A);
tic
Y2=triu(repmat(cellfun(@(x) var(x(isfinite(x))), num2cell(repmat(A,n,1) + triu(NaN(n, n),1),2)),1,n));
toc
%% Comparison
isequal(Y, Y2, Y3)

You can incorporate these functions and vectorize your code to improve it's performance.