Posts: 9
Threads: 1
Joined: Nov 2017
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?
Posts: 606
Threads: 3
Joined: Nov 2016
Nov-26-2017, 04:31 AM
(This post was last modified: Nov-26-2017, 04:49 AM by heiner55.)
I have not done it yet, because you coded it yourself.
I don't know Scilab or Freefem++.
Posts: 9
Threads: 1
Joined: Nov 2017
Well that's the whole point of the thread, I haven't been able to code this matrix...
Posts: 606
Threads: 3
Joined: Nov 2016
Nov-26-2017, 05:24 AM
(This post was last modified: Nov-26-2017, 05:24 AM by heiner55.)
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 ?
Posts: 9
Threads: 1
Joined: Nov 2017
Nov-26-2017, 08:13 AM
(This post was last modified: Nov-26-2017, 08:13 AM by MaximeDt.)
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!
Posts: 606
Threads: 3
Joined: Nov 2016
Nov-26-2017, 09:03 AM
(This post was last modified: Nov-26-2017, 09:03 AM by heiner55.)
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)
Posts: 606
Threads: 3
Joined: Nov 2016
Nov-26-2017, 12:50 PM
(This post was last modified: Nov-26-2017, 12:50 PM by heiner55.)
My largest Ns is 152, then I become a MemoryError.
I need more memory for my PC
Posts: 9
Threads: 1
Joined: Nov 2017
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!
Posts: 606
Threads: 3
Joined: Nov 2016
Nov-26-2017, 01:23 PM
(This post was last modified: Nov-26-2017, 01:23 PM by heiner55.)
Thanks, nice plot.
For Ns=100, it looks like a sinus.
|