Problemas de Contorno I:
Métodos de tiro 
y de diferencias finitas

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 la sección 3 y, más concretamente, en el algoritmo descrito en la misma con el nombre de 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 

		y(
          x) = (Subscript[y, 1])(
            x) + ((β - (Subscript[y, 1])(b))/((Subscript[y, 2])(b))) (Subscript[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]DSolve[{y''[x] == p[x] y'[x] + q[x] y[x] + r[x], 
        y[a] == α, y[b] == β}, y[x], x]

solex = y[x] /. %[[1]]

grex = Plot[solex, {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.

Clear[x, X, f]\(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}], 
        PlotStyle -> PointSize[0.01]];

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]

sol1 = 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.

Clear[x, X, f]\(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}], 
        PlotStyle -> PointSize[0.01]];

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]

sol2 = y[x] /. %[[1]]

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

Mostramos las gráficas de la solución exacta y la aproximada de este segundo 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:

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

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

graprox = 
      ListPlot[Table[{nodo[[i]], y[[i]]}, {i, 1, n + 1}], 
        PlotStyle -> PointSize[0.015]];

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

Show[grex, graprox];

Método de tiro para problemas de contorno no lineales

Consideramos ahora el caso de un problema de contorno de segundo orden, no necesariamente lineal. Es decir, se trata del  problema de contorno 

		y'' = f(x, y, y'),   x ∈ [a, b], 
      
y(a) = α,   y(b) = β.
					
Para su resolución numérica vamos a implementar el algoritmo que aparece en el capítulo con el título de Algoritmo del método de tiro no lineal (secante).

Como ejemplo, vamos a resolver el problema de contorno no lineal

		y'' = Log[y] - y' + 2  (e^x) - x,   
      x ∈ [0, 1], 
y(0) = 1,   
      y(1) = e, 
		
cuya solución exacta es y = (e^x).

Por tanto, la implementación del método de tiro para el problema de contorno no lineal dado es la siguiente.

Entrada de datos

\(Clear["Global`*"]
\)f[x_, X_] := {X[[2]], 
        Log[X[[1]]] - X[[2]] + 2.  (E^x) - x}
        
\(a = 0; b = 1; 

α = 1; \ 
    β = E; 

n = 100; 

m[0] = 0.; m[1] = 0.1; 


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


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

niter = 10; \ tol = 10\^\(-6\); 
    \)

Resolvemos, en primer lugar el sistema de ecuaciones diferenciales asociado a este problema de contorno de segundo orden para la pendiente m[0].

Primera aproximación

\(valor = {{α, m[0]}}; 


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

    
\)yy[0] = Table[valor[[i, 1]], {i, n + 1}];

Segunda aproximación

Del mismo modo, resolvemos el sistema de ecuaciones diferenciales asociado al problema para la pendiente m[1].

\(valor = {{α, m[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)\)]
\ \ \ \ \ \ \ ]
      
\)yy[1] = Table[valor[[i, 1]], {i, n + 1}];

Aproximaciones sucesivas

Por último, mediante el método de la secante, se calculan nuevas pendientes y nuevas soluciones numéricas del problema.

\(For[k = 2, k ≤ niter, \(k++\), 
      
\ \ \ \ \ \ \ 
      m[k] = m[k - 1] - 
          \(\((\(yy[k - 1]\)[\([n + 1]\)] - β)\) 
              \((m[k - 1] - m[k - 2])\)\)\/\(\(yy[k - 1]\)[\([n + 1]\)] - 
              \(yy[k - 2]\)[\([n + 1]\)]\); 
\ \ \ \ \ \ \ 
      valor = {{α, m[k]}}; 
\ \ \ \ \ \ \ \ \ 
      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)\)]

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
\ \ 
      yy[k] = Table[valor[\([i, 1]\)], {i, n + 1}]; 
      
\ \ \ \ \ \ \ 
      If[Abs[\(yy[k]\)[\([n + 1]\)] - N[β]] < tol, Break[]]]; 

    


    If[k == niter + 1, yaprox = yy[n + 1], yaprox = yy[k]]; \)

graprox = 
      ListPlot[Table[{nodo[[i]], yaprox[[i]]}, {i, 1, n + 1}], 
        PlotStyle -> PointSize[0.015]];

La gráfica de la solución exacta es

grex = Plot[(E^x), {x, a, b}];

Show[graprox, grex];

Téngase en cuenta, a la hora de aplicar este método, que no hemos hecho un estudio de las condiciones de convergencia del método de la secante aplicado, así como  la gran inestabilidad del método de tiro en sí, por lo que aconsejamos que, si ocurren errores en la resolución de otro problema que se proponga, revise los valores iniciales de la pendiente, es decir, los valores de m[0] y m[1].

Método de diferencias finitas para problemas lineales

Cosideremos 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 en el capítulo bajo el título de 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

Clear["Global`*"]

