Introducción al M.E.F. en una variable

Resolución aproximada de problemas de segundo orden

Aplicaremos a continuación el conocido método de Galerkin en el  espacio de elementos finitos que corresponde al elemento finito de Lagrange, 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 y p es de clase (C^1)([a, b]), con p(x) > 0). Para ello emplearemos la correspondiente formulación débil del mismo. Como ya sabemos, una vez fijada una partición concreta en el intervalo [a, b], dicho espacio de funciones aproximantes coincide con el de los splines lineales relativos a la partición, resultando equivalente este método con el de la poligonal de Euler aplicado al funcional correspondiente.

Ejemplo

Resolveremos aproximadamente el siguiente problema de contorno autoadjunto,

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

Clear["Global`*"]

Introducción de los datos del problema

Funciones que definen el problema autoadjunto de segundo orden

\(p[x_] := 1
\)\(q[x_] := 1
\)f[x_] := -Sin[x]

Dominio y condiciones de contorno

\(a = 0; \ \ y\_a = 0; 

b = π; \ y\_b = 0; \)

Cambio de variable dependiente para homogeneizar las condiciones de contorno de tipo Dirichlet

\(α[x_] := y\_a + \((x - a)\) \(y\_b - y\_a\)\/\(b - a\)
\)\(f\&^[x_] := 
      f[x] - q[x]\ α[x] + \(p'\)[x]\ \ \(α'\)[x]\)

Formulación débil del problema

