- Code: Tout sélectionner
# taille du domaine [0,L]x[0,H]
L=2.;
H=3.;
#nombre de points de discretisation et pas de maillage
N=100;
M=180;
#
l=L/(N+1);
h=H/(M+1);
# vecteur discretisation selon x
x=linspace(l,L-l,N)
# vecteur discretisation selon y
y=linspace(h,H-h,M)
#construction de la matrice A par ses diagonales
#
D0=(2/(l**2)+2/(h**2))*ones(N*M)# diagonale principale
D1=-1/l**2*ones(N*M)# surdiagonale
D1[N::N]=0.#correction de la surdiagonale (voisin de droite n existe pas au bord droit)
DM1=-1/l**2*ones(N*M)# sousdiagonale
DM1[N-1::N]=0.#correction de la sousdiagonale (voisin de gauche n existe pas au bord gauche)
DN=-1/h**2*ones(N*M)
A=spdiags(D0,[0],N*M,N*M)+spdiags(D1,[1],N*M,N*M)+spdiags(DM1,[-1],N*M,N*M)
A=A+spdiags(DN,[N],N*M,N*M)+spdiags(DN,[-N],N*M,N*M)# A est pentadiagonale
#construction d'un second membre trivial
f=ones(N*M)
#
#resolution du systeme lineaire creux
#
t=clock()
u=spsolve(A,f)
t2=clock()
print('resolution systeme done in %5.4f s' % (t2 - t))
#uu=matrix(u,N,M)# permet de reconstruire uu comme un tableau a 2 indices,
#uu(i,j)=u(k), avec k la numerotation de TD
# representation graphique de la solution
uu=reshape(u,(M,N))
uubord=0*ones((M+2,N+2));
uubord[1:M+1,1:N+1]=uu;
#print(A.toarray())
from mpl_toolkits import mplot3d
ax = plt.axes(projection='3d')
X, Y = meshgrid(linspace(0,L,N+2),linspace(0,H,M+2))
ax = plt.axes(projection='3d')
#ax.contour3D(X, Y, uubord, 50, cmap='binary')
ax.plot_surface(X, Y, uubord, rstride=1, cstride=1,
cmap='viridis', edgecolor='none');
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u');
J'ai fais:
- Code: Tout sélectionner
import numpy as np
import matplotlib.pyplot as plt
import pylab
import scipy as sp
from time import clock
from scipy import sparse
from scipy.sparse.linalg import spsolve
from scipy.sparse import spdiags
u=lambda x, y : cos(x) + sin(y)
f=lambda x, y : -cos(x) - sin(y)
L=2.;
H=3.;
N=20;
M=15;
def erreur_Dirchlet2(N, M, u, f, p):
l=L/(N+1);
h=H/(M+1);
x=linspace(l,L-l,N)
y=linspace(h,H-h,M)
xx = []
yy = y*ones((1,M))
a = ones((N,1))
for i in range(len(a)):
xx.append([x[i]*a[i][0]])
ff = f(xx,yy)
F = ff.reshape(N*M)
D0=(2/(l**2)+2/(h**2))*ones(N*M)
D1=-1/l**2*ones(N*M)
D1[N::N]=0.
DM1=-1/l**2*ones(N*M)
DM1[N-1::N]=0.
DN=-1/h**2*ones(N*M)
A=spdiags(D0,[0],N*M,N*M)+spdiags(D1,[1],N*M,N*M)+spdiags(DM1,[-1],N*M,N*M)
A=A+spdiags(DN,[N],N*M,N*M)+spdiags(DN,[-N],N*M,N*M)
U=spsolve(A,F)
err=abs(U-u(xx,yy))
normerr=sum([ei**p for ei in err])
normerr=normerr**(1/p)
normeU=sum([abs(ui)**p for ui in U])
normeU=normeU**(1/p)
return normerr/normeU
J'avance pas du tout, si quelqu'un peut aider, je remercie d'avance.
