Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.

Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación

Pedro C. Espinoza H.

CAPÍTULO I
MÉTODOS DE ELEMENTOS FINITOS EN LA SOLUCIÓN DE
LAS ECUACIONES DIFERENCIALES PARCIALES ELÍPTICAS

1. ELEMENTOS FINITOS SOBRE UN MALLADO TRIANGULAR

Consideremos la región del plano Ω=[0 , 1]x[0 , 1]. Tomando una partición
uniforme de 6 puntos en [0 , 1] (h=1/5), se genera un mallado con 36 nodos en el
cuadrado Ω. 16 son nodos centrales y el resto está en la frontera (Fig. 1).
Fig. 1 Fig. 2
1 1

0.9 0.9

0.8 0.8

0.7 0.7

0.6 0.6

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Subdividiendo cada cuadrado del mallado en dos triángulos (Fig. 2) se tiene un

mallado triangular de Ω. Consideremos en cada nodo central los 6 triángulos
concurrentes, por ejemplo en el nodo q1  (h, h)  (0.2,0.2) se tienen 6 triángulos
concurrentes (Fig. 3), que lo denotamos EF(1). Enumerado los triángulos
concurrente en el sentido horario: T(1,1) o brevemente (1,1), luego sigue
T(1,2)=(1,2), hasta T(1,6)= (1,6). (Fig. 4)
Fig. 3 Fig. 4
1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
0 0.2 0.4 0.6 0.8 1

En cada uno de estos hexágonos se define un Elemento Finito Lineal (EFL) k que

es una función combinada o ensamblada por porciones de planos, esto es funciones

1
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.

de la forma z  k ( x, y )  ax  by  c . Estas funciones valen 1 en el nodo central y en

los dos nodos restantes se anulan.

Por ejemplo en el triángulo T(1,1) (Fig. 4) cuyos vértices son (h,0), (0,h) y el nodo
central q1  (h, h) , se define el EF 1 ( x, y )  ax  by  c que vale 1 en q1 y en los
otros se anula. Esto conduce al siguiente sistema de ecuaciones:
1 (h,0)  ah  c  0

(1) 1 (0, h)  bh  c  0
 (h, h)  ah  bh  c  1
 1
x y
Cuya solución es a=1/h, b=1/h, c= -1, de este modo 1 ( x, y ) 
  1 en T(1,1).
h h
En general si se tienen 3 puntos no colineales p1=(x1,y1), p2=(x2,y2), p3=(x3,y3)
que representan los vértices de cada uno de los triángulos descritos, lo que se
busca es resolver el siguiente sistema de ecuaciones:
1 ( x1, y1 )  ax1  by1  c  z1

(2) 1 ( x2 , y2 )  ax2  by2  c  z2
 ( x , y )  ax  by  c  z
 1 3 3 3 3 3

La solución de este sistema se hará en general con la ayuda de Matlab.

Programa 1
syms x1 x2 x3 y1 y2 y3 z1 z2 z3 a b c M Ma Mb Mc D
p1=[x1, y1]; p2=[ x2, y2]; p3=[ x3, y3]; z=[z1;z2;z3];
M=[p1, 1; p2, 1 ; p3, 1];
Ma=M; Ma(:,1)=z;
Mb=M; Mb(:,2)=z;
Mc=M; Mc(:,3)=z;
D=det(M);
a=det(Ma)/D; b=det(Mb)/D; c=det(Mc)/D;

Ejemplo 1
En el caso del EF 1 ( x, y ) , cuyo soporte y gráfica se muestran abajo
Fig. 5 Fig. 6
ELEMENTO FINITO EF(1)

0.5

0
1
0.8 1

0.6 0.8
0.6
0.4
0.4
0.2 0.2
0 0

2
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.

Los valores de a, b y c , serán determinados con el Programa 1, dándole valores a las

