Métodos de Ritz y Galerkin para problemas de contorno

Utilizaremos el método de Ritz, con funciones coordenadas de tipo polinómico o trigonométrico,  para aproximar las extremales de un funcional

		\(TraditionalForm
      \`F[y] = ∫\_a\%b\( F(x, y(x), y' \((x)\))\) dx\)
		
que satisfacen ciertas condiciones de contorno de tipo Dirichlet, no necesariamente homogéneas, en los extremos del intervalo: y(a) = (Subscript[y, a]),  y(b) = (Subscript[y, b]).                                      
Cuando esta cuestión de Cálculo de Variaciones proviene de la búsqueda de soluciones de un problema de contorno autoadjunto

		\(TraditionalForm
      \`\(-\((\(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), entonces también es posible utilizar el método de Galerkin, que es un método general que puede ser aplicado cuando el problema no es autoadjunto. Por lo tanto, se puede aplicar este método de Galerkin al problema

		\(TraditionalForm
      \`L[y] = \(\(\(a\_0\)(x)\) \(y'\)' \((x)\) + 
            \(\(a\_1\)(x)\) y' \((x)\) + \(\(a\_2\)(x)\) \(y(x)\) = f(x)\), \ 
      x ∈ \([a, b]\)\)
y(a) = (Subscript[y, a]),  y(b) = (Subscript[y, b])
		
Se supone que las funciones (Subscript[a, 0]),  (Subscript[a, 1]) y (Subscript[a, 2]) son continuas en [a, b].

Ejemplo

Resolveremos a continuación  el  problema de las extremales para el siguiente funcional cuadrático

		F[
          y] = (∫_a^b)(y' (((x))^2) + ((y(x))^2)) dx
		
satisfaciendo las  condiciones de contorno y(0) = 1,  y(1) = 3.

Clear["Global`*"]

Definición de las condiciones de frontera tipo Dirichlet

a = 0.;  ya = 0.;  b = π // N;   yb = 0.;

Elección del tipo de funciones coordenadas

Para elegir un tipo concreto bastará con ejecutar sólo una de las celdillas que a continuación definen las  funciones coordenadas, polinómicas o de tipo trigonométrico, satisfaciendo unas condiciones de tipo Dirichlet homogéneas.

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
          ))

Método de Ritz

Definición del funcional asociado

Clear[x, y, dy, F, FF]F[x_, y_, dy_] := (dy^2) + (y^2) - 2  Sin[x] yFF[y_] := 
      (∫_a^b) F[x, y, (∂_x) y] \[DifferentialD]x

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

Entrada del orden de aproximación

m = 3;

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

Paso a condiciones de frontera homogéneas

