Resolución Numérica de
Problemas de Contorno

Método de tiro para problemas de contorno lineales

Consideremos el problema de contorno lineal

		y'' = p(x) y' + q(x) y + r(x),  
      x ∈ [a, b], 
y(a) = α,  
      y(b) = β.
                                          
 Como se indica en el Algoritmo del método de tiro lineal, la resolución de este problema de contorno se reduce a la resolución de dos problemas de valores iniciales y la construcción de una combinación lineal de las mismas. En concreto, se trata de resolver los problemas de valores inciales:

		y'' = p(x) y' + q(x) y + r(x),  
      x ∈ [a, b], 
y(a) = α,  
      y' (a) = 0, 
e
		y'' = p(x) y' + q(x) y,  
      x ∈ [a, b], 
y(a) = 0,  
      y' (a) = 1.
                                                
Si (Subscript[y, 1]) e (Subscript[y, 2]) son las soluciones de estos problemas, respectivamente, entonces se demuestra que la solución del problema inicial es 

		\(TraditionalForm
      \`y(x) = \(y\_1\)(x) + 
          \(\(β - \(y\_1\)(b)\)\/\(\(y\_2\)(b)\)\) \(\(y\_2\)(x).\)\)
 
Ahora bien, las soluciones exactas de los problemas de valores iniciales anteriores sólo pueden obtenerse en algunas ocasiones, de tal forma que en la mayoría de los casos es necesario aplicar un método numérico para obtener una aproximación de cada una. Apliquemos este algoritmo a un ejemplo concreto. Consideremos el problema de contorno lineal:

		y'' = y + y' + sen x,  
      x ∈ [0, π], 
y(0) = 0,  
      y(π) = 1.

En este caso

Clear["Global`*"]

\(a = 0; \ b = π; \ α = 0; \ β = 1; 
p[x_] := 1; \ 
    q[x_] := 1; \ r[x_] := Sin[x]; \)

A continuación se calcula la solución exacta del problema con Mathematica , ayudados por la orden DSolve. En la mayoría de los casos este cálculo deberá omitirse, ya que no será posible llevarlo a cabo de forma exacta.

Clear[y]solexacta = 
      DSolve[{y''[x] == p[x] y'[x] + q[x] y[x] + r[x], 
          y[a] == α, y[b] == β}, y[x], x]grexacta = 
      Plot[y[x] /. solexacta[[1]], {x, a, b}, PlotRange -> All];

Resolución del primer problema de valores iniciales

Comenzamos el método de tiro resolviendo el problema de valores iniciales

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

Para ello, consideramos el problema de valores iniciales para el sistema de ecuaciones diferenciales de primer orden asociado es este problema de segundo orden, es decir, el sistema

		y' = z,   y(a) = α, 
      
z' = p(x) z + q(x) y + r(x),   
      z(a) = 0.
					
Haciendo X=(({{y}, {z}})),  y  f(X)=(({{z}, {p(x) z + q(x) y + r(x)}})), este sistema se puede escribir como el siguiente problema de valores iniciales de orden 1 para la función vectorial X:

		X' = f(X),  X(a)=(({{α}, {0}})).

Una solución numérica se obtiene con sólo aplicar uno de los métodos de resolución aproximada de problemas de valores iniciales estudiados, por ejemplo, el método de Runge-Kutta clásico.

\(f[x_, X_] := {X[\([2]\)], p[x] X[\([2]\)] + q[x] X[\([1]\)] + r[x]}; 

    valor = {{α, 0}}; 
n = 100; 
h = \(b - a\)\/n // N; 

    nodo = Table[a + j\ h, {j, 0, n}]; 


    For[j = 2, j ≤ n + 1, \(j++\), 
	\ \ \ \ \ \ \ 
      k1 = f[nodo[[j - 1]], 
          valor[[j - 1]]]; 

      	\ \ \ \ \ \ \ 
      k2 = f[nodo[[j - 1]] + h\/2, 
          valor[[j - 1]] + 
            \(\(h\ \)\/2\) k1]; 
	\ \ \ \ \ \ \ 
      k3 = f[nodo[[j - 1]] + h\/2, 
          valor[[j - 1]] + h\/2\ k2]; 

      	\ \ \ \ \ \ \ 
      k4 = f[nodo[[j - 1]] + h, 
          valor[[j - 1]] + h\ k3]; 

      	\ \ \ \ \ \ \ 
      AppendTo[valor, 
        valor[[j - 1]] + 
          1\/6\ h\ \((k1 + 2\ k2 + 2\ k3 + k4)\)]
\ \ \ \ \ \ \ ]; \)

Los valores aproximados de la solución y1 en los nodos considerados son:

y1 = Table[valor[[i, 1]], {i, 1, n + 1}];

La gráfica de la solución numérica de este problema de valores inciales es

gry1 = ListPlot[Table[{nodo[[i]], y1[[i]]}, {i, n + 1}]];

Calculamos y representamos a continuación la solución exacta del primer problema de valores iniciales tratado.

Clear[y]DSolve[{y''[x] == p[x] y'[x] + q[x] y[x] + r[x], 
          y[a] == α, y'[a] == 0}, y[x], x] // Chopsol1 = y[x] /. %[[1]];

grsol1 = Plot[sol1, {x, a, b}, PlotRange -> All];

Mostramos conjuntamente las gráficas de la solución exacta y la aproximada de este primer problema de valores iniciales

Show[gry1, grsol1, PlotRange -> All];

Resolución del segundo problema de valores iniciales

Resolvamos ahora el segundo el problema de valores iniciales

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

Para ello, consideramos de nuevo el problema de valores iniciales para el sistema de ecuaciones diferenciales de primer orden asociado es este problema de segundo orden, es decir, el sistema

		y' = z,   y(a) = 0, 
      
z' = p(x) z + q(x) y,   z(a) = 1.
					
Llamando X=(({{y}, {z}})),  y  f(X)=(({{z}, {p(x) z + q(x) y}})), este sistema se puede escribir como el siguiente problema de valores iniciales de orden 1 para la función vectorial X:

		X' = f(X),  X(a)=(({{0}, {1}})).

Del mismo modo,  una solución numérica de este problema se obtiene aplicando uno de los métodos de resolución aproximada de problemas de valores iniciales estudiados, por ejemplo, el método de Runge-Kutta clásico.

\(f[x_, X_] := {X[\([2]\)], p[x] X[\([2]\)] + q[x] X[\([1]\)]}; 

    valor = {{0, 1}}; 


    For[j = 2, j ≤ n + 1, \(j++\), 
	\ \ \ \ \ \ \ 
      k1 = f[nodo[[j - 1]], 
          valor[[j - 1]]]; 

      	\ \ \ \ \ \ \ 
      k2 = f[nodo[[j - 1]] + h\/2, 
          valor[[j - 1]] + 
            \(\(h\ \)\/2\) k1]; 
	\ \ \ \ \ \ \ 
      k3 = f[nodo[[j - 1]] + h\/2, 
          valor[[j - 1]] + h\/2\ k2]; 

      	\ \ \ \ \ \ \ 
      k4 = f[nodo[[j - 1]] + h, 
          valor[[j - 1]] + h\ k3]; 

      	\ \ \ \ \ \ \ 
      AppendTo[valor, 
        valor[[j - 1]] + 
          1\/6\ h\ \((k1 + 2\ k2 + 2\ k3 + k4)\)]
\ \ \ \ \ \ \ ]; \)

Los valores aproximados de la solución y2 en los nodos considerados son:

y2 = Table[valor[[i, 1]], {i, 1, n + 1}];

La gráfica de la solución numérica de este problema de valores inciales es

gry2 = ListPlot[Table[{nodo[[i]], y2[[i]]}, {i, n + 1}]];

Calculamos y representamos a continuación la solución exacta del segundo problema de valores iniciales tratado.

Clear[y]DSolve[{y''[x] == p[x] y'[x] + q[x] y[x], 
          y[a] == 0, y'[a] == 1}, y[x], x] // Chopsol2 = y[x] /. %[[1]];

grsol2 = Plot[sol2, {x, a, b}, PlotRange -> All];

Mostramos conjuntas las gráficas de la solución exacta y la aproximada de este primer problema de valores iniciales

Show[gry2, grsol2, PlotRange -> All];

Solución del problema

Según se indica en este capítulo, la solución numérica del problema de contorno incial se obtiene como combinación de las soluciones de los problemas de valores iniciales anteriormente considerados, según la siguiente expresión:

yaprox = y1 + ((β - y1[[n + 1]])/y2)[[n + 1]] y2;

La gráfica de esta solución numérica es

graprox = 
      ListPlot[Table[{nodo[[i]], yaprox[[i]]}, {i, 1, n + 1}]];

Comparativa gráfica y numérica con la solución exacta

A continuación mostramos simultáneamente la gráfica de la solución exacta y la solución numérica calculada.

Show[grexacta, graprox];

yaprox

\(yexacta = 
      Table[
	N[\((y[x] /. solexacta[\([1]\)])\) /. x -> nodo[\([k]\)]] // 
          Chop
	, {k, n + 1}]\)

Finalmente, calculamos el máximo de los errores en valor absoluto entre la solución exacta y la solución numérica calculada en los nodos considerados.

errormaxTiroLineal = Abs[yaprox - yexacta] // Max

Método de diferencias finitas para problemas lineales

Consideremos el problema de contorno lineal

		y'' = p(x) y' + q(x) y + r(x),  
      x ∈ [a, b], 
y(a) = α,  
      y(b) = β.

Apliquemos el algoritmo descrito como Algoritmo del método de diferencias finitas lineal al problema de contorno lineal resuelto en la primera sección de este notebook, es decir, al problema:

		y'' = y + y' + sen x,  
      x ∈ [0, π], 
y(0) = 0,  
      y(π) = 1.

Entrada de datos

\(a = 0.; \ b = Pi; \ α = 0.; \ β = 1.; 
p[x_] := 1; \ 
    q[x_] := 1; \ r[x_] := Sin[x]; 
n = 100; 
h = \(b - a\)\/n; 

    nodo = Table[a + j\ h, {j, 0, n}]; \)

Matriz de coeficientes

\(coef = Table[0., {i, n + 1}, {j, n + 1}]; 
coef[\([1, 1]\)] = 1.; \ 
    coef[\([n + 1, n + 1]\)] = 1.; 


    For[i = 2, i ≤ n, \(i++\), 
\ \ \ \ \ \ \ 
      coef[\([i, i]\)] = \(-2\) - \(h\^2\) q[nodo[\([i]\)]]; 
      
\ \ \ \ \ \ 
      coef[\([i, i + 1]\)] = 1 - \(h\/2\) p[nodo[\([i]\)]]; 
      
\ \ \ \ \ 
      coef[\([i, i - 1]\)] = 1 + \(h\/2\) p[nodo[\([i]\)]]]\)

Términos independientes

\(ind = Table[\(h\^2\) r[nodo[\([i]\)]], {i, n + 1}]; 

    ind[\([1]\)] = α; 
ind[\([n + 1]\)] = β; \)

Resolución del sistema

yaproxDF = LinearSolve[coef, ind];

graproxDF = 
      ListPlot[Table[{nodo[[i]], yaproxDF[[i]]}, {i, n + 1}]];

Solución exacta

Clear[y]DSolve[{y''[x] == p[x] y'[x] + q[x] y[x] + r[x], 
          y[a] == α, y[b] == β}, y[x], x] // Chopsolexacta = y[x] /. %[[1]];

grexacta = Plot[solexacta, {x, a, b}, PlotRange -> All];

\(yexacta = 
      Table[
	N[solexacta /. x -> nodo[\([k]\)]] // Chop
	, {k, n + 1}]
        \)

Comparativa gráfica y numérica con la solución exacta

A continuación mostramos simultáneamente la gráfica de la solución exacta y la solución numérica calculada.

Show[grexacta, graproxDF];

También comparamos las aproximaciones obtenidas mediante el método de tiro lineal y las que se acaban de obtener mediante el método de diferencias finitas.

yaprox - yaproxDF

Finalmente, calculamos el máximo de los errores en valor absoluto entre la solución exacta y la solución numérica calculada en los nodos considerados.

errormaxDF = Abs[yaproxDF - yexacta] // Max

Ejercicios

1.- Consideremos el siguiente problema de contorno:

		y'' + 2  y' - 3  y = 9  x,   
          x ∈ (0, 1)
y(0) = 0,    y(1) = 2.

a) Halle su solución exacta mediante la orden DSolve o a mano.
b) Resuélvalo numéricamente mediante el método de tiro lineal a partir de una partición uniforme del intervalo [0, 1] con 11 nodos, y calcule el error de tal aproximación como el máximo de los errores absolutos en los nodos.
c) Resuelva las mismas cuestiones aplicando el método de diferencias finitas.
d) ¿Qué método ofrece mejores resultados? Realice un constraste gráfico de ambos métodos entre sí y de éstos con la solución exacta.