Discussion :

[ansarollah]

Bonjour,
Je suis un nouvel utilisateur de matlab.
Merci pour votre aide en tous cas.
je suis en train de faire un calcul numérique d’un EDP. et comme j’ai la solution exacte j’essaie de le comparer avec la solution exacte mais elle me donne toujours un erreur très grand, je sais pas ou est l’erreur, car l’erreur entre la solution exacte et la solution approchée doit être très petit.
le problème que je cherche à résoudre est dans le fichier joint.

Et le code que j’ai fait :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
gd =[4 4;
0 0;
0 0;
1 2;
1 2;
0 0];
ns =[99 99;
49 50 ];
sf='c2-c1';
c=10^10;
[dl,bt]=decsg(gd,sf,ns);
[p,e,t]=initmesh(dl,'Hmax',0.083);
N = size(p,2); % number of nodes
u1=zeros(N,1);
k=inline('1','x','y');
%dir
dir=inline('r^(-1)*cos(teth)','teth','r');
%neumann
neu=inline('-r^(-2)*cos(teth)','teth','r');
result=1;
% solution exacte:
u_ex=inline('r^(-1)*cos(teth)','teth','r');
r=zeros(N,1);
[A,R,b]=assema(p,t,1,0,0);
%condition of dirichlet
for E = 1:size(e,2)
    if (e(6,E)==1)
       % node numbers of boundary edge E
       node = e(1,E);
       2
       nodes = e(1:2,E);
       x = p(1,nodes); y= p(2,nodes);
       [teth,rr]=cart2pol(x(1),y(1));
       A(node,node)=c;
       z=feval(dir,teth,rr);
       b(node)=c*z;
   end %end if
   %condition of Neu
   if (e(6,E)==0)
       nodes = e(1:2,E); % node numbers of boundary edge E
       x = p(1,nodes); y = p(2,nodes);
       ds = sqrt((x(1)-x(2))^2+(y(1)-y(2))^2);
       k_ = feval(k,mean(x),mean(y));
       [teth,rr]=cart2pol(mean(x),mean(y));
       g_=feval(neu,teth,rr);
       r(nodes) = r(nodes) + k_*g_*[1; 1]*ds/2;
   end %end if
end
u=(A+R)\(b+r);
%calcule la solution exacte
for j=1:size(p,2)
   [teth,rr]=cart2pol(p(1,j),p(2,j));
   u_ext(j)=feval(u_ex,teth,rr); end
   %calculer l’erreur
    err=norm(u_ext'-u)

j’attends votre aide.