variables x1, x2, x3, y1, y2, y3, z1, z2, z3, con las coordenadas de los nodos que forman
cada uno de los 6 triángulos que forman su soporte.
Para esto ingresamos previamente:
h=1/5; G1=zeros(6,2);
Tabla1
EN EL TRIANGULO (1,1) EN EL TRIANGULO (1,2) EN EL TRIANGULO (1,3)
x1=h; y1=0; z1=0; x1=0; y1=h; z1=0; x1=0; y1=2*h; z1=0;
x2=0; y2=h; z2=0; x2=0; y2=2*h; z2=0; x2=h; y2=2*h; z2=0;
x3=h; y3=h; z3=1; x3=h; y3=h; z3=1; x3=h; y3=h; z3=1;
% Coeficientes % Coeficientes % Coeficientes
C11=eval([a,b,c]); C12=eval([a,b,c]); C13=eval([a,b,c]);
%C11 = 5 5 -1 %C12 = 5 0 0 %C13 = 0 -5 2
% Gradiente % Gradiente % Gradiente
G1(1,:)=[C11(1), C11(2)]; G1(2,:)=[C12(1), C12(2)]; G1(3,:)=[C13(1), C13(2)];
%G11 = 5 5 %G12 = 5 0 %G13 = 0 -5

Tabla2
EN EL TRIANGULO (1,4) EN EL TRIANGULO (1,5) EN EL TRIANGULO (1,6)
x1=h; y1=2*h; z1=0; x1=2*h; y1=0; z1=0; x1=h; y1=0; z1=0;
x2=2*h; y2=h; z2=0; x2=2*h; y2=h; z2=0; x2=2*h; y2=0; z2=0;
x3=h; y3=h; z3=1; x3=h; y3=h; z3=1; x3=h; y3=h; z3=1;
% Coeficientes % Coeficientes % Coeficientes
C14=eval([a,b,c]); C15=eval([a,b,c]); C16=eval([a,b,c]);
%C14 = -5 -5 3 %C15 = -5 0 2 %C16 = 0 5 0
% Gradiente % Gradiente % Gradiente
G1(4,:)=[C14(1), C14(2)]; G1(5,:)=[C15(1), C15(2)]; G1(6,:)=[C16(1), C16(2)];
%G14 = -5 -5 %G15 = -5 0 %G16 = 0 5

La definición precisa de 1 en el triángulo T (1,1) es:

1 ( x, y )  5 x  5 y  1; si ( x, y )  T (1,1) ,
y la descripción del triángulo: T (1,1)   ( x, y ) / 0  x  h; h  x  y  h 

2. UNA ECUACIÓN ELIPTICA DE TIPO DIRICHLET

  u  u  f en 
(1) 
 u  0 en 

Forma débil del modelo (1)

(1)  u

 v  uv dxdy   fvdxdy para todo v  H 01 ()


donde Ω=[0 , 1]x[0 , 1].

a) La forma débil del modelo discreto

3
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.
16
Supongamos que la solución aproximada de (1) es u h   x i i
entonces se puede
i 1

plantear la forma débil del modelo (1) como sigue

(2)  uh   k  uh k dxdy   fk dxdy para cada k 1, 2,....,16
 

b) La matriz de entrada F
Es el lado derecho de (2) con Fk   f dxdy , donde S (k )  Sop(
k k )
S (k )

Estos términos se pueden aproximar tomando qk el nodo central del Sop( k ) de y

evaluando en este punto la función f , la integral se aproxima por
Fk  f (qk )   dxdy
k
S (k )
6 6
Sea c k   k dxdy , pero S (k )   T(k, i) entonces ck  
S (k ) i 1 i 1
  dxdy
k
T(k, i)
luego Fk  f (qk )ck con k  1, 2...,16 y F  [ Fk ] es la matriz columna.

c) De la forma débil al sistema de ecuaciones

16
Reemplazando u h   x i i
en (2) se pasa a un sistema de 16 ecuaciones lineales
i 1

16 16
(3)  xi  i  k dxdy   xi    dxdy  f (q )   dxdy
i k k k ,  k 1, 2,....,16
i 1 S (k ) i 1 S(k ) S (k )

Con las notaciones:

ai , k    i   k dxdy ; bi , k     dxdy ;
i k Fk  f (qk )   k dxdy
S (k ) S (k ) S(k )

El sistema (3) toma la forma:

a1,1 x1  ...  a1,16 x16  b1,1 x1  ...  b1,16 x16  F1

a2,1 x1  ...  a2,16 x16  b2,1 x1  ...  b1,16 x16  F2
(4) 
..............................................................
a16,1 x1  ..  a16,16 x16  b16 ,1 x1  ..  b16,16 x16  F16

que matricialmente es Ax  Bx  F o ( A  B) x  F

donde A  [ai , k ], B  [bi , k ] y F  [ Fk ]

3. CÁLCULO DE LAS COMPONENTES DE LAS MATRICES A, B y F

