Resolución aproximada
de problemas de Sturm-Liouville mediante
métodos variacionales

Se considera el problema de valores propios con ecuación diferencial

		(p(x) y' (x))' - q(x) y(
              x) + λr(x) y(x) = 
        0,  x∈[a, b]

y condiciones de contorno homogéneas de uno de los siguientes tipos:
   - Dirichlet:     y(a) = y(b) = 0.
   - Neumann:   y' (a) = y' (b) = 0.
   - Mixtas:       y(a) = y' (b) = 0 o y' (a) = y(b) = 0.

Suponemos que las funciones q y r son al menos continuas y p derivable con continuidad en el intervalo [a, b], con p(x),  r(x) > 0, aunque estas hipótesis se pueden debilitar.
        
Aproximaremos numéricamente los primeros valores propios, junto con las funciones propias asociadas, mediante los conocidos métodos de Ritz y Galerkin utilizando funciones coordenadas de tipo polinómico, o bien trigonométrico, en el intervalo [a, b].

A continuación también veremos la aplicación de estos mismos métodos, pero con funciones splines lineales, que son funciones  polinómicas de grado  uno  a trozos y continuas.

Ejemplo

Resolveremos numéricamente el problema de valores propios

		y'' (x) + λy(x) = 0,  
      x ∈ [0, π]
y(0) = y(π) = 0

cuyos valores propios exactos también se pueden calcular, y son (n^2), n = 1, 2,  ... ; las funciones propias normalizadas, es decir, las que verifican que (∫_0^π) r(x)  (((Subscript[φ, n])(x))^2)dx=1, son

		(Subscript[φ, n])(
          x) = (Sqrt[(2/π)]) sen(n x),  n ≥ 1.

Clear["Global`*"]

Definición de los extremos del intervalo

a = 0; b = π;

Definición del funcional asociado

\(Clear[x, p, q, r]
\)p[x_] := 1q[x_] := 0r[x_] := 1

\(Clear[y, dy, F, FF, FR, RR]
\)\(F[x_, y_, dy_] := p[x]\ dy\^2 + q[x]\ y\^2
\)FF[y_] := 
      (∫_a^b) F[x, y, (∂_x) y] \[DifferentialD]x

Métodos de Ritz y Galerkin con funciones polinómicas o trigonométricas

Definición de las funciones coordenadas

Para este ejemplo, utilizaremos funciones coordenadas de tipo trigonométrico, ya que como veremos, son las más apropiadas en este caso. Se incluye el código correspondiente a las de tipo polinómico para cada uno de los tipos de condiciones de cortono.

Que satisfacen unas condiciones de frontera homogéneas de tipo Dirichlet: y(a) = y(b) = 0

de tipo polinómico

\(Clear[x, \ k, \ ϕ, \ dϕ]
\)\(ϕ[x_, k_] := \ \((x - a)\) \((b - x)\) x\^\(k - 1\)
\)dϕ[x_, k_] := 
      (b - x) (x^(-1 + k)) + 
        (-1 + k) (b - x) (x^(-2 + k)) (-a + x) - 
        (x^(-1 + k)) (-a + x)

de tipo trigonométrico

\(Clear[x, \ k, \ ϕ, \ dϕ]
\)\(ϕ[x_, k_] := Sin[k\ π\/\(b - a\)\ \((x - a)\)]
\)dϕ[x_, k_] := 
      ((k π 
          Cos[((k π (-a + x))/(-a + b))])/(-a + b
          ))

Que satisfacen las condiciones de frontera mixtas y' (a) = y(b) = 0

\(Clear[x, k, ϕ, dϕ]
\)\(ϕ[x_, k_] := \((b - a)\)\^\(k + 1\) - \((x - a)\)\^\(k + 1\)
\)dϕ[x_, k_] = -((1 + k) (((-a + x))^k));

Que satisfacen las condiciones de frontera mixtas del y(a) = y' (b) = 0

\(Clear[x, k, ϕ, dϕ]
\)\(ϕ[x_, k_] := \((b - a)\)\^\(k + 1\) - \((b - x)\)\^\(k + 1\)
\)dϕ[x_, k_] = (1 + k) (((b - x))^k);

Método de Ritz

Entrada del orden de aproximación

m = 3;

Resolución aproximada utilizando un orden de aproximación cualquiera

\(Clear[c]
\)(Subscript[y, m])[x_] := (∑_(i = 1)^m) c[i] ϕ[x, i]

Definición de la función de varias variables a la que calcularemos sus puntos críticos, junto con la ligadura de tipo integral

\(Clear[G, R]
\)\(cc = Table[c[i], {i, m}]; 

G = Collect[FF[y\_m[x]], cc]
\)R = 
      Collect[(∫_a^b) r[x] (Subscript[y, m])[x]^2 \[DifferentialD]x, cc]

\(Clear[H, λ]
\)\(H = G - λ\ R; 


    ec = Table[∂\_\(c[i]\)H == 0, {i, m}]; 


    ecs = AppendTo[ec, R == 1]; 

var = Table[c[i], {i, m}]; 


    vars = AppendTo[var, λ]; 


    ptoscrit = Chop[Solve[ecs, vars]]\)

Nótese que cada valor propio aparece dos veces, asociado a dos funciones propias normalizadas, siendo una la opuesta de la otra. Esto es debido a que los funcionales que intervienen son todos cuadráticos, por lo que esta función propia normalizada estará univocamente determinada salvo el signo.

Hay que hacer notar también que obtendremos las soluciones  en forma simbólica, si empleamos la sentencia  Solve, en vez de  NSolve, que nos la proporcionará en forma decimal. De todas formas habrá ocasiones en las que nos interese más esto, ya que el tiempo de ejecución será menor en un caso que en otro. Por ello, bastará con que cambiemos las entencia correspondiente, según nos convenga.

Collect[(Subscript[y, m])[x], x] /. ptoscrit

vaprox = N[FF[(Subscript[y, m])[x]] /. ptoscrit]

Método de Galerkin

Planteamiento del sistema lineal para (Subscript[c, 1]), (Subscript[c, 2]) , ..., (Subscript[c, m])

\(Clear[m, λ]
\)m = 5;

Aquí vamos a aprovechar el hecho de que se trata de un problema autoadjunto y, por tanto la matriz que obtendremos deberá ser simétrica. Así, no vamos a tener que calcular todos los elementos de la misma, sino sólo los que están a un lado de la diagonal principal. Con esto se consigue cierto ahorro en el tiempo de cálculo, sobre todo cuando el orden de la matriz es grande.

\(\(matA = Table[0, {i, m}, {j, m}]; 


    For[i = 1, i ≤ m, \(i++\), 
      
matA[[i, i]] = 
        ∫\_a\%b\((p[x]\ dϕ[x, i]\^2 + q[x]\ ϕ[x, i]\^2)\) 
              \[DifferentialD]x - 
          λ\ 
            \(∫\_a\%b r[x]\ \(ϕ[x, i]\^2\) 
                \[DifferentialD]x\); 
	
      For[j = 1, j < i, \(j++\), 
		
        matA[[i, j]] = 
          ∫\_a\%b
                \((p[x]\ dϕ[x, i]\ dϕ[x, j] + 
                    q[x]\ ϕ[x, i]\ ϕ[x, j])\) 
                \[DifferentialD]x
                  
						\ \ \ \ \ \ \ \  - 
            λ\ 
              \(∫\_a\%b r[x]\ ϕ[x, i]\ ϕ[x, j] 
                  \[DifferentialD]x\); 
									
		
        matA[[j, i]] = 
          matA[[i, j
            ]]
	]
]
      															\)\)

Además, si no necesitamos obtener el valor exacto de los elementos de dicha matriz, de forma simbólica, o si la memoria de nuestro ordenador no nos lo permitiera, para un orden de aproximación elevado; lo más conveniente sería cambiar las sentencias de integración simbólica  Integrate, por las correspondientes fórmulas de integración numérica  NIntegrate. Aunque en este último caso también se puede llegar a producir ciertas imprecisiones debido a las grandes oscilaciones que pueden tener las funciones del integrando. De todas formas, aunque aparezcan varios mensajes en rojo advirtiéndonos de esta situación, los resultados obtenidos suelen ser aceptables.

MatrixForm[matA]

det = Det[matA]\(
\)NSolve[det == 0, λ]

Vemos que en este caso, y de forma excepcional, hemos obtenido los valores propios exactos al emplear funciones coordenadas de tipo trigonométrico. Ésto, como es normal, no va a ocurrir en general, y es debido a que en este ejemplo las funciones propias coinciden precisamente con estas funciones coordenadas.

Empleo de funciones splines lineales

Veremos a continuación la resolución de este mismo problema del ejemplo anterior, empleando ahora funciones splines lineales; que equivale al método de la poligonal de Euler cuando se emplee la técnica de Ritz, o incluso al método de los elementos finitos con elemento finito de Lagrange; también resulta muy parecido al método de diferencias finitas usual, cuando se aproxima la derivada segunda por la correspondiente diferencia dividida centrada.

Determinación de la partición uniforme

\(Clear[n]
\)n = 6

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

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

Definición de las funciones de base del elemento finito de Lagrange

\(p\_1[x_] := 1 - x
\)(Subscript[p, 2])[x_] := x

\(Clear[k, t, ϕ, dϕ]
\)\(ϕ[t_, k_] := 
      Which[x[k - 1] ≤ t ≤ x[k], \ 
        p\_2[\(t - x[k - 1]\)\/\(x[k] - x[k - 1]\)]\ , 
        x[k] ≤ t ≤ x[k + 1], \ 
        p\_1[\(t - x[k]\)\/\(x[k + 1] - x[k]\)], True, 0]
\)dϕ[t_, k_] := 
      Which[x[k - 1] ≤ t ≤ x[k], 
        -  (1/(x[k] - x[k - 1])), 
        x[k] ≤ t ≤ x[k + 1],  (1/(x[k + 1] - x[k])), 
        True, 0]

Método de los elementos finitos (Galerkin) para el cálculo aproximado de los valores propios

\(Clear[λ]
\)\(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] - λ\ r[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] - λ\ r[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] - λ\ r[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; 

MatrixForm[A\_rig]\)

\(det = Det[A\_rig]; 


    valpropios = NSolve[det == 0, λ]\)

Método de la poligonal de Euler (técnica de Ritz)  para el cálculo aproximado de las funciones propias

Definición de las condiciones de frontera de tipo Dirichlet homogéneas:

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

y[n] = 0; \)