\(a = 0; \ b = π; \ α = 0; \ β = 1; 

p[x_] := 1; \ 
    q[x_] := 1; \ r[x_] := Sin[x]; 

n = 100; 

h = \(b - a\)\/n // 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

yaprox = LinearSolve[coef, ind];

Contraste gráfico

graprox = 
      ListPlot[Table[{nodo[[i]], yaprox[[i]]}, {i, n + 1}], 
        PlotStyle -> PointSize[0.015]];

Clear[y]DSolve[{y''[x] == p[x] y'[x] + q[x] y[x] + r[x], 
        y[a] == α, y[b] == β}, y[x], x]

solex = y[x] /. %[[1]]

grex = Plot[solex, {x, a, b}, PlotRange -> All];

Show[grex, graprox];

Método de diferencias finitas para problemas no lineales

Consideremos el problema de contorno no lineal

		y'' = f(x, y, y'),   x ∈ [a, b], 
      
y(a) = α,   y(b) = β.
					
El método de diferencias finitas aplicado a este problema se reduce a resolver el sistema de ecuaciones:

		(Subscript[y, 0]) = α, 
      
(((Subscript[y, i + 1]) - 2  (Subscript[y, i]) + (Subscript[y, i - 1]))/(h^2)) = 
        f((Subscript[x, i]), (Subscript[y, i]), (((Subscript[y, i + 1]) - (Subscript[y, i - 1]))/(2  h))),   
      i = 1, ..., n, 
(Subscript[y, n]) = β, 

siendo a = (Subscript[x, 0]) < (Subscript[x, 1]) < ... < (Subscript[x, n]) = b una partición uniforme del intervalo [a, b] de paso h.

Apliquemos el método al ejemplo de la sección 2 de este notebook,

		y'' = Log[y] - y' + 2  (e^x) - x,   
      x ∈ [0, 1], 
y(0) = 1,   
      y(1) = e, 
		
cuya solución es y = (e^x).

Entrada de datos

\(f[x_, y_, z_] := Log[y] - z + 2  E\^x - x; 

a = 0; \ 
    b = 1; 

α = 1; β = E; 
    

n = 100; 

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

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

niter = 10; \)

Resolución del sistema por el método de Newton-Raphson

Aplicamos el método de Newton-Raphson para sistemas de ecuaciones lineales cuyo algoritmo se puede encontrar en el Capítulo 6 y la práctica correspondiente.
Comenzamos introduciendo en la variable lista las funciones coordenadas de la que define el sistema de ecuaciones no lineales que se obtiene en este caso.

\(lista = 
      Table[\(y[i + 1] - 2  y[i] + y[i - 1]\)\/h\^2 - 
          f[nodo[\([i]\)], y[i], \(y[i + 1] - y[i - 1]\)\/\(2  h\)], {i, 2, 
          n}]; 

PrependTo[lista, y[1] - α]; 


    AppendTo[lista, y[n + 1] - β]; \)

A continuación se calcula el jacobiano de esta función y se aplica el algoritmo de Newton-Raphson considerando como criterios de parada, o bien el número de interaciones niter, o bien que la diferencia de dos aproximaciones sucesivas tenga norma menor que una cierta tolerancia.

\(j = Table[D[lista[\([i]\)], y[j]], {i, n + 1}, {j, n + 1}]; 
    

c[0] = Table[1., {i, n + 1}]; 


    tolerancia = 10^\((\(-10\))\); 

iter = 0; 


    For[i = 1, i ≤ niter, \(i++\), 
\ \ \ \ 
      d = lista /. Table[y[j] -> \(c[i - 1]\)[\([j]\)], {j, n + 1}]; 
      
\ \ \ \ \ \ \ If[\@\(d.d\) < tolerancia, Break[]]; 
      
\ \ \ \ \ \ \ 
      u[i] = LinearSolve[
          j /. Table[y[j] -> \(c[i - 1]\)[\([j]\)], {j, n + 1}], 
          
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \(-lista\) /. 
            Table[y[j] -> \(c[i - 1]\)[\([j]\)], {j, n + 1}]]; 
      
\ \ \ \ \ \ \ c[i] = c[i - 1] + u[i]; 
      iter = iter + 1]
\)\(Print["Solución:", c[iter]]
\)Print["Número de iteraciones:    ", iter]

Contraste gráfico

graprox = 
      ListPlot[Table[{nodo[[i]], c[iter][[i]]}, {i, n + 1}], 
        PlotStyle -> PointSize[0.015]];

grex = Plot[(E^x), {x, a, b}];

Show[grex, graprox];

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) Resuelvalo 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) Con las datos del apartado anterior, aplique el método de tiro para problemas no lineales al problema anterior y calcule, de nuevo, el error cometido.
d) Resuelva las mismas cuestiones aplicando el método de diferencias finitas.
e) ¿Qué método ofrece mejores resultados? Realice un constraste gráfico de los tres métodos entre sí y de éstos con la solución exacta.

2.- Dado el problema de contorno

		\(TraditionalForm\`\(y'\)' = y'\^2 + y\^2 - sen\ x\  - 1, \ \ 
      x ∈ \((0, π\/2)\), 
y(0) = 0, \ \ \ 
      y(π\/2) = 1.\)

resuélvalo numéricamente mediante dos métodos diferentes, a partir de una partición uniforme del intervalo [0,(π/2)] de paso (π/32). Contraste la solución obtenida con la solución exacta del problema y estime el error cometido.