Poligonal de Euler.
Relación con el elemento finito de Lagrange

Se utilizará en esta ocasión el método de la poligonal de Euler para aproximar las  soluciones de un problema general del Cálculo de Variaciones que proviene de la búsqueda de las extremales de un funcional del tipo

		F[
          y] = (∫_a^b) F(x, y(x), y' (x)) dx

satisfaciendo ciertas condiciones de contorno de tipo Dirichlet no necesariamente homogéneas: y(a) = (Subscript[y, a]),  y(b) = (Subscript[y, b]).
Veremos también cómo esta poligonal puede ser construida como combinación lineal de un cierto número de funciones de base, que coinciden precisamente con las del espacio de elementos finitos correspondiente al elemento finito de Lagrange. 
Así pues, cuando se aplique el conocido método de Galerkin en el correspondiente espacio de dimensión finita generado por estas funciones de base para la resolución aproximada de un problema de contorno autoadjunto,

		-(p(x) y' (x))' + q(x) y(
              x) = f(x),  x ∈ [a, b]
y(a) = (Subscript[y, a]),  y(b) = (Subscript[y, b])
		
(q y f son continuas en [a, b], p de clase (C^1) y p positiva). Ambas técnicas resultarán equivalentes, como pasaba al aplicar indistintamente el método de Ritz o el de Galerkin con funciones polinómicas o trigonométricas en [a, b] para un problema autoadjunto.

Clear["Global`*"]

Ejemplo

Resolveremos  el problema consistente en encontrar la curva plana de menor longitud que une los puntos (a, (Subscript[y, a])) y (b, (Subscript[y, b])), cuya solución es el segmento de recta definido por ambos puntos, que tiene por ecuación y(x) = 
        (Subscript[y, a]) + (((Subscript[y, b]) - (Subscript[y, a]))/(b - a)) (x - a), x ∈ [a, b]. La obtenendremos como extremal del funcional longitud, definido por

		L[y]=(∫_a^b)(Sqrt[ 1 + y' (((x))^2)]) dx

Definición del funcional asociado

\(Clear[x, y1, p1, F, FF]
 (*\ 
      denotaremos\ \ y ≡ y1, \ y' ≡ p1, \ 
      \(y'\)' ≡ p2\ y\ así\ sucesivamente\ *) 
\)\(F[x_, y1_, p1_] := \@\(1 + p1\^2\)
\)FF[y_] := 
      (∫_a^b) F[x, y, (∂_x) y] \[DifferentialD]x

Introducción de los extremos del intervalo y las condiciones de contorno de tipo Dirichlet de nuestro problema

\(a = 0.; b = 5.; 
ya = 3.; \ yb = \(-2.\); \)

Construcción de la partición del intervalo

Introducimos el número total de subintervalos, de igual longitud, en los que se subdivide el intervalo [a, b]. Pero en este caso concreto ya que sabemos que la solución debe ser una línea recta. Por tanto,  bastará con tomar el propio intervalo sin subdividir para obtener la solución exacta, como veremos más adelante.

n = 3;

Y posteriormente se definen los nodos de la correspondiente partición uniforme del intervalo.

\(Clear[x, h, k]
\)\(h = \(b - a\)\/n; 

x[i_] := a + i\ h\)

Y se implementan en los nodos extremos las condiciones de contorno de tipo Dirichlet introducidas previamente.

\(Clear[y]
\)\(y[0] = ya; 

y[n] = yb; \)

Cálculo de la función de varias variables que se va a optimizar y de sus puntos críticos

\(Clear[G]
\)\(G = 
      Simplify[∑\+\(i = 1\)\%n
            
∫\_\(x[i - 1]\)\%\(x[i]\)F[x, 
                y[i - 1] + \(\((x - x[i - 1])\)\ \((y[i] - y[i - 1])\)\)\/h, 
                \(y[i] - y[i - 1]\)\/h] \[DifferentialD]x]; 


    ptoscrit1 = 
      Chop[NSolve[Table[∂\_\(y[i]\)G == 0, {i, n - 1}], 
          Table[y[i], {i, n - 1}]]]; \)

yy1 = 
      Table[y[i], {i, 0, n}] /. 
        ptoscrit1[[1]]

Construcción de la poligonal correspondiente

\(yy1 = 
      Table[y[i], {i, 0, n}] /. 
        ptoscrit1[[1]]; 


    g1 = ListPlot[
        Table[{x[i], yy1[[i + 1]]}, {i, 
            0, n}], PlotJoined -> True, 
        DisplayFunction -> Identity]; \)

Valor aproximado del extremo del funcional

G /. ptoscrit1[[1]]

Evaluación aproximada de la integral

Puede resultar muy conveniente usar una regla de cuadratura apropiada para evaluar mucho más fácilmente las integrales de la fórmula anterior. En este caso hemos implementado la conocida regla de los trapecios, obteniendo -como veremos- valores muy aceptables, con un esfuerzo computacional pequeño.

\(Clear[GG]
\)\(GG = 
      ∑\+\(i = 1\)\%n 1\/2\ h\ 
          \((F[x[i - 1], y[i - 1], \(y[i] - y[i - 1]\)\/h] + 
              F[x[i], y[i], \(y[i] - y[i - 1]\)\/h])\); 


    ptoscrit2 = 
      NSolve[Table[∂\_\(y[i]\)GG == 0, {i, n - 1}], 
        Table[y[i], {i, n - 1}]]; \)

yy2 = 
        Table[y[i], {i, 0, n}] /. 
          ptoscrit2[[1]];
g2 = 
        ListPlot[
          Table[{x[i], 
              yy2[[i + 1]]}, {i, 0, n}], 
          PlotJoined -> True, 
          DisplayFunction -> Identity];

N[GG /. ptoscrit2[[1]]]

Solución exacta

Si conocemos la solución exacta del  problema de extremales, calculada resolviendo la  correspondiente ecuación de Euler, podemos compararla con las distintas soluciones aproximadas.

y[x_] := (x - a) ((yb - ya)/(b - a)) + ya

Repasamos antes de continuar, las distintas definiciones que tenemos implementadas relacionadas con el símbolo "y", para comprobar que no habrá incompatibilidades.

Information["y", LongForm -> False]

Definimos también la gráfica correspondiente a la solución exacta, aunque no la dibujemos todavía, ya que la utilizaremos posteriormente para realizar comparaciones con las distintas soluciones aproximadas obtenidas.

gexact = 
        Plot[y[x], {x, a, b}, DisplayFunction -> Identity];

El valor exacto que alcanza el funcional  en esta extremal es el siguiente:

FF[y[x]]

Comparación de los valores obtenidos para la solución exacta y las aproximaciones

\(TableForm[
      Table[{x[i], y[x[i]], 
          yy1[[i + 1]], 
          yy2[[i + 1]]}, {i, 0, 
          n}], 
TableHeadings -> {None, 
{"x[i] 
", "y[x[i]]  
", "y[i]  
", "y[i](trap.)
"}}]\)

Vemos que en este caso particular ambas aproximaciones coinciden y también a su vez  con la solución exacta, que es bien conocida. Hay que hacer notar, sin embargo, que esto es un hecho totalmente excepcional.

Contraste gráfico

Show[g1, g2, gexact, 
        DisplayFunction -> $DisplayFunction];

Comparamos en este gráfico las representaciones de la solución exacta, y la que proporcion el método de la poligonal de Euler, junto con la que se obtiene utilizando la regla de los trapecios para aproximar las integrales correspondientes.

Relación con el elemento finito de Lagrange

Expondremos a continuación la idea subyacente en el conocido como elemento finito de Lagrange, que es el ejemplo más simple que se puede estudiar. En este elemento finito  se emplea sólamente el valor en ambos extremos del intervalo de la función que se quiere aproximar, obteniéndose así la simple interpolación de Lagrange de la misma mediante polinomios de grado 1. La idea consiste en utilizar este elemento finito para construir espacios de funciones aproximantes, según las particiones que tomemos del intervalo considerado, pegando trocitos de segmentos de recta, resultando así  funciones lineales a trozos pero que al menos son continuas. De esta manera volvemos a aproximar la verdadera solución del problema mediante una poligonal, como ocurría con la poligonal de Euler; resultando ambos métodos equivalentes cuando se aplican a problemas de contorno de segundo orden autoadjuntos.
El único problema, como sabemos, es que no podremos aplicar este método para aproximar directamente la ecuación de segundo orden de partida, ya que estas funciones poligonales no tienen la regularidad mínima requerida; no sólo no son de clase 2, sino que ni tan siquiera son de clase 1, ya que sólo lo son a trozos. Pero esto bastará para poder aplicarselo a la correspondiente formulación débil del problema de segundo orden de partida, ya que estas funciones splines lineales están incluidas en el espacio (H^1)(a, b)  de las funciones continuas con derivada continua a trozos y acotada (que es uno de los llamados espacios de Sobolev en dimensión uno que tanto se emplean cuando se habla de la formulación débil de distintos problemas variacionales y de contorno).

Base del elemento finito de Lagrange en un intervalo de referencia y afinidad a cualquier otro intervalo

Así pues, la idea del elemento finito de Lagrange, definido sobre un intervalo  de referencia, por ejemplo [0, 1], consiste en obtener el único polinomio de grado 1, p(x) = (Subscript[c, 0]) + (Subscript[c, 1]) x, que tome ciertos valores dados, (Subscript[α, 1]) y (Subscript[α, 2]), en los extremos de dicho intervalo. Es inmediato comprobar que dicho polinomio se expresa como

		p(
          x) = (Subscript[α, 1]) (Subscript[p, 1])(x) + (Subscript[α, 2]) (Subscript[p, 2])(
              x),
						
donde (Subscript[p, 1])(x) = 1 - x y (Subscript[p, 2])(x) = x (se cumple que(Subscript[p, 1])(0) = 1, (Subscript[p, 1])(1)=0 y (Subscript[p, 2])(0) = 0, (Subscript[p, 2])(1) = 1).
Finalmente, si en vez de trabajar en el intervalo de referencia [0, 1], tuviésemos  un intervalo cualquiera [a, b], bastaría con considerar la afinidad x -> ((x - a)/(b - a)) entre ellos , escribiéndose el polinomio como

		\(TraditionalForm\`p(
          x) = \(α\_1\) \(\(p\_1\)(\(x - a\)\/\(b - a\))\) + \(α\_2\) \(\(p\_2\)(\(x - a\)\/\(b - a\))\)\).

Veamos todo esto, con un ejemplo en concreto. Para ello, empezaremos por definir los que se denominan polinomios de base de dicho elemento finito y que no son más que  (Subscript[p, 1])(x) y (Subscript[p, 2])(x).

\(p\_1[x_] := 1 - x
\)\(p\_2[x_] := x
\)Plot[{(Subscript[p, 1])[x],  (Subscript[p, 2])[x]}, {x, 0, 1}];

Ahora si queremos obtener, por ejemplo la solución obtenida en el problema resuelto mediante el método de la poligonal de Euler en el que calculábamos la recta que pasa por los puntos (a, (Subscript[y, a])) y (b, (Subscript[y, b])), nos bastará emplear el correspondiente operador de interpolación lagrangiana:

ya (Subscript[p, 1])[((x - a)/(b - a))] + yb (Subscript[p, 2])[((x - a)/(b - a))] // 
      Simplify

No es más que la ecuación de la recta que pasa por los puntos (a, (Subscript[y, a])) y (b, (Subscript[y, b])), como comprobamos fácilmente:

ya (Subscript[p, 1])[((x - a)/(b - a))] + yb (Subscript[p, 2])[((x - a)/(b - a))] == 
        (x - a) ((yb - ya)/(b - a)) + ya // Simplify

Discretización por elementos finitos

Sean  (Subscript[S, 1])((Subscript[x, 0]),…,(Subscript[x, n])) el espacio de funciones splines lineales continuas en el intervalo [a, b] y \(\^0\)(Subscript[S, 1])((Subscript[x, 0]),…,(Subscript[x, n])) de aquellas que se anulan  en los puntos de la frontera. Entonces, una base de  \(\^0\)(Subscript[S, 1])((Subscript[x, 0]),…,(Subscript[x, n])) es el conjunto {(Subscript[w, 1]), ..., (Subscript[w, n - 1])} de funciones lineales a trozos tales que (Subscript[w, i])((Subscript[x, j])) = 1 si i = j y (Subscript[w, i])((Subscript[x, j])) = 0 si i ≠ j, que definimos a continuación:

Clear[i, t, w]\(\(w[i_]\)[t_] := 
      Which[x[i - 1] ≤ t ≤ x[i], \ 
        p\_2[\(t - x[i - 1]\)\/\(x[i] - x[i - 1]\)]\ , 
        x[i] ≤ t ≤ x[i + 1], \ 
        p\_1[\(t - x[i]\)\/\(x[i + 1] - x[i]\)], True, 0]
\)\(\(w[0]\)[t_] := 
      Which[x[0] ≤ t ≤ x[1], \ 
        p\_1[\(t - x[0]\)\/\(x[1] - x[0]\)]\ , True, 0]
\)w[n][t_] := 
      Which[x[n - 1] ≤ t ≤ x[n],  
        (Subscript[p, 2])[((t - x[n - 1])/(x[n] - x[n - 1]))] , True, 0]

Las gráficas de estas funciones se muestran a continuación.

\(\(Plot[\ \(w[0]\)[t], {t, a, b}, PlotRange -> {0, 1}, 
        AspectRatio -> Automatic, 
        PlotLabel -> \*""\!\(w\_0\)""\ ];\)\)

\(\(Plot[\ \(w[1]\)[t], {t, a, b}, PlotRange -> {0, 1}, 
        AspectRatio -> Automatic\ , 
        PlotLabel -> \*""\!\(w\_i\)""];\)\)

\(\(Plot[\(w[n]\)[t], {t, a, b}, PlotRange -> {0, 1}, 
        AspectRatio -> Automatic\ , 
        PlotLabel -> \*""\!\(w\_n\)""];\)\)

Construcción del correspondiente spline lineal a partir de la base del espacio de elementos finitos

La idea es que cualquier spline lineal puede obtenerse mediante la correspondiente combinación lineal de estas funciones de base, lineales a trozos.

Clear[c, zm]\(y\&_[t_] := 
      ya\ \(w[0]\)[t] + ∑\+\(i = 1\)\%\(n - 1\)\ y[i]\ \ \(w[i]\)[t] + 
        yb\ \(w[n]\)[t]\)

De manera que al evaluar esta función en cada uno de los nodos de la partición considerada  obtendremos el correspondiente valor asignado. Vamos a probar, con algunos nodos concretos.

\(y\&_[t] /. t -> x[0]\)

\(y\&_[t] /. t -> x[4]\)

\(y\&_[t] /. t -> x[n]\)

Vemos así que el método de la poligonal de Euler no es más que la aplicación del conocido como método de Ritz al funcional que se quiere optimizar, empleando como funciones "coordenadas" precisamente estas funciones lineales a trozos que se toman como funciones de base para los splines lineales de clase (C^0).

Ejercicios

1.-  Aproximar, mediante el método de la poligonal de Euler, con pasos h = 1, 2,  ...  las extremales del funcional

		F[
          y] = (∫_1^2) (x y' (((x))^2)  + 
                (4/x) ((y(x))^2) + (8/x) y(x) ) dx
		
que se quiere optimizar en el dominio de funciones

		D={y∈(ℭ^1)([1, 2]):y(1)=2,y(2) = 3    }
		
y comparar con la  solución exacta (que se obtiene verificando que las funciones del tipo y(x) = (Ax^2) + (B/(x^2)) - 1\(TraditionalForm\`\(\(\ \ \)\(\(verifican\)\(\ \ \)\(la\)\(\ \)\)\)\) ecuación de Euler asociada a este funcional).

2.-  Hallar  también una aproximación numérica de la solución del problema de contorno de segundo orden lineal y no homogéneo

		y'' (x) + y(x) = 2  x,  
      x ∈ [0, 1]
y(0) = 0,  y(1) = 1

mediante el método de la poligonal de Euler, empleando  particiones uniformes con 10 y 20 subintervalos. Comparar los resultados  con la solución exacta.

3.- Resolver aproximadamente el problema de contorno autoadjunto

		y'' (x) - y(x) = sen x,  
      x ∈ [0, π]
y(0) = 0,  y(π) = 0

empleando el método de la poligonal de Euler con 10 y 20 subintervalos. Comparar los resultados  con la solución exacta y con los del ejemplo de la práctica correspondiente al método de los elementos finitos.

4.-  Para encontrar ahora una aproximación de la curva plana de longitud fija  l > (Sqrt[(((b - a))^2) + ((((Subscript[y, b]) - (Subscript[y, a])))^2)])  que une los puntos (a, (Subscript[y, a])) y (b, (Subscript[y, b])) y que minimiza el momento estático dado por el  funcional

		F[
          y] = (∫_0^1) y(x)  (Sqrt[1 + y' (((x))^2)]) 
            dx             
                                                               
satisfaciendo las condiciones de contorno y(a) = (Subscript[y, a]) e y(b) = (Subscript[y, b]), programar el método de la poligonal de Euler  junto con una ligadura de tipo integral, dada por el hecho de que la curva plana tiene una longitud fijada de antemano. Así pues, habrá que minimizar una función de varias variables con una restricción, para lo que se empleará el método de los multiplicadores de Lagrange, para resolver el sistema no lineal resultante. Utilizar para ello una partición  uniforme del intervalo con  2 y 3  subintervalos respectivamente, y  comparar estos resultados con la solución exacta que se obtiene resolviendo el correspondiente problema de tipo "isoperimétrico" (con ligadura de tipo integral), que debe de ser una catenaria.

(Nota.- Para realizar este programa, el lector puede inspirarse en el método de los multiplicadores de Lagrange implementado más adelante para la aproximación numérica de problemas de Sturm-Liouville; ver la práctica correspondiente).