Cálculo de la función de varias variables con ligaduras que se optimizará

\(Clear[G, R]
\)\(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]
\)R = 
      Simplify[(∑_(i = 1)^n)
            (∫_(x[i - 1])^(x[i]))r[x] 
              (((y[i - 1] + 
                      (((x - x[i - 1]) (y[i] - y[i - 1]))/h)))^2)
                   \[DifferentialD]x]

Empleo del método de los multiplicadores de Lagrange

\(Clear[H, λ]
\)H = G - λ R

\(ec = Table[∂\_\(y[i]\)H == 0, {i, n - 1}]; 


    ecs = AppendTo[ec, R == 1]; 

var = Table[y[i], {i, n - 1}]; \)

\(For[k = 1, k ≤ Length[valpropios], \ \(k++\), 
      
	
      λ = λ /. \(valpropios[\([k]\)]\)[\([1]\)]; 
      
funcprop[k] = \ Chop[NSolve[ecs, var]]; 
      
	]\)

Comparación de alguna de estas funciones propias aproximadas con las exactas (en caso de conocer éstas)

Veamos a continuación la comparación gráfica, con la solución exacta, de alguna de las funciones propias aproximadas calculadas anteriormente, y que se encuentran almacenadas en las listas  funcprop[k].

\(\(k = 2; 


    y\&^[t_] := 
      ∑\+\(i = 1\)\%\(n - 1\)\ y[i]\ \ ϕ[t, i] /. 
        \(funcprop[k]\)[\([1]\)]; 


    gaprox = Plot[\(-y\&^[t]\), {t, a, b}, DisplayFunction -> Identity]; 
    

gexact = 
      Plot[\@\(2\/π\)\ Sin[k\ x], {x, a, b}, 
        DisplayFunction -> Identity]; 


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

Nótese que, en general, cuando se emplean splines lineales se obtienen aproximaciones mucho mejores que cuando se toman funciones coordenadas de tipo polinómico para este mismo tipo de problemas; no obstante, estas aproximaciones son bastante más aceptables para los primeros valores propios y funciones propias que para los siguientes, degradándose los resultados a medida que comparamos las sucesivas aproximaciones con las soluciones exactas. También hay que hacer notar que, como para el problema de Dirichlet,  las funciones propias normalizadas asociadas a cada uno de los valores propios vienen determinadas unívocamente salvo el signo, alguna de las aproximaciones de las funciones propias corresponde a uno de los signos de la función propia exacta, mientras que otras lo son para el signo opuesto de otra de las funciones propias del problema, cuestión que tendrá fácil arreglo sin más que cambiar de signo los resultados obtenidos según cual de las dos funciones propias queramos aproximar, una o su opuesta.

Ejercicios

1.- Usar el método de Ritz, con funciones coordenadas de tipo polinómico,  para aproximar los dos primeros valores propios y funciones propias asociadas,  de los siguientes problemas de Sturm-Liouville:
a)   \(TraditionalForm
      \`\(y'\)' \((x)\) + \(λ(1 + x\^2)\) \(y(x)\) = 0, \ 
      x ∈ \([0, 1]\), \ y(0) = \(y(1) = 0\)\).
b)  (x y' (x))' + λ y(x) = 0,   
      x ∈ [0, 1],  y' (0) = y(1) = 0.
c) y'' (x) + λ (x^2) y(x) = 0,   
      x ∈ [0, 1],  y(0) = y' (1) = 0.

2.-  Seguidamente, utilizar también el método de Galerkin  para aproximar los cuatro primeros valores propios de los problemas de Sturm-Liouville anteriores.

3.- Para los siguientes problemas de valores propios, calcular también lo que se obtendría utilizando funciones coordenadas de tipo trigonométrico, para poder compararlos con los calculados mediante funciones de tipo polinómico.
 a)   y'' (x) + λ x y(x) = 0,  
      x ∈ [0, 1],  y(0) = y(1) = 0.
 b)  (x^2) y'' (x) + λ y(x) = 0, 
       x ∈ [1, (e^2)],  y(1) = y((e^π)) = 0.
 (Hallar también los valores exactos en este caso).

4.- Implementar ahora el método de los trapecios para el cálculo aproximado de las integrales involucradas en el cálculo de la función con ligaduras que se va a optimizar, cuando se quiere aproximar las funciones propias asociadas a cada uno de los valores propios (para ello el lector se puede inspirar en la práctica correspondiente al método de la poligonal de Euler).