Python Forum
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
A matrix of matrices
#11
I am coding it with numpy.
How are you going to code this matrix?
Speaking of which, I am rather new with Python (I used to use Scilab or Freefem++).
Is this language really suitable for maths? And with the extension numpy, does it really make a difference at the end of the day?
Reply
#12
I have not done it yet, because you coded it yourself.
I don't know Scilab or Freefem++.
Reply
#13
Well that's the whole point of the thread, I haven't been able to code this matrix...
Reply
#14
Sorry. I thought you have done it yourself.

(It is like image processing)
(replace a rectangle in an image)

If I have time tomorrow, I will do it.
Or if you have luck, someone is faster.
Or try to code it yourself.

Meanwhile you could read this or that:
https://stackoverflow.com/questions/3566...age-matrix
https://stackoverflow.com/questions/2650...x-in-numpy

Only I am curious:
How large can N be in your N²xN² matrix ?
Reply
#15
This matrix represents a square surface of side N². The utltimate goal of my exercise is to approximate the solution of the heat equation on this surface, and for this I am using the finite difference method. This whole part is actually rather easy. A couple of loops in and the job is done. Only thing is that, as I said, I am quite new at Python and there stil are stuff that confuse me (like first element of a matrix starts at 0 and not 1). So I'm getting there but it takes time. Python on top of it seems more rigid than the other languages I used to use.
Maybe in the end I only have syntax issues. I tried to concatenate the matrices. That sort of works but it only gives me a "line of matrices", and not a proper array. I also tried to think of it as a list of lists, which totally make sense, but as I said in my first post, I have no idea how to initialize such array; filling it would be no problem.
So maybe there are some more dedicated functions and tools in numpy that I'm not aware of. Thus my thread.
To answer your question, there is "theoretically" no limits on N. And given my exercise, it's going to be rather large.
Cheers!

And... Actually thanks to your help and to the links you posted, I finally found the solution!

u=-1*ones(Ns-1)
D=4*eye(Ns)+diag(u,1)+ diag(u,-1)
E=-eye(Ns)
N=Ns*Ns
A=array([0]*N*N).reshape(N,N)
for i in range (0,N,Ns):
      A[i:i+Ns,i:i+Ns]=D
      A[i:i+Ns,i+Ns:i+2*Ns]=E
      A[i+Ns:i+2*Ns,i:i+Ns]=E
A[N-Ns:N,N-Ns:N]=D
print('A = ', A)
For Ns=3:
Output:
A = [ [ 4 -1 0 -1 0 0 0 0 0] [-1 4 -1 0 -1 0 0 0 0] [ 0 -1 4 0 0 -1 0 0 0] [-1 0 0 4 -1 0 -1 0 0] [ 0 -1 0 -1 4 -1 0 -1 0] [ 0 0 -1 0 -1 4 0 0 -1] [ 0 0 0 -1 0 0 4 -1 0] [ 0 0 0 0 -1 0 -1 4 -1] [ 0 0 0 0 0 -1 0 -1 4]]
for Ns=4
Output:
A = [ [ 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0] [-1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0] [ 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0] [ 0 0 -1 4 0 0 0 -1 0 0 0 0 0 0 0 0] [-1 0 0 0 4 -1 0 0 -1 0 0 0 0 0 0 0] [ 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 0] [ 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0] [ 0 0 0 -1 0 0 -1 4 0 0 0 -1 0 0 0 0] [ 0 0 0 0 -1 0 0 0 4 -1 0 0 -1 0 0 0] [ 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0] [ 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0] [ 0 0 0 0 0 0 0 -1 0 0 -1 4 0 0 0 -1] [ 0 0 0 0 0 0 0 0 -1 0 0 0 4 -1 0 0] [ 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 0] [ 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1] [ 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4]]
Thanks a lot!
Reply
#16
Great work, but I get an error in line 8 with your program.
Also it would be good to post the whole program with imports,...

Meanwhile I corrected it myself (line 14):

#!/usr/bin/python3
from numpy import *

Ns= 4
u = -1*ones(Ns-1)
D =  4*eye(Ns) + diag(u,1) + diag(u,-1)
E = -1*eye(Ns)
print(D)
print()

