1
\$\begingroup\$

I am solving a simple 1D steady heat question using spectral method. I am a long time MATLAB and Mathematica user and trying to learn Python. I compare only linear solving times and Python is way slower than MATLAB which doesn't feel right. Can you please tell me why Python is an order of magnitude slower for only even linear solve?

MATLAB code

  N=5000;
% HERE I AM JUST CREATING THE MATRIX FOR d/dx
 Dhat=zeros(N+1,N+1);

 for j=1:N;
 for i=mod(j-1,2):2:j;
 Dhat(i+1,j+1)=2*j;
 end
 end
 Dhat(1,:)=Dhat(1,:)/2;
 cbar=[2; ones(N-1,1); 2];
 p=0:1:N; 
 pn=cos((pi/N)*kron(p',p));
 I=ones(N+1,1);
 II=(-1).^(0:N);
 G=kron(I,II);
 x=-cos(pi*(0:N)/N)'; 
 T=G.*pn;
 Tinv=(2/N)*(G'./kron(cbar,cbar')).*pn;

 D=T*Dhat*Tinv;     % derivative matrix
 % We have d/dx matrix

 D2=D*D;            % d2/dx2 to solve d^2T/dx^2 = F

 T0= 0;             % Temp at left boundary
 Tn= 10;            % Temp at right boundary
 F= zeros(N+1,1);   % source term
 F(1)=0;            % Boundary condition
  F(end)=10;         % Boundary condition

 D2(1,:)=zeros;   
 D2(end,:)=zeros;
 D2(1,1)=1;          % At the left boundary T = 0 
 D2(end,end)=1;      % At the right boundary T = 10
 tic
 T=D2/F';
 toc

Python code

def chebgl(N):
    import numpy as np
    cbar = np.array([1] * (N - 1))
    cbar = np.insert(cbar, 0, 2)
    cbar = np.insert(cbar, N, 2)

    p = np.array(range(N + 1))
    pkp = np.kron(p, p)
    pkp = pkp.reshape(N + 1, N + 1)
    pn = np.cos((np.pi / N) * pkp)

    I = np.ones(N + 1)
    Ia = -1 * np.ones(N + 1)
    II = np.power(Ia, range(N + 1))
    G = np.kron(I, II)
    G = G.reshape(N + 1, N + 1)

    T = np.multiply(G, pn)
    GTrans = np.transpose(G)

    cbarK = np.kron(cbar, cbar)
    cbarK = cbarK.reshape(N + 1, N + 1)

    Tinv = (2 / N) * np.multiply((np.divide(GTrans, cbarK)), pn)

    dhat = np.zeros((N + 1, N + 1))

    for j in range(1, N + 1, 1):

         for i in range((j - 1) % 2, j, 2):
             dhat[i, j] = 2 * j

    dhat[0] = dhat[0] / 2
     # This is the operator for d/dx

     D = np.matmul(np.matmul(T, dhat), Tinv)  
     # predefined x-locations
     x = -np.cos(np.pi * (p / N))

     return D, x

Main code

import numpy as np
import chebGL
import matplotlib.pyplot as plt
from time import process_time

# Because I want to see all double precision digits
np.set_printoptions(precision=15)
# Number of points (d^2/dx^2) T = F(x) is going to be solved  
N=5000  

D, x = chebGL.chebgl(N)

T0 = 0        # Left boundary condition
Tend = 10     # Right boundary condition
F = 0 * np.ones(N + 1)  # Source Term
F[0] = T0     # inputting the boundary condition
F[N] = Tend  # inputting the boundary condition


D2T = np.matmul(D,D)     # Creating d^2/dx^2 operator
D2T[0] = abs(0 * D2T[0])
D2T[N] = abs(0 * D2T[0])

D2T[0,0] = 1     # inputting the left boundary condition
D2T[N,N] = 1    # inputting the right boundary condition

t1_start = process_time()

T = np.linalg.solve(D2T, F)
t1_stop = process_time()

print("Elapsed time during the whole program in seconds:", t1_stop- 
t1_start)
\$\endgroup\$
  • 1
    \$\begingroup\$ Thank you, I edit my question. \$\endgroup\$ – Erdem Sep 15 at 15:26
  • 1
    \$\begingroup\$ @AJNeufeld No problem, if you can reproduce it, perhaps you can take it to meta :) But the saga continues. Now OP is editing from older revision. Oh boy, I bet this question will end up with more edits than wished for. \$\endgroup\$ – dfhwze Sep 15 at 15:40
  • \$\begingroup\$ I'm not sure this post asks for a code review. \$\endgroup\$ – IEatBagels Sep 16 at 14:29
  • 1
    \$\begingroup\$ @IEatBagels Depending on which revision you look at, it is. It's simply phrased a little iffy, but that's ok. There have been enough edits. Consider it a performance question, we've done plenty of those. \$\endgroup\$ – Mast Sep 16 at 18:25