Utilizamos el cambio \(TraditionalForm
      \`z(x) = y(x) - \((y\_a + \((x - a)\) y\_b - y\_a\/\(b - a\))\)\).

α[x_] := ya + (x - a) ((yb - ya)/(b - a))(Subscript[y, m])[x_] :=  (Subscript[z, m])[x]  + α[x]

Función cuyos puntos críticos se calcularán

G = Simplify[FF[(Subscript[y, m])[x]]]

ptoscrit = 
      NSolve[Table[(∂_(c[i]))G == 0, {i, m}], Table[c[i], {i, m}]]

(Subscript[z, m])[x_] := 
      (∑_(i = 1)^m) c[i] ϕ[x, i] /. 
        ptoscrit[[1]](Subscript[y, m])[x_] := (Subscript[z, m])[x] + α[x]

Collect[(Subscript[y, m])[x], x]

vaprox1 = FF[(Subscript[y, m])[x]] // Chop

Método de Galerkin para problemas autoadjuntos

Definición de las funciones que determinan el problema de contorno

Clear[x, p, q, f]p[x_] := 1;  q[x_] := 1;  f[x_] := -Sin[x];

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

Paso a condiciones de frontera homogéneas y definición del nuevo término independiente

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

Planteamiento y resolución del sistema lineal

m = 3;

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++, 
          
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; 
            
matA[[j, i
                ]] = 
              matA[[i, j]]]; 
          
matA[[i, i
              ]] = 
            (∫_a^b)
                (p[x] dϕ[x, i]^2 + q[x] ϕ[x, i]^2) 
                \[DifferentialD]x];

MatrixForm[matA]

\(\(matd = 
      Table[∫\_a\%b f\&^[x]\ ϕ[x, i] \[DifferentialD]x, {i, 
          m}]; \)\)

c=LinearSolve[matA,matd]

Construcción de la solución aproximada como combinación lineal de las funciones coordenadas

(Subscript[z, m])[x_] := 
      (∑_(i = 1)^m) c[[i]] 
          ϕ[x, i](Subscript[y, m])[x_] := (Subscript[z, m])[x] + α[x]

(Subscript[y, m])[x]

Collect[(Subscript[y, m])[x], x]

vaprox2 = FF[(Subscript[y, m])[x]]

% // N

Solución exacta

Si conocemos la solución exacta del  problema de extremales, calculada resolviendo la  correspondiente ecuación de Euler, podemos ahora compararla con las distintas soluciones aproximadas. Compruebe que la extremal que satisface las condiciones de contorno especificadas en el ejemplo propuesto es la siguiente:

\(y[x_] := 
      \(\((1 - 3\ E)\)\ E\^x + 
          \((3\ E - E\^2)\)\ 
            E\^\(-x\)\)\/\(1 - E\^2\)
\)vexacto = Simplify[FF[y[x]]]

Comparación de los valores obtenidos para las distintas aproximaciones y el valor exacto

SetPrecision[{vaprox1, vaprox2, vexacto}, 10]

Método de Galerkin para problemas no autoadjuntos

Definición del operador diferencial asociado

\(Clear[x, y, a0, a1, a2, f]

\)a0[x_] := 1a1[x_] := 0a2[x_] := -1f[x_] := 0\(L[y_] := 
      a0[x]\ ∂\_{x, 2}y + a1[x]\ ∂\_x y + a2[x]\ y\)

Paso a condiciones de frontera homogéneas y definición del nuevo término independiente

\(α[x_] := ya + \((x - a)\) \(yb - ya\)\/\(b - a\)
\)\(f\&^[x_] := f[x] - L[α[x]]\)

Planteamiento del sistema lineal

m = 3;

matA = 
      Table[(∫_a^b) L[ϕ[x, i]] ϕ[x, j] 
            \[DifferentialD]x, {i, m}, {j, m}];

\(\(matd = 
      Table[∫\_a\%b f\&^[x]\ ϕ[x, i] \[DifferentialD]x, {i, 
          m}]; \)\)

c = LinearSolve[matA, matd]

Construcción de la solución aproximada como combinación lineal de las funciones coordenadas

(Subscript[z, m])[x_] := 
      (∑_(i = 1)^m) c[[i]] 
          ϕ[x, i](Subscript[y, m])[x_] := (Subscript[z, m])[x] + α[x]

(Subscript[y, m])[x]

Collect[(Subscript[y, m])[x], x]

N[%]

Ejercicios

1.-  Resolver aproximadamente mediante alguno de los métodos de Ritz o de Galerkin, con dos funciones coordenadas de tipo polinómico, el problema de optimización asociado al funcional

 		\(TraditionalForm
      \`F[y] = ∫\_0\%1
            \((y' \((x)\)\^2 - \(y(x)\)\^2 - 2  x\ \ \(y(x)\))\) dx\)

junto con las condiciones de contorno homogéneas y(0) = y(1) = 0. Comparar  también estas soluciones aproximadas con la extremal que se  obiene resolviendo la correspondiente ecuación de Euler.

2.-  Resolver también numéricamente, con dos y tres  funciones coordenadas de tipo polinómico, el problema de optimización asociado al funcional

 		\(TraditionalForm
      \`F[y] = ∫\_0\%1\((
            y' \((x)\)\^2 + \(y(x)\)\^2 + 2  x\ \ \(y(x)\))\)\)

junto con condiciones de contorno no homogéneas y(0) = 1,  y(1) = 2. Comparar  también estas soluciones aproximadas con la solución exacta.

3.-  Usar  indistintamente los métodos de Ritz o de Galerkin con dos y tres funciones coordenadas de tipo polinómico para resolver de forma aproximada los siguientes problemas de contorno autoadjuntos:
	
(a)		y'' (x) - (1 + x) y(x) = 1,  
      x ∈ [0, 1]
y(0) = 1,  y(1) = 2
	
(b)		\(TraditionalForm
      \`\((\((x\^2 + 1)\) y' \((x)\))\)' - 
          \(e\^\(2  x\)\) \(y(x)\) = 1, \ 
      x ∈ \([0, 1]\)\)
y(0) = 0,  y(1) = 1

4.-  Resolver también aproximadamente los siguientes problemas de contorno, dados en  forma no autoadjunta:				
	
(a)		\(TraditionalForm
      \`3  x\ \(y'\)' \((x)\) + 3  x\ \(y'\)' \((x)\) + 3  y' \((x)\) - 
          7 \( x\^3\) \(y(x)\) - 4  x\^2 = 0, \ x ∈ \([0, 1]\)\)
y(0) = 1,  y(1) = 2

(b)		\(TraditionalForm
      \`\(x\^2\) \(y'\)' \((x)\) + y' \((x)\) - x\ \(y(x)\) = 0, \ 
      x ∈ \([2, 5]\)\)
y(2) = 1,  y(5) = 10

5.- Realizar las adaptaciones necesarias en el código de Mathematica anterior para poder también resolver problemas de contorno con condiciones de contorno de tipo Neumann o mixtas, no necesariamente homogéneas, y, a continuación, resolver numéricamente los siguientes problemas:
	
(a)		y'' (x) - (1 + x) y(x) = 1,  
      x ∈ [0, 1]
y' (0) = 1,  y' (1) = 2

(b)		\(TraditionalForm
      \`\(x\^2\) \(y'\)' \((x)\) + y' \((x)\) - x\ \(y(x)\) = 0, 
      x ∈ \([2, 5]\)\)
y(2) = 1,  y' (5) = 0