Clase I

1. Aproximaciones de derivadas por diferencias finitas

Trataremos al problema de manera unidimensional para facilitar el entendimiento. Sea \(x_0\) un número real cualquiera y \(h\) un número positivo. Se defien una malla de paso \(h\) asociado a \(x_0\) como el conjunto de puntos.

\[x_i=x_0+ih, \hspace{1cm} i=1,2,3,...,N\]

Se aproxima una función \(y=y(x)\), solución exacta de la ecuación diferencial en estudio, asi como sus derivadas en los puntos de la malla. La herramienta básica en el cálculo de aproximaciones para las derivadas es la serie de Taylor que relaciona valores de la función y sus derivadas, en el punto \(x\), con los valores de esa misma función en una vecindad de \(x\), es decir, con \(y(x+h)\). Si \(y(x)\) tiene derivadas hasta el orden \(n+1\) en \(x\), podemos escribirla como:

\[\begin{equation} y(x+h)=y(x) + hy'(x)h + \frac{h^2}{2!}y''(x) + ... + \frac{h^n}{n!}y^{(n)}(x) + \frac{h^{n+1}}{(n+1)!}y^{(n+1)}(\xi) \hspace{1cm} x<\xi<x+h \hspace{1cm} (1) \end{equation}\]

El útlimo término de la expresión \((1)\) representa el error de la aproximación de \(y(x+h)\) por el polinomio evaluado en la variable h de grado \(n\).

\[P_n(h)=y(x)+hy'(x)+\dfrac{h^2}{2!}y''(x)+...+\dfrac{h^n}{n!}y^{(n)}(x)\]

Si el valor de \(n=1\) se obtiene una aproximación para la derivada \(y'(x)\), conocida como la fórmula progresiva

\[ y(x+h)=y(x)+hy'(x)+\frac{h^2}{2!}y''(\xi)\\ y(x+h)=y(x)+hy'(x)+\frac{h^2}{2}y''(\xi)\\ y(x+h)-y(x)-\frac{h^2}{2}y''(\xi)=hy'(x)\\ y'(x)=\frac{y(x+h)-y(x)-\frac{h^2}{2}y''(\xi)}{h}\\ y'(x)=\frac{y(x+h)-y(x)}{h}-\frac{h}{2}y''(\xi) \] Luego utilizamos la siguiente notación para representar la diferencia progresiva de \(y(x)\)

\[ \Delta y(x)=y(x+h)-y(x) \] El término \(-\frac{h}{2}y''(\xi)\) representa el error de esa aproximación. Por tanto se tiene

\[ y'(x)=\frac{1}{h}\Delta y(x) - \frac{h}{2}y''(\xi) \]

De modo semejante sustituimos \(h\) por \(-h\) en la ecuación \((1)\), aún con \(n=1\), obtenemos la fórmula regresiva

\[ y(x-h)=y(x)-hy'(x)+\frac{h^2}{2!}y''(\xi)\\ y(x-h)=y(x)-hy'(x)+\frac{h^2}{2}y''(\xi)\\ hy'(x)=y(x)-y(x-h)+\frac{h^2}{2}y''(\xi)\\ y'(x)=\frac{y(x)-y(x-h)+\frac{h^2}{2}y''(\xi)}{h}\\ y'(x)=\frac{y(x)-y(x-h)}{h}+\frac{h}{2}y''(\xi)\\ \] Luego utilizamos la siguiente notación para representar la diferencia regresiva de \(y(x)\)

\[ \Delta y(x)=y(x)-y(x-h) \] El término \(\frac{h}{2}y''(\xi)\) representa el error de esa aproximación. Por tanto se tiene

\[ y'(x)=\frac{1}{h}\Delta y(x)+\frac{h}{2}y''(\xi) \]

Luego si tomamos para \(n=2\) en la ecuación \((1)\) para \(h\) y \(-h\), respectivamente tenemos

\[ y(x+h)=y(x)+hy'(x)+\frac{h^2}{2!}y''(x)+\frac{h^3}{3!}y'''(\xi_1) \]

\[ y(x-h)=y(x)-hy'(x)+\frac{h^2}{2!}y''(x)-\frac{h^3}{3!}y'''(\xi_2) \]

ahora realizamos la diferencia de \(y(x+h) - y(x-h)\) y se tiene

\[ y(x+h)-y(x-h)=2hy'(x)+\frac{h^3}{3!}y'''(\xi)\\ 2hy'(x)=y(x+h)-y(x-h)-\frac{h^3}{6}y'''(\xi)\\ y'(x)=\frac{y(x+h)-y(x-h)-\frac{h^3}{6}y'''(\xi)}{2h} \] la fórmula central o fórmula de diferencia centrada

\[ y'(x)=\frac{y(x+h)-y(x-h)}{2h}-\frac{h^2}{12}y'''(\xi) \hspace{1cm} (2) \] donde \(\xi \in (x-h,x+h)\) y fue utilizado el teorema de valor intermedio valido para funciones continua:

\[ \frac{1}{2}(y'''(\xi_1)+y'''(\xi_2))=y'''(\xi) \] Para algún \(\xi \in [mín\{\xi_1,xi_2\}, máx\{\xi_1,\xi_2\}]\). Denotando la diferencia central por \(\delta_h y(x)\) entonces de \((2)\) tenemos

\[ y'(x)=\frac{1}{2h}\delta_h y(x)-\frac{h^2}{12}y'''(\xi) \]

2. Error y orden de aproximación de una fórmula de diferencia

Sea \(F(x,h)\) una fórmula de diferencia para la aproximacción de la derivada de orden \(q\) de una función \(y(x)\) con error \(E(x,h)\), entonces:

\[y^{(q)}(x)=F(x,h)+E(x,h)\]

Decimos que la fórmula \(F(x,h)\) es de orden \(p\) si \(E(x,h)=h^pR(x)\), donde \(R(x)\) no depende de \(h\). En ese caso usamos la notación \(E(x,h)=O(h^p)\). Esa notación significa que \(\lim\limits_{h \rightarrow 0} \dfrac{E(x,h)}{h^p}\) es una constante finita.

Por ejemplo, en el caso de la fórmula centrada tenemos que:

\[y'(x)=F(x,h)=\dfrac{y(x+h)-y(x-h)}{2h} \hspace{1cm} \text{y} \hspace{1cm} E(x,h)=-\dfrac{h^2}{12}y'''(\xi) \hspace{1cm} (3)\] Entonces se tiene que la fórmula es de segundo orden, \(E(x,h)=O(h^2)\)

Ahora hallamos la aproximación para las segunda derivada. Haciendo \(n=3\) en \((1)\).

\[ y(x+h)=y(x)+hy'(x)+\frac{h^2}{2!}y''(x)+\frac{h^3}{3!}y'''(x)+\frac{h^4}{4!}y^{(4)}(\xi_1) \]

\[ y(x-h)=y(x)-hy'(x)+\frac{h^2}{2!}y''(x)-\frac{h^3}{3!}y'''(x)+\frac{h^4}{4!}y^{(4)}(\xi_2) \]

ahora realizamos la suma de \(y(x+h) + y(x-h)\) y se tiene

\[ y(x+h)+y(x-h)=2y(x)+h^2y''(x)+\frac{h^4}{12}y^{(4)}(\xi)\\ h^2y''(x) = y(x+h)+y(x-h)-2y(x)-\frac{h^4}{12}y^{(4)}(\xi)\\ y''(x)=\frac{y(x+h)+y(x-h)-2y(x)-\frac{h^4}{12}y^{(4)}(\xi)}{h^2}\\ y''(x)=\frac{y(x+h)-2y(x)+y(x-h)}{h^2}-\frac{h^2}{12}y^{(4)}(\xi)\\ \text{Denotando al diferencia centrada por} \hspace{1cm} \delta_{\frac{h}{2}}y(x)\\ y''(x)=\frac{1}{h^2}\delta_{\frac{h}{2}}y(x)-\frac{h^2}{12}y^{(4)}(\xi) \] donde \(\xi \in (x+h,x+h)\). En este caso

\[ F(x,h)=\frac{y(x+h)-2y(x)+y(x-h)}{h^2} \hspace{1cm} \text{y } \hspace{1cm} E(x,h)=-\frac{h^2}{12}y^{(4)}(\xi) \hspace{1cm} (4) \]

Entonces se tiene que la fórmula es de orden \(2\).

3. Problema de valor incial en ecuaciones ordinarias

Un problema de valor inicial (PVI) es un problema de evolución, en el cual la información inicial (conocida), es propagada para el interior del dominio por la ecuacion diferencial. Expresamos matemáticamente el PVI en la forma:

\[ y' = f(x,y) \hspace{4cm} \\ y(a)=\alpha \hspace{4.8cm} (5) \] Donde \(f:ℝ^2 \longrightarrow ℝ^n\) es una función continua. La función \(y=y(x) \hspace{0.5cm} (x \geq a)\) es la solución y \(\alpha\) es el valor inicial en \(a\). Ecuaciones ordinarias de mayor orden pueden ser transformadas a sistemas de ecuaciones de la forma \((5)\). Por ejemplo, la ecuación de segundo orden,

\[ y'' = f(x,y,y') \hspace{4.15cm} \\ y(a)=\alpha \hspace{6cm} \\ y'(a) = \beta \hspace{5.7cm} (6) \] con solución \(y(x)\) y los valores iniciales \(\alpha\) y \(\beta\) dados en el punto a, puede ser reescrita en forma de un sistema

\[ y' = v, \hspace{2.3cm} y(a) = \alpha \\ v' = f(x, y, v), \hspace{1cm} v(a) = \beta \] Para mostrar la existencia y unicidad supondremos que la función \(f(x, y)\) satisfacen las condiciones:

  1. \(f : E \longrightarrow ℝ^n\) es continua, donde \(E =\left\lbrace{(x, y), x \in [a, b], y \in ℝ^n }\right\rbrace.\)

  2. existe una constante \(L\) tal que para todo \(x \in [a, b]\) y cualesquiera dos vectores \(y\) y \(y^*\) en \(ℝ^n\) \[ |f(x, y) − f(x, y^∗)| \leq L |y − y^∗| \hspace{2cm} (7) \]

4. Ejemplos de métodos numéricos

Euler en 1768 usando diferencias progresivas, desarrollo una aproximación para \((5)\), método conocido como el Método de Euler.

4.1 Método de Euler

El intervalo \([a,b]\) es dividido en \(N\) partes iguales, cada uno de longitud \(h\) formando un conjunto discreto, con \(x_0=a\) y \(x_N=b\), \(R_h=\left\lbrace{x_i=a+ih,i=0,1,2,...,N}\right\rbrace\). Aquí \(h=(b-a)/N\). Sea \(y_i \approx y(x_i)\) aproximaciones para \(y(x_i), i=1,2,...,N\) y \(y_0=\alpha\). En cada uno de los puntos \(x_0,x_1,...,N-1\), aproximamos la ecuación diferencial \((5)\) por

\[ \Delta y_i = f(x_i,y_i), \hspace{1cm} i=0,1,2,...,N-1 \] \[ y_{i+1}=y_i+hf(x_i,y_i), \hspace{1cm} i=0,1,2,...,N-1 \] Este es el Método de Euler que es un método explicito. Observe que el valor de \(y_{i+1}\) es calculado explicitamente en función de \(y_i\), mismo cuando la función es no lineal. Si en lugar de \(\Delta\) fuese usado \(\nabla\) tendríamos la versión implícita del Método de Euler.

\[ \nabla y_i = f(x_i,y_i), \hspace{1cm} i=1,2,...,N \] \[ y_i=y_{i-1}+hf(x_i,y_i), \hspace{1cm} i=1,2,...,N \] \[ y_{i+1}=y_i+hf(x_{i+1},y_{i+1}), \hspace{1cm} i=0,1,2,...,N-1 \hspace{1cm} (8) \] Note que en el caso en que f es no lineal, la ecuación \((8)\) define el vector \(y_{i+1}\) implícitamente y por tanto su calculo exige la utilización métodos para la solución de sistemas de ecuaciones algebraicas no lineales.

Observe también que la ecuación diferencial en el proceso de discretización, fue substituida por otra ecuación llamada ecuaciones en diferencias que tiene su propia solución independiente de la solución de aquella ecuación diferencial.

Ejemplo

Considere el PVI \[ y'=-y+x \] \[ y(0)=1 \hspace{1cm} (9) \] Cuya solución excata es \(y(x)=x-1+2e^{-x}\). Encuentre las aproximaciones del problema \((9)\) usando los Métodos de Euler explícito, Euler implícito y de los Trapecios cuando \(h=0.1\) y \(h=0.2\). Además grafique las aproximaciones junto con la solución exacta para cada \(h\).

Método Euler explícito Calculamos \(N\), para ello usamos la fórmula \(h=(b-a)/N\)

h = 0.1
a = 0
b = 1
N = (b-a)/h
sprintf("El intervalo es dividido en %s partes iguales.",N)
## [1] "El intervalo es dividido en 10 partes iguales."
Rh<-c()
for(i in 0:(N)){
  xi = a + i*h
  Rh<-c(Rh,xi)
}
print(Rh)
##  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

El método explícito es el siguiente \(y_{i+1}=y_i+h(-y_i+x_i)\) Se tiene como dato que \(y_0=\alpha\) entonces \(y_0=1\)

y<-c(1)
for(i in 0:(N)){
  yi = y[i]+h*(-y[i]+Rh[i])
  y<-c(y,yi)
}
print(y)
##  [1] 1.0000000 0.9000000 0.8200000 0.7580000 0.7122000 0.6809800 0.6628820
##  [8] 0.6565938 0.6609344 0.6748410 0.6973569

Calculamos la solución exacta, para ello definimos la función.

funcion<-function(x){
  return(y = x-1+2*exp(-x))
}

sol_exacta<-c()

for(i in 0:(N+1)){
  sol = funcion(Rh[i])
  sol_exacta<-c(sol_exacta,sol)
}
print(sol_exacta)
##  [1] 1.0000000 0.9096748 0.8374615 0.7816364 0.7406401 0.7130613 0.6976233
##  [8] 0.6931706 0.6986579 0.7131393 0.7357589

Realizamos las graficas de la aproximacion y solución exacta cuando \(h=0.1\)

plot(Rh,y, type="b", cex=2, pch=21, bg="blue",col="black", main="Solución exacta y aproximada cuando h = 0.1")
lines(Rh,sol_exacta, type="b", cex=2, pch=21, bg="green",col="red")
legend("topright",legend=c("Aproximada","Exacta"),fill = c("blue", "green"), lty = 1, col = 1:2)