4
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.

3.1 Cálculo de las componentes de la Matriz de Rigidez A

a1,1   1  1dxdy
SOP
6
Sop (i )   T (i, k )
k 1
6
a1,1     1 1dxdy
k 1 T (1,k )

1  G1,k en T (1, k ) , como 1 es lineal a trozos 1  G1, k es un vector constante.
6
a1, 1   (G1, k  G1,k )  dxdy
k 1 T (1, k )
6 6
2
 dxdy  area(T (1, k ))  1 h  50  a1,1   (G1,k  G1, k )area(T (1, k ))  501  (G1, k  G1,k )
T (1,k )
2 k 1 k 1

En las tablas 1 y 2 se tienen los valores de los gradientes G1, k . Entonces como en

Matlab G1, k  G1,k = dot(G1(k,:), G1(k,:)) se tiene en valor de la suma a1,1 , corriendo

el siguiente programa:
Suma=0;
for i=1:6
Suma=Suma+ dot(G1(i,:), G1(i,:));
end
Sale Suma=200. Finalmente a1,1 =200/50=4.

En forma similar se calculan las otras componentes de la matriz A.

3.2 Cálculo de las componentes de la Matriz de Masa B
Para calcular la primera componente de la matriz B  [bi , j ] esto es b1,1     dxdy
1 1
SOP
6
se tiene que integrar en cada triángulo del soporte de 11 = SOP   T (i, k ) .
k 1
2
a) La componente b1,1     dxdy   ( )
1 1 1 dxdy
SOP SOP

La integral de (1 ) 2 en el triángulo T (1,1)

La descripción de T (1,1) como región del plano XY es
T (1,1) : 0  x  h; h  x  y  h
La función 1 definida en el triángulo T (1,1) es
z  1 ( x, y )  5 x  5 y  1
h h
2
Luego J 11   11dxdy    (5x  5 y  1) dydx
T (1,1) 0 hx

Integrando con la ayuda de Matlab:

syms x y
J11=int(int((5*x+5*y-1)^2,'y',h-x,h),'x',0,h)

5
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.

J11=1/300

b) La integral de (1 ) 2 en el triángulo T (1,2)

T (1,2) : 0  x  h; h  y  2h  x ; z  1 ( x, y )  5 x
h 2h x
2
J 12   11dxdy  
T (1, 2 ) 0
 (5 x )
h
dydx

J12=int(int((5*x)^2,'y',h,2*h-x),'x',0,h)
J12=1/300
En forma similar se calculan las integrales de (1 ) 2 en los triángulos restantes,
obteniendo finalmente: b1,1     dxdy =6/300=1/50
1 1
SOP

3.3 Cálculo de las componentes de la Matriz de Entrada F

Esto es inmediato tomando en cuenta que Fk  f (qk )   k dxdy y ck    dxdyk
S(k ) S (k )

y aprovechando las integrales calculadas en 3.2

c11=int(int((5*x+5*y-1),'y',h-x,h),'x',0,h)
c11=1/150
c12=int(int(5*x,'y',h,2*h-x),'x',0,h)
c11=1/150
En forma similar se calculan las integrales de 1 en los triángulos restantes. Un
argumento que se puede usar es que las integrales sobre cada triángulo son los
volúmenes de prismas, donde todas son simplemente traslaciones y/o rotaciones
del primer prisma cuya base es el triángulo T(1,1)
Finalmente: c11   1dxdy =6/150=1/25=h^2
SOP

4. IMPLEMENTACIÓN DEL MÉTODO DE LOS ELEMENTOS FINITOS

Ejemplo 1
Consideremos un problema de Dirichlet donde se conoce su solución analítica o
exacta ue  ( x  x 2 )( y  y 2 )

 u  u  u  ( x  x 2 )( y  y 2 )  2( x  x 2  y  y 2 ) en 
(2) 
 u  0 en 
4.1 Implementación de las matrices y solución del sistema (4)
h=1/5;
x=[h, 2*h, 3*h, 4*h, h, 2*h, 3*h, 4*h, h, 2*h, 3*h, 4*h, h, 2*h, 3*h, 4*h];
y=[ h, h, h, h, 2*h, 2*h, 2*h, 2*h, 3*h, 3*h, 3*h, 3*h, 4*h, 4*h, 4*h, 4*h];
%Matriz de Entrada
F0=(x-x.^2).* (y-y.^2)+2*( x-x.^2+y-y.^2);
C= h^2*ones(16,1);
F=F0'.*C;