Sea w una función test, en el espacio adecuado,  verificando las condiciones de contorno homogéneas w(a) = w(b) = 0. Entonces se tiene que

		(∫_a^b)(-(p(x) y' (x))'  +  
                    q(x) y(x))  w(x)dx=(∫_a^b) f(x) w(x) dx
		
y, aplicando la fórmula de integración por partes en el primer término de la integral de la izquierda, se obtiene que

		𝔄( y, w) = (Subscript[𝔉, f])( w ),
siendo

		𝔄( y, w ) = (∫_a^b)p(x) y' (x) ϕ' (x) + q(x) y(x)  w(x)  dx
(Subscript[𝔉, f])( w) = (∫_a^b) f(x) w(x) dx

El problema variacional o débil asociado al problema de contorno inicial consiste en hallar  una función y(x) continua y con derivada primera continua a trozos y acotada, satisfaciendo las condiciones de contorno de tipo Dirichlet dadas y tal que

		(Vd)       𝔄( y,  w ) = (Subscript[𝔉, f])( w ) 
                           
para toda función test w en  el espacio de Sobolev  (Subscript[H, 0]^1) (a, b) que en el caso unidimensional coincide precisamente con el de las funciones continuas con derivada primera continua a trozos y acotada que se anulan en los puntos de la frontera.

Construcción de la partición del intervalo, con nodos equidistantes

Por simplicidad, podemos considerar  una partición del dominio con n + 1 nodos equidistantes, es decir,  a =(Subscript[x, 0]) < (Subscript[x, 1]) < ... < (Subscript[x, n]) = b, donde (Subscript[x, i]) = a + i b,  para i = 0, 1, ..., n,  siendo h=((b - a)/n), de manera que se subdivide el intervalo de partida en n subintervalos de igual longitud.

\(Clear[n, h]
\)\(n = 7; 

h = \(b - a\)\/n; \)

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

Base del correspondiente espacio de funciones aproximantes a partir del elemento finito de refererncia

Ahora, dada una partición a =(Subscript[x, 0]) < (Subscript[x, 1]) < ... < (Subscript[x, n]) = b del intervalo [a, b], de lo que se trata es de construir un spline lineal que tome ciertos valores en los nodos de dicha partición y de manera que aproxime a la solución de cierto problema de contorno. La  base de este espacio  (^0)(Subscript[S, 1])((Subscript[x, 0]), ..., 
              (Subscript[x, n]))  estará formada por las funciones lineales a trozos (Subscript[w, i]), tales que (Subscript[w, i])((Subscript[x, j])) = (Subscript[δ, ij]), i, j = 0, 1, ..., n.

Clear[i, t, w]
\(p\_1[x_] := 1 - x
\)\(p\_2[x_] := x
\)\(\(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]

Matrices de rigidez y de carga

Sea y la solución del problema (Vd), entonces,

		\(TraditionalForm
      \`\(y\&^\)(x) = 
        \(y\_a\) \(\(w\_0\)(x)\) + 
          ∑\+\(i = 1\)\%\(n - 1\)\(c\_i\) \(\(w\_i\)(x)\) + 
          \(y\_b\) \(\(w\_n\)(x)\)\)

y, por consiguiente, el problema (Vd) es equivalente a la resolución del sistema de ecuaciones lineales

		(L)      (Subscript[A, rig]) c = (Subscript[d, car])

donde (Subscript[A, rig])  se denomina matriz de rigidez, (Subscript[d, car])  vector de cargas, y  (Subscript[c, i]) = (Subscript[y, i]), para i = 1, 2, ..., n - 1.
Los valores que toman estas  expresiones  son los siguientes:

		(Subscript[A, rig]) =(Subscript[(𝔄( (Subscript[w, i]),(Subscript[w, j]) )), 1 ≤ i, 
            j ≤ n - 1])
(Subscript[d, car]) = (Subscript[(((Subscript[d, j]))), 1 ≤ j ≤ n - 1])

donde (Subscript[d, j])=(Subscript[𝔉, f])( (Subscript[w, j])) -(Subscript[y, a]) 𝔄( (Subscript[w, 0]),(Subscript[w, j]) )-(Subscript[y, b]) 𝔄( (Subscript[w, n]),(Subscript[w, j]) ).

No obstante lo mejor es homegeneizar previamente las correspondientes condiciones de contorno de tipo Dirichlet, mediante el cambio de variable dependiente dado por y(x) = z(x) + α(x) con \(TraditionalForm
      \`α(x) = y\_a + \(\(y\_b - y\_a\)\/\(b - a\)\) \((x - a)\)\), y, de esta manera, aproximaríamos mediante \(TraditionalForm
      \`\(z\&^\)(x) = 
        ∑\+\(i = 1\)\%\(n - 1\)\(c\_i\) \(\(w\_i\)(x)\)\),  para después hallar 

		\(TraditionalForm\`\(y\&^\)(x) = \(z\&^\)(x) + α(x)\).

Este cambio de variable dependiente sólo afectaría al segundo término de la ecuación de segundo orden, quedando \(TraditionalForm
      \`\(f\&^\)(x) = 
        f(x) - \(q(x)\) \(α(x)\) + p' \((x)\) α' \((x)\)\) en lugar de f(x).

Cálculo de la matriz de rigidez

\(A = Table[0, {i, 1, n - 1}, {j, 1, n - 1}]; 


    For[i = 1, i ≤ n - 1, \(i++\), 
      A[\([i, i]\)] = 
        ∫\_\(x[i - 1]\)\%\(x[i]\)\((
                p[x]\ \((
                        \(p\_2'\)[\(x - x[i - 1]\)\/\(x[i] - x[i - 1]\)]
                            \/\(x[i] - x[i - 1]\))\)\^2 + 
                  q[x] \((p\_2[\(x - x[i - 1]\)\/\(x[i] - x[i - 1]\)])\)\^2)
                \) \[DifferentialD]x + 
          
∫\_\(x[i]\)\%\(x[i + 1]\)\((
                p[x]\ \((
                        \(p\_1'\)[\(x - x[i]\)\/\(x[i + 1] - x[i]\)]
                            \/\(x[i + 1] - x[i]\))\)\^2 + 
                  q[x] \((p\_1[\(x - x[i]\)\/\(x[i + 1] - x[i]\)])\)\^2)\) 
              \[DifferentialD]x
	]; 



    For[i = 1, i ≤ n - 2, \(i++\), 
      A[\([i, i + 1]\)] = 
        
∫\_\(x[i]\)\%\(x[i + 1]\)\((
              p[x]\ \((
                    \(p\_1'\)[\(x - x[i]\)\/\(x[i + 1] - x[i]\)]
                        \/\(x[i + 1] - x[i]\))\) 
                  \((\(p\_2'\)[\(x - x[i]\)\/\(x[i + 1] - x[i]\)]
                        \/\(x[i + 1] - x[i]\))\)
 + 
                q[x] \((p\_1[\(x - x[i]\)\/\(x[i + 1] - x[i]\)])\) 
                  \((p\_2[\(x - x[i]\)\/\(x[i + 1] - x[i]\)])\))\) 
            \[DifferentialD]x; 
		
      A[\([i + 1, i]\)] = A[\([i, i + 1]\)]
	]\)\(

\(A\_rig = A; \)\)

Simplify[(Subscript[A, rig])] // MatrixForm

Obsérvese que sólo se han calculado los elementos  a priori no nulos de la matriz de rigidez, teniendo en cuenta la expresión de cada uno de los tramos de las funciones de base. De hecho cuando se aplique este método para un número elevado de subintervalos convendrá almacenar tan sólo los valores de la diagonal principal y la superior, ya que por la simetría del problema la diagonal inferior coincide con aquella, y emplear algún algoritmo especial de resolución de sistemas tridiagonales y simétricas, que además será diagonalmente dominante y/o definida positiva. Nótese también que, cuando la ecuación diferencial sea de coeficientes constantes, entonces también van a coincidir todos los valores de la diagonal principal y todos los de las dos diagonales adyacentes, con lo cual nos podríamos ahorrar todos estos cálculos y efectuarlos sólo una vez.

Cálculo del vector de carga

\(d = Table[0, {i, 1, \ n - 1}]; 


    For[i = 1, i ≤ n - 1, \(i++\), 
      
d[[i]] = 
        ∫\_\(x[i - 1]\)\%\(x[i]\)f\&^[x]\ \ 
                p\_2[\(x - x[i - 1]\)\/\(x[i] - x[i - 1]\)] 
                \[DifferentialD]x + 
            
∫\_\(x[i]\)\%\(x[i + 1]\)f\&^[x]\ \ 
                p\_1[\(x - x[i]\)\/\(x[i + 1] - x[i]\)] \[DifferentialD]x // N
          

]
\)\(d\_car = d; 

MatrixForm[d\_car]\)

Obsérvese que, cuando se quiera agilizar el cálculo efectivo de todas estas integrales, puede ser conveniente el empleo de la sentencia  NIntegrate, que proporcionaría una aproximación numérica de las mismas, lo cual puede ser más que suficiente en muchos casos, evitando así la posible dificultad que pudiera tener el cálculo simbólico de las mismas,  según la expresión del término independiente de la ecuación diferencial.

Resolución del sistema lineal resultante

c = LinearSolve[(Subscript[A, rig]),  (Subscript[d, car])]

Construcción de la solución aproximada y representación gráfica de la misma

\(y\&_[t_] := 
      ∑\+\(i = 1\)\%\(n - 1\)\ c[\([i]\)]\ \ \(w[i]\)[t] + α[t]

    \)\(\(gaprox = Plot[y\&_[t], {t, a, b}]; \)\)

¿Solución exacta?

Si conocemos cuál es la solución exacta del  problema de extremales, calculada resolviendo la  correspondiente ecuación de Euler. Podemos ahora compararla con las distintas soluciones aproximadas.

Clear[y]
DSolve[{y''[x]-y[x]==Sin[x],y[a]==0,y[b]==0},y[x],x]

y[x_] = y[x] /. Chop[%][[1]]

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

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

\(TableForm[Table[{x[i], y[x[i]] // N, y\&_[x[i]]}, {i, 0, n}], 
      
TableHeadings 
        -> {
          None, {"x[i] 
", "y[x[i]]  
", "yy[x[i]]  
"}}]\)

Contraste gráfico

Show[ gaprox, gexact, DisplayFunction -> $DisplayFunction];

Comparamos en este gráfico las representaciones de la solución exacta, junto con la que se obtiene de forma aproximada utilizando el método de los elementos finitos.

Resolución aproximada de problemas de cuarto orden

Se trata de la resolución del correspondiente problema unidimensional de cuarto orden mediante el  método de los elementos finitos.

		(p(x) y'' (x))'' - (q(x) 
                y' (x))' + r(x) y(x) = f(x),  
      x ∈ [a, b]
y(a) = (Subscript[y, a]),  
      y(b) = (Subscript[y, b]);     y' (a) = ((Subscript[y, a])^1),  y' (b) = ((Subscript[y, b])^1)

(r y f son continuas y p y q son de clase (C^2)([a, b]) y (C^1)([a, b]) respectivamente). Para ello emplearemos la correspondiente formulación débil del mismo. 
 Usaremos el conocido elemento finito de Hermite  definido en el intervalo de referencia  I≡[0,1] mediante los datos de la función y la derivada primera en los extremos.  Para que se verifique la necesaria propiedad de unisolvencia utilizaremos el espacio de polinomios  de grado inferior o igual a 3. Así pues trabajaremos con aproximaciones polinómicas de grado 3 a trozos y de clase (C^1); ya que le pediremos a estos trozos de cúbicas que "peguen bien" en cada extremo de los subintervalos a considerar, en el sentido de que coincidan en valor tanto de la función como de su derivada primera. A este tipo de funciones cúbicas a trozos se le suele denominar spline cúbico.

Clear["Global`*"]

Definición de las funciones de base del elemento finito de Hermite en el intervalo de referencia [0,1]

\(q\_\(0, 0\)[t_] := \((\(-1\) + t)\)\^2\ \((1 + 2\ t)\)
\)\(q\_\(1, 0\)[t_] := \(-t\^2\)\ \((\(-3\) + 2\ t)\)
\)\(q\_\(0, 1\)[t_] := \((\(-1\) + t)\)\^2\ t
\)(Subscript[q, 1, 1])[t_] := (-1 + t) (t^2)

Problema de contorno de cuarto orden a resolver

Se desea aproximar, mediante el método de los elementos finitos, la solución del siguiente problema de contorno de cuarto orden, relativo al problema de la viga con función de carga sinusoidal:

		(y^(iv  ))) = -sen x,  
      x ϵ [0, π]
y(0) = y(π) = 0
y' (0) = y' (π) = 0

Introducción de los extremos del intervalo y definición de los nodos de la partición

Dominio

Clear[a, b]\(a = 0.; 
b = N[π]; \)

Término independiente

f[x_] := -Sin[x]

Nodos de la partición

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

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

Construcción de las funciones de base del espacio de elementos finitos en un intervalo genérico con nodos (Subscript[x, i])

\(\(w\_0[i_]\)[t_] := 
      N[Which[x[i - 1] ≤ t ≤ x[i], 
          q\_\(1, 0\)[\(t - x[i - 1]\)\/h], 
			
          x[i] ≤ t ≤ x[i + 1], 
          q\_\(0, 0\)[\(t - x[i]\)\/h], True, 0]]
\)\(\(w\_1[i_]\)[t_] := 
      N[Which[x[i - 1] ≤ t ≤ x[i], 
          h\ q\_\(1, 1\)[\(t - x[i - 1]\)\/h], 
			
          x[i] ≤ t ≤ x[i + 1], 
          h\ q\_\(0, 1\)[\(t - x[i]\)\/h], True, 0]]\)

Matriz de rigidez

A = Table[0, {i, 1, 2 (n - 1)}, {j, 1, 2 ( n - 1)}];

For[i = 1,i ≤ n - 1,i++,
	A[[2 i - 1, 2 i - 1]]=(∫_(x[i - 1])^(x[i]))(((((Subscript[q, 1, 0]))''[((x - x[i - 1])/h)])^2)/((h^4)))\[DifferentialD]x+(∫_(x[i])^(x[i + 1]))(((((Subscript[q, 0, 0]))''[((x - x[i])/h)])^2)/((h^4)))\[DifferentialD]x;
	A[[2 i, 2 i]]=
		(∫_(x[i - 1])^(x[i]))(((((Subscript[q, 1, 1]))''[((x - x[i - 1])/h)])^2)/((h^2)))\[DifferentialD]x+(∫_(x[i])^(x[i + 1]))(((((Subscript[q, 0, 1]))''[((x - x[i])/h)])^2)/((h^2)))\[DifferentialD]x;
	A[[2 i - 1, 2 i]]=
		(∫_(x[i - 1])^(x[i]))((((Subscript[q, 1, 0]))''[((x - x[i - 1])/h)] ((Subscript[q, 1, 1]))''[((x - x[i - 1])/h)])/((h^3)))\[DifferentialD]x+(∫_(x[i])^(x[i + 1]))((((Subscript[q, 0, 0]))''[((x - x[i])/h)] ((Subscript[q, 0, 1]))''[((x - x[i])/h)])/((h^3)))\[DifferentialD]x;
A[[2 i, 2 i - 1]] = 
              A[[2 i - 1, 2 i]];
	
	If[i > 1,
		A[[2 i - 1, 
                      2 i - 3]]=(∫_(x[i - 1])^(x[i]))((((Subscript[q, 1, 0]))''[((x - x[i - 1])/h)] ((Subscript[q, 0, 0]))''[((x - x[i - 1])/h)])/((h^4)))\[DifferentialD]x;
		A[[2 i - 1, 
                      2 i - 2]]=(∫_(x[i - 1])^(x[i]))((((Subscript[q, 1, 0]))''[((x - x[i - 1])/h)] ((Subscript[q, 0, 1]))''[((x - x[i - 1])/h)])/((h^3)))\[DifferentialD]x;
	A[[2 i, 
                      2 i - 3]]=(∫_(x[i - 1])^(x[i]))((((Subscript[q, 1, 1]))''[((x - x[i - 1])/h)] ((Subscript[q, 0, 0]))''[((x - x[i - 1])/h)])/((h^3)))\[DifferentialD]x;
	A[[2 i, 
                      2 i - 2]]=(∫_(x[i - 1])^(x[i]))((((Subscript[q, 1, 1]))''[((x - x[i - 1])/h)] ((Subscript[q, 0, 1]))''[((x - x[i - 1])/h)])/((h^2)))\[DifferentialD]x];
	If[i < n - 1,
		A[[2 i - 1, 
                      2 i + 1]]=(∫_(x[i])^(x[i + 1]))((((Subscript[q, 0, 0]))''[((x - x[i])/h)] ((Subscript[q, 1, 0]))''[((x - x[i])/h)])/((h^4)))\[DifferentialD]x;
	A[[2 i - 1, 
                      2 i + 2]]=(∫_(x[i])^(x[i + 1]))((((Subscript[q, 0, 0]))''[((x - x[i])/h)] ((Subscript[q, 1, 1]))''[((x - x[i])/h)])/((h^3)))\[DifferentialD]x;
	A[[2 i, 
                      2 i + 1]]=(∫_(x[i])^(x[i + 1]))((((Subscript[q, 0, 1]))''[((x - x[i])/h)] ((Subscript[q, 1, 0]))''[((x - x[i])/h)])/((h^3)))\[DifferentialD]x;
	A[[2 i, 
                      2 i + 2]]=(∫_(x[i])^(x[i + 1]))((((Subscript[q, 0, 1]))''[((x - x[i])/h)] ((Subscript[q, 1, 1]))''[((x - x[i])/h)])/((h^2)))\[DifferentialD]x]];

(Subscript[A, rig]) = A;

Simplify[(Subscript[A, rig])] // MatrixForm

Obsérvese que sólo se han calculado los elementos que se sabe pueden ser no nulos de la matriz de rigidez, teniendo en cuenta la expresión de cada uno de los tramos de las funciones de base.

Vector de carga

\(d = Table[0, {i, 1, 2\ n - 2}]; 


    For[i = 1, i ≤ n - 1, \(i++\), 
      d[[2\ i - 1]] = 
        NIntegrate[
            f[t]\ q\_\(1, 0\)[\(t - x[i - 1]\)\/h], {t, x[i - 1], x[i]}] + 
          NIntegrate[
            f[t]\ q\_\(0, 0\)[\(t - x[i]\)\/h], {t, x[i], x[i + 1]}]; 
      d[[2\ i]] = 
        NIntegrate[
            f[t]\ h\ q\_\(1, 1\)[\(t - x[i - 1]\)\/h], {t, x[i - 1], x[i]}] + 
          NIntegrate[
            f[t]\ h\ q\_\(0, 1\)[\(t - x[i]\)\/h], {t, x[i], x[i + 1]}]]; 


    d\_car = d; 

MatrixForm[d\_car] // Chop\)

Obsérvese que, en este caso, se ha utilizado la aproximación numérica NIntegrate para el cálculo de las integrales, previendo la dificultad que pudieran tener según la expresión del término independiente de la ecuación diferencial.

Resolución del sistema

c = LinearSolve[(Subscript[A, rig]),  (Subscript[d, car])] // Chop

Solución aproximada

Construcción de la solución aproximada

yy[x_] := 
      (∑_(i = 1)^(n - 1))(
          c[[2 i - 1]] (Subscript[w, 0])[i][x] + 
            c[[2 i]] (Subscript[w, 1])[i][x])

Representación

ga = Plot[yy[t], {t, a, b}, AspectRatio -> 0.25];

Solución exacta: contraste gráfico

Cálculo de la solución exacta

Téngase en cuenta que este bloque sólo se podrá ejecutar si el problema de contorno inicial tiene solución mediante los métodos directos implementados en Mathematica.

Clear[y]
DSolve[{Derivative[4][y][x]==-Sin[x],y[a] == 0,y[b] == 0,y'[a]==0,y'[b]==0},y[x],x]

y[x_] = 
      y[x] /. Chop[%][[1]]

Gráfica de la solución exacta

gb = 
      Plot[y[x], {x, a, b}, PlotStyle -> {Dashing[{0.02, 0.02}]}, 
        AspectRatio -> 0.25];

Contraste gráfico

Show[ga, gb];

Contraste numérico

m = 12;

TableForm[
      Table[{a + (((b - a) i)/m), y[a + (((b - a) i)/m)], 
          yy[a + (((b - a) i)/m)]}, {i, 0, m}], 
      TableHeadings 
        -> {
          None, {"x=a+(b-a)/m i \n", "y[x] (ex.) \n", 
            "yy[x] (apr.) \n"}}]

Ejercicios

1.-  Hallar  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 del elemento finito de Lagrange, empleando  particiones uniformes con 10 y 20 subintervalos, respectivamente. Comparar los resultados obtenidos  con la solución exacta y con los que ya se obtuvieron al resolver este mismo problema mediante la poligonal de Euler.

2.-  Para el problema anterior, aprovechar el hecho de que la ecuación es de coeficientes constantes, y programar el cálculo de tan sólo los elementos de la matriz de rigidez que sean distintos entre sí, así como la resolución del sistema lineal tridiagonal resultante, mediante algún algoritmo específico para este tipo de matrices.

3.-  Programar ahora el mismo proceso pero con el elemento finito de Hermite de clase (𝒞^1). Concretamente, aplicarlos a la resolución aproximada de problemas de cuarto orden como el siguiente:

		\(TraditionalForm
      \`\(\((\(p(x)\) \(y'\)' \((x)\))\)'\)' - \((\(q(x)\) y' \((x)\))\)' + 
          \(r(x)\) \(y(x)\) = f(x), \ x ∈ \([a, b]\)\)
y(a) = (Subscript[y, a]),  y(b) = (Subscript[y, b]),  
      y' (a) = (Subscript[y, a]^((1))),  y' (b) = (Subscript[y, b]^((1)))

considerando los coeficientes p,  q,  r y f suficientemente regulares en el intervalo [a, b].

4.-  Resolver también un problema como el enunciado anteriormente, en el cual las condiciones de contorno son las correspondientes    a  una viga
   - simplemente apoyada: y(a) = y(b) = 0,  
      y'' (a) = y'' (b) = 0
   - empotrada en un extremo y simplemente apoyada en el otro: y(a) = y(b) = 0,  
      y' (a) = y' (b) = 0
   - en voladizo: y(a) = y' (a) = 0,  
      y'' (b) = y''' (b) = 0
empleando también el método de los elementos finitos, con elemento de Hermite de clase (𝒞^1) y una partición uniforme del intervalo con n + 1 nodos equidistantes en total. 
La idea es que ahora, aunque ciertas condiciones de contorno no podrán satisfacerse por nuestra aproximación que sólo será de clase (𝒞^1), aparecerán nuevos grados de libertad en el problema que obligan a considerar más funciones de base en los extremos del intervalo, con lo cual la matriz de rigidez se verá ampliada con nuevas filas y columnas.