Viga de un solo vano: Script en Máxima

En la viga de un solo vano empezamos a tratar Máxima como una calculadora. Cada vez que nos daba una respuesta, le pedíamos que trabajase con ella. Ahora, lo que vamos a intentar es modificar lo que hicimos para que sirva para casi cualquier cosa. Es decir, hacer un script más genérico...

Usaremos

Lo que ya aprendimos en la viga de un solo vano: Integrate, subst, solve...

• first(A): Para acceder al primer elemento de la lista A (igual funcionan second(A), third(A)...)
  
• part(A,i): Para acceder al elemento en posición i de la lista A
• rhs: accede a lo que hay a la derecha del igual en una ecuación
• map(función, lista): Aplica una determinada función a cada elemento de la lista
  

Para mayor claridad (y para tener que subir menos imágenes) vamos a enviar las órdenes en grupos.

Integrando

Esta parte es casi igual que en el artículo anterior, pero asignamos e integramos en una sola orden:

v:integrate(q,x)+v0;
m:integrate(v,x)+m0;
c:m/EI;
g:integrate(c,x)+g0;
f:integrate(g,x)+f0;

Podríamos modificar q, y en vez de ser un valor constante, podría ser cualquier función continua.

Ecuaciones de borde

Aquí si hay modificaciones y aclaraciones:
La respuesta de solve() es una lista con todas las soluciones posibles. Es decir, si el sistema no fuera lineal habría más de una solución posible. Cada solución sería un elemento de la lista.
Y si pedimos que nos resuelva más de una incógnita (como es el caso) entonces cada elemento de la lista de soluciones sería, a su vez, una lista de incógnitas resueltas. Solve nos da una lista de listas...
Sabemos que la lista de soluciones sólo tendrá una solución (es un sistema lineal), así que sólo nos interesa el primer elemento: Asignamos el primer elemento (first) a la variable soluciones:

/*Viga biempotrada: los apoyos ni bajan ni giran: 4 ecuaciones, 4 incógnitas */
ec1:subst(x=0,f)=0;
ec2:subst(x=l,f)=0;
ec3:subst(x=0,g)=0;
ec4:subst(x=l,g)=0;
solve([ec1,ec2,ec3,ec4], [v0,m0,g0,f0]);
soluciones:first(%);

Definición de las funciones de la viga

Una vez asignadas los distintos valores de las distintas incógnitas a la lista "soluciones", las metemos en las funciones con la orden subst. Comento la primera línea, el resto son similares:
Definimos la funcion v sustituyendo (subst) v0 por el valor que hemos obtenido:
De la lista "soluciones" necesitamos el primer valor y lo cogemos con part(soluciones,1).
Todos los valores de la lista "soluciones" son igualdades. Lo que queremos es lo que está a la derecha del signo "=" y lo obtenemos con la función rhs().
En el resto de sustituciones, hay que sustituir más de un valor, por lo que metemos la lista de valores a sustituir entre corchetes (por que es una lista)

/*Conocidas las constantes de integracion, escribo las ecuaciones completas */
"Cortante";
v:subst(v0=rhs(part(soluciones,1)),v);
"Momento";
m:subst([v0=rhs(part(soluciones,1)), m0=rhs(part(soluciones,2))],m);
"Giro";
g:subst([v0=rhs(part(soluciones,1)), m0=rhs(part(soluciones,2)), g0=rhs(part(soluciones,3))],g);
"Flecha";
f:subst([v0=rhs(part(soluciones,1)), m0=rhs(part(soluciones,2)), g0=rhs(part(soluciones,3)), f0=rhs(part(soluciones,4))],f);

Análisis de las funciones

Una vez definidas las funciones, podemos buscar los ceros, los valores en los extremos y los valores máximos. Como cada función es la integral de la anterior, los ceros de unas nos indicarán los máximos de otras.
Buscamos las x que hacen que la función sea cero con solve. y las guardamos en una lista: xv0 sería la lista de valores de x que hacen que el cortante sea cero.
Ahí donde el cortante es cero, el momento será máximo.
Y para localizar los valores máximos usamos la función map, que aplicará la función sobre todos los elementos de la lista:
La función a aplicar lambda la definimos con la lista de parámetros que necesita ([i]) y como se aplica a la finción.
Al estar anidadas es un poco más confuso de ver, pero si nos fijamos donde se abren y se cierran los paréntesis, es muy aclaratorio.
Definimos listas con nombre .ext para los valores en los extremos, y listas con nombre .max para valores máximos.

/*busco los ceros y miro los máximos y los valores en los extremos*/
"Cortante extremos";
vext:map(lambda([i], subst(x=i,v)),[0,l]);
"Cortante=0 en"; xv0:solve(v=0,x);
"Momento max, Momento extremos";
mmax:map(lambda([i], subst(x=rhs(i),m)),xv0);
mext:map(lambda([i], subst(x=i,m)),[0,l]);
"Momento=0 en"; xm0:solve(m=0,x);
"Giro max, giro extremo";
gmax:map(lambda([i], subst(x=rhs(i),g)),xm0);
gext:map(lambda([i], subst(x=i,g)),[0,l]);
"Giro=0 en"; xg0:solve(g=0,x);
"Flecha max, flecha en extremos";
fmax:map(lambda([i], subst(x=rhs(i),f)),xg0);
fext:map(lambda([i], subst(x=i,f)),[0,l]);
"Flecha=0 en"; xf0:solve(f=0,x);

(Pulsa en la imagen, porque si no, te quedarás miope)

Varios finales alternativos

¿Para que vale el script?
Para adaptarlo a los diferentes casos que se nos ocurran:

1. Podemos modificar la ley de cargas, y estudiar casos en los que la carga no es constante.
2. Podemos modificar las condiciones de borde, y estudiar sin mucho trabajo otros casos típicos o no tan típicos
3. Podemos convertir algunas condiciones de borde en compatibilidad de deformaciones y equilibrio.

Lo más difícil ya está hecho. Ahora sólo se trata de jugar con ello... (otro día)