6
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.

% Matriz de Rigidez
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];
%Matriz de Maza
B0=[6 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
1 6 1 0 1 1 0 0 0 0 0 0 0 0 0 0;
0 1 6 1 0 1 1 0 0 0 0 0 0 0 0 0;
0 0 1 6 0 0 1 1 0 0 0 0 0 0 0 0;
1 1 0 0 6 1 0 0 1 0 0 0 0 0 0 0;
0 1 1 0 1 6 1 0 1 1 0 0 0 0 0 0;
0 0 1 1 0 1 6 1 0 1 1 0 0 0 0 0;
0 0 0 1 0 0 1 6 0 0 1 1 0 0 0 0;
0 0 0 0 1 1 0 0 6 1 0 0 1 0 0 0;
0 0 0 0 0 1 1 0 1 6 1 0 1 1 0 0;
0 0 0 0 0 0 1 1 0 1 6 1 0 0 1 1;
0 0 0 0 0 0 0 1 0 0 1 6 0 0 1 1;
0 0 0 0 0 0 0 0 1 1 0 0 6 1 0 0;
0 0 0 0 0 0 0 0 0 1 1 0 1 6 1 0;
0 0 0 0 0 0 0 0 0 0 1 1 0 1 6 1;
0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 6] ;
B=(h^2/12)*B0;
%Solución del sistema de ecuaciones (4): la solución aproximada uh
uh=inv(A+B)*F;
De otro lado la solución analítica es ue  ( x  x 2 )( y  y 2 ) , entonces haciendo
ue=(x-x.^2).*(y-y.^2); se tiene los valores de ue en cada nodo de la partición.
La siguiente tabla muestra los nodos y los correspondientes valores de ue y uh.
format short
TABLA=[x' y' ue' uh]
x y u uh
0.2000 0.8000 0.0256 0.0258
0.4000 0.8000 0.0384 0.0386
0.6000 0.8000 0.0384 0.0386

7
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.

0.8000 0.8000 0.0256 0.0257

0.2000 0.6000 0.0384 0.0386
0.4000 0.6000 0.0576 0.0580
0.6000 0.6000 0.0576 0.0579
0.8000 0.6000 0.0384 0.0386
0.2000 0.4000 0.0384 0.0386
0.4000 0.4000 0.0576 0.0579
0.6000 0.4000 0.0576 0.0580
0.8000 0.4000 0.0384 0.0387
0.2000 0.2000 0.0256 0.0257
0.4000 0.2000 0.0384 0.0386
0.6000 0.2000 0.0384 0.0387
0.8000 0.2000 0.0256 0.0258

4.2 Gráfica de la solución aproximada uh

a) Nodos de la partición

h=1/5;
xx=[0,h,2*h,3*h,4*h,1,0,h,2*h,3*h,4*h,1,0,h,2*h,3*h,4*h,1,0,h,2*h,3*h,4*h,1,0,h,2*h,3*h,4*h,1,0,h,2*h,3*h,4*h,1];
yy=[0,0,0,0,0,0,h,h,h,h,h,h,2*h,2*h,2*h,2*h,2*h,2*h, 3*h, 3*h,3*h, 3*h, 3*h, 3*h,4*h,4*h,4*h, 4*h,4*h, 4*h, 1, 1, 1,
1, 1,1];