N = Ns*Ns
A = array([0]*N*N).reshape(N,N)

for i in range(0,N-Ns,Ns):
    A[i:i+Ns,      i:i+Ns]      = D
    A[i:i+Ns,      i+Ns:i+2*Ns] = E
    A[i+Ns:i+2*Ns, i:i+Ns]      = E
A[N-Ns:N, N-Ns:N] = D

print(A)
Reply
#17
My largest Ns is 152, then I become a MemoryError.
I need more memory for my PC
Reply
#18
Since you helped me so much, you deserve my full code! haha
(I'm not sure that it's fully correct (that is, whether it actually solve the exercise I was trying to do), but it does work and produces something!)
Enjoy:
#Fichier TP3
#Méthode des différences finies 
#Equation de la chaleur

from numpy import *
from numpy.linalg import *
from matplotlib.pyplot import *
from mpl_toolkits.mplot3d import Axes3D

#Définition des fonctions
def uex(x,y):
    z=sin(2*pi*x)*sinh(2*pi*y)
    return z

def ud(x):
    y=sin(2*pi*x)*sinh(2*pi)
    return y

#Maillage 
print('Choix du nombre Ns de points intérieurs du maillage')
Ns = int(input('Ns = '))
h = 1/(Ns+2)


X=linspace(0,1,Ns+2)   #maillage complet


#Définition de la matrice A et du vecteur B
u=-1*ones(Ns-1)
D=4*eye(Ns)+diag(u,1)+ diag(u,-1)
E=-eye(Ns)
N=Ns*Ns
A=array([0]*N*N).reshape(N,N)
for i in range (0,N,Ns):
    for j in range (i,N-Ns,Ns):
        A[i:i+Ns,i:i+Ns]=D
        A[i:i+Ns,i+Ns:i+2*Ns]=E
        A[i+Ns:i+2*Ns,i:i+Ns]=E
A[N-Ns:N,N-Ns:N]=D  
print('A = ', A)

B=zeros((N,1))
C=zeros((Ns,1))
for i in range (0,Ns):
    C[i]=ud(X[i])
B[N-Ns:N]=C


#Définition de Uh
Uh=dot(inv(A),B)

#Définition de Zh
Zh=Uh.reshape(Ns,Ns)
print("Zh = ",Zh)       

#Rectification des bords de Zh
for i in range (0,Ns):
    for j in range (0,Ns):
        Zh[i,Ns-1]=ud(X[i])
        Zh[0,j]=0
        Zh[Ns-1,j]=0
        Zh[i,0]=0
        
print ('Zh = ',Zh)
        
#Définition de la solution exacte
Uex=zeros((Ns,Ns))
for i in range (0,Ns):
    for j in range (0,Ns):
        Uex[i,j]=uex(X[i],X[j])
print('Uex = ', Uex)

#Calcul de l'erreur
err=amax(Zh-Uex)
print("Erreur = ", err)

#Affichage des résultats        
coordX=zeros(Ns)
coordY=zeros(Ns) 
for k in range (0,Ns): 
    for i in range (0,Ns):
        coordX[k]=X[k] 
        coordY[i]=X[i]        
fig = figure() 
ax = Axes3D(fig) 
coordX, coordY = np.meshgrid(coordX,coordY)
#ax.plot_surface(coordX, coordY, abs(Zh-Uex), rstride=1, cstride=1, cmap='hot')
ax.plot_surface(coordX, coordY, Zh, rstride=1, cstride=1, cmap='hot')
#ax.plot_surface(coordX, coordY, Uex,rstride=1, cstride=1, cmap='hot')
show()
     
I let you see that output yourself!
Reply
#19
Thanks, nice plot.
For Ns=100, it looks like a sinus.
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  Estimating transition matrices in python Shreya10o 1 2,067 May-14-2020, 10:41 PM
Last Post: Larz60+
  How to optimize between two matrices, one of the matrices is derived from Excel data niloufar 0 1,630 Nov-05-2019, 06:28 PM
Last Post: niloufar
  Pandas dataframe: sum of exponentially weighted correlation matrices per row vvvcvvcv 1 3,206 May-29-2018, 01:09 AM
Last Post: scidam

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020