Hola a todos

Tengo una duda que ha tenido entretenido todo el fin de semana, y
todavÃ***a no sé qué estoy haciendo mal.

Para un programa de elementos finitos, estoy definiendo una serie de
elementos con formulación isoparamétrica, para determinar las
integrales de superficie de cada elemento mediante integración numérica.
Supongamos un cuadrilátero, con funciones de interpolación lineales, en
coordenadas naturales (notación de maxima, pero es muy parecida a la de
maple):

N41-xi)*(1+eta)/4;
N31+xi)*(1+eta)/4;
N21+xi)*(1-eta)/4;
N11-xi)*(1-eta)/4;

N:matrix([N1,N2,N3,N4]);

La matriz de derivadas parciales respecto a xi y eta:

DN:diff(N,xi);
DN:addrow(DN,diff(N,eta));

Supongamos que tenemos un elemento cuadrado de lado 1, con una esquina
en (0,0). Sus coordenadas en forma matricial son:

xy:matrix([0,0],[1,0],[1,1],[0,1]);

la matriz jacobiana de la transformación, el determinante y la inversa:

JN.xy;
detJ:determinant(J);
invJ:invert(J);

Tenemos una matrix B y su transpuesta tal que:

B:invJ.DN;
BT:transpose(B);

Se trata de determinar la integral

ke:integrate(integrate(BT.B*detJ,eta,-1,1),xi,-1,1);

[ 2 1 1 1 ]
[ - - - - - - - ]
[ 3 6 3 6 ]
[ ]
[ 1 2 1 1 ]
[ - - - - - - - ]
[ 6 3 6 3 ]
(%o23) [ ]
[ 1 1 2 1 ]
[ - - - - - - - ]
[ 3 6 3 6 ]
[ ]
[ 1 1 1 2 ]
[ - - - - - - - ]
[ 6 3 6 3 ]
Si lo hago lo mismo mediante integración numérica, se supone que la
integral de área definida antes es igual al sumatorio de del valor de
la función (BT.B*detJ) evaluada en los puntos de Gauss y ponderada por
el peso de cada uno de ellos. El número de puntos depende del grado del
polinomio en cada variable, en este caso 1*1=1 pto. de Gauss, con un
peso de Wi*Wj=2*2=4. El punto en este caso deberÃ***a ser xi=eta=0.

Copio y pego el procedimiento en maxima:

(%i28) pgw;
(%o28) [[0, 0, 4]]
(%i29) ke:zeromatrix(4,4)$

(%i30) declare(i,integer);
(%o30) done
(%i31) K:BT.B*detJ$

(%i32) for i thru length(pgw) do ke:ke+subst([xi = pgw[i][1],eta =
pgw[i][2]],K)*pgw[i][3];
(%o32) done
(%i33) ke;
[ 1 1 ]
[ - 0 - - 0 ]
[ 2 2 ]
[ ]
[ 1 1 ]
[ 0 - 0 - - ]
[ 2 2 ]
(%o33) [ ]
[ 1 1 ]
[ - - 0 - 0 ]
[ 2 2 ]
[ ]
[ 1 1 ]
[ 0 - - 0 - ]
[ 2 2 ]


Que obviamente, está mal. Si empleo cuatro puntos de Gauss, en lugar de
uno;

(%i34) pgwgw4;
1 1 1 1
(%o34) [[- -------, - -------, 1], [- -------, -------, 1],
sqrt(3) sqrt(3) sqrt(3) sqrt(3)
1 1 1 1
[-------, -------, 1], [-------, -
-------, 1]] sqrt(3) sqrt(3) sqrt(3) sqrt(3)
(%i35) ke:zeromatrix(4,4)$

(%i36) for i thru length(pgw) do ke:ke+subst([xi = pgw[i][1],eta =
pgw[i][2]],K)*pgw[i][3]; (%o36) done

(%i38) fullratsimp(%);
[ 2 1 1 1 ]
[ - - - - - - - ]
[ 3 6 3 6 ]
[ ]
[ 1 2 1 1 ]
[ - - - - - - - ]
[ 6 3 6 3 ]
(%o38) [ ]
[ 1 1 2 1 ]
[ - - - - - - - ]
[ 3 6 3 6 ]
[ ]
[ 1 1 1 2 ]
[ - - - - - - - ]
[ 6 3 6 3 ]


que es, como era de esperar, la solución exacta. Hasta donde yo sé, si
las funciones de interpolación son lineales, bastarÃ***a con un punto de
integración; si son cuadráticas, son al menos necesarios 4 (dos en cada
coordenada), si son cúbicas, nueve puntos, y asÃ*** sucesivamente.

La pregunta del millón de euros es ¿me estoy colando con el ordende
integración? Si no es asÃ*** ¿dónde estoy metiendo la pata?


Un saludo a todos y gracias si habéis leÃ***do hasta aquÃ***.
--
Un Saludo

Han Solo
The Rebel Alliance

Emacs is not on every system

So what? [...] Do you tell your administrative people to stick with
notepad.exe? Do you tell your fat kids they can only have the crummy
games that come with their video games or plain dress that comes with
Barbie?