uh=uh';
%b) Valores de la solución aproximada uh
N1=1:1:6; N2=7:1:12; N3=13:1:18; N4=19:1:24; N5=25:1:30; N6=31:1:36;
xx1=xx(N1) ;xx2= xx(N2) ;xx3= xx(N3) ; xx4= xx(N4); xx5= xx(N5);xx6= xx(N6);
yy1=yy(N1);[0,0,0,0,0,0];yy2=yy(N2);yy3=yy(N3); yy4=yy(N4); yy5=yy(N5);yy6=yy(N6);
%c) Puntos para las trayectorias transversales
bb=xx; aa=yy;
aa1=aa(N1);aa2= aa(N2);aa3= aa(N3); aa4= aa(N4); aa5= aa(N5);aa6= aa(N6);
bb1=bb(N1);bb2= bb(N2);bb3= bb(N3); bb4= bb(N4); bb5= bb(N5);bb6= bb(N6);
J1=1:1:4; J2=5:1:8; J3=9:1:12; J4=13:1:16;
%d) Extensión de uh a los nodos de la frontera
uh1=[0,0*uh(J1),0];
uh2=[0,uh(J1),0];
uh3=[0,uh(J2),0];
uh4=[0,uh(J3),0];
uh5=[0,uh(J4),0];
uh6=[0,0*uh(J4),0];
%e)Ploteando los puntos de la gráfica de uh
plot3(xx1,yy1, uh1, 'or', 'MarkerSize',8,'MarkerFaceColor',[0 1 0]);
hold on
grid
plot3(xx2,yy2, uh2, 'or', 'MarkerSize',8,'MarkerFaceColor',[0 1 0]);
plot3(xx3,yy3, uh3, 'or', 'MarkerSize',8,'MarkerFaceColor',[0 1 0]);
plot3(xx4,yy4, uh4, 'or', 'MarkerSize',8,'MarkerFaceColor',[0 1 0]);
plot3(xx5,yy5, uh5, 'or', 'MarkerSize',8,'MarkerFaceColor',[0 1 0]);
plot3(xx6,yy6, uh6, 'or', 'MarkerSize',8,'MarkerFaceColor',[0 1 0]);
plot3(xx1,yy1,uh1,'r', 'LineWidth',2);

8
Universidad Nacional Mayor de San Marcos. Facultad de Ciencias Matemática.
Maestría en Matemática Aplicada. Curso: Modelaje Gráfico y Simulación
Pedro C. Espinoza H.

plot3(xx2,yy2, uh2,'r', 'LineWidth',2);

plot3(xx3,yy3, uh3,'r', 'LineWidth',2);
plot3(xx4,yy4, uh4,'r', 'LineWidth',2);
plot3(xx5,yy5, uh5,'r', 'LineWidth',2);
plot3(xx6,yy6, uh6,'r', 'LineWidth',2);
%f)Trayectorias transversales
plot3(aa1,bb1,uh1,'b', 'LineWidth',2);
plot3(aa2,bb2, uh2,'b', 'LineWidth',2);
plot3(aa3,bb3, uh3,'b', 'LineWidth',2);
plot3(aa4,bb4, uh4,'b', 'LineWidth',2);
plot3(aa5,bb5, uh5,'b', 'LineWidth',2);
plot3(aa6,bb6, uh6,'b', 'LineWidth',2);
legend('= uh');
title('ELEMENTOS FINITOS – UPG– FCM – UNMSM–2017–1');

ELEMENTOS FINITOS – UPG– FCM – UNMSM–2016–2

= uh

0.06

0.05

0.04

0.03

0.02

0.01

0
1
0.8 1
0.6 0.8
0.4 0.6
0.4
0.2 0.2
0 0

4.3 Gráfica de la solución exacta con Matlab

h=1/5;
ax=0:h:1;by=ax;
[x,y]=meshgrid(ax,by);
z=(x-x.^2).*(y-y.^2);
surf(x,y,z);
colormap(hsv)

0.06

0.04

0.02

0
1
1
0.5
0.5

0 0

a) Valores de la solución analítica u en los puntos del mallado

0 0 0 0 0 0

9
a) Valores de la solución analítica u en los puntos del mallado

0 0 0 0 0 0

0 0.0256 0.0384 0.0384 0.0256 0

0 0.0384 0.0576 0.0576 0.0384 0
0 0.0384 0.0576 0.0576 0.0384 0
0 0.0256 0.0384 0.0384 0.0256 0
0 0 0 0 0 0

b) Valores de la partición del cuadrado = [0,1]x[0,1].

x=
0 0.2000 0.4000 0.6000 0.8000 1.0000

0 0.2000 0.4000 0.6000 0.8000 1.0000

0 0.2000 0.4000 0.6000 0.8000 1.0000

0 0.2000 0.4000 0.6000 0.8000 1.0000

0 0.2000 0.4000 0.6000 0.8000 1.0000

0 0.2000 0.4000 0.6000 0.8000 1.0000

y=
0 0 0 0 0 0
0.2000 0.2000 0.2000 0.2000 0.2000 0.2000

0.4000 0.4000 0.4000 0.4000 0.4000 0.4000

0.6000 0.6000 0.6000 0.6000 0.6000 0.6000

0.8000 0.8000 0.8000 0. 8000 0. 8000 0. 8000

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000