s3Eva_IIT2019_T2 Diferenciación, valor en frontera

La ecuación de problema de valor de frontera, con h = 1/4 = 0.25:

y'' = -(x+1)y' + 2y + (1-x^2) e^{-x}

0 ≤ x ≤ 1
y(0) = -1
y(1) = 0

Se interpreta en la tabla como los valores que faltan por encontrar:

i 0 1 2 3 4
xi 0 1/4 1/2 3/4 1
yi -1 0

Por ejemplo, se usan las derivadas en diferencias finitas divididas centradas para segunda derivada y hacia adelante para primera derivada.

Semejante al procedimiento aplicado para métodos con EDO.

f''(x_i) = \frac{f(x_{i+1})-2f(x_{i})+f(x_{i-1})}{h^2} + O(h^2) f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} + O(h)

Se formula entonces la ecuación en forma discreta, usando solo los índices para los puntos yi:

\frac{y[i+1]-2y[i]+y[i-1]}{h^2} = -(x_i +1) \frac{y[i+1]-y[i]}{h} + 2y[i] + (1-x_i ^2) e^{-x_i}

se multiplica todo por h2

y[i+1]-2y[i]+y[i-1] = -h(x_i +1)(y[i+1]-y[i]) + + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1]-2y[i]+y[i-1] = -h(x_i +1)y[i+1] + + h(x_i +1)y[i] + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1](1 +h(x_i +1)) + y[i](-2 - h(x_i +1)-2 h^2)+y[i-1] = = h^2 (1-x_i ^2) e^{-x_i}

Ecuación que se aplica en cada uno de los puntos desconocidos con i =1,2,3


i = 1

y[2](1 +h(x_1 +1)) + y[1](-2 - h(x_1 +1)-2 h^2)+y[0] = = h^2 (1-x_1 ^2) e^{-x_1} y[2]\Big(1 +\frac{1}{4}\Big(\frac{1}{4} +1\Big)\Big) + y[1]\Big(-2 - \frac{1}{4}\Big(\frac{1}{4} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big) -1 = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{4}\Big) ^2\Big) e^{-1/4} y[2]\Big(1 +\frac{5}{16} \Big) + y[1]\Big(-2 - \frac{5}{16}- \frac{2}{16}\Big) = = 1+ \frac{1}{16} \Big(1-\frac{1}{16}\Big) e^{-1/4} \frac{21}{16} y[2]- \frac{39}{16}y[1] = 1+ \frac{15}{16^2} e^{-1/4} - \frac{39}{16}y[1] + \frac{21}{16} y[2] = 1+ \frac{15}{16^2} e^{-1/4}

multiplicando ambos lados de la ecuacion por 16 y reordenando

- 39 y[1] + 21 y[2] = 16+ \frac{15}{16} e^{-1/4}

i = 2

y[3](1 +h(x_2 +1)) + y[2](-2 - h(x_2 +1)-2 h^2)+y[1] = = h^2 (1-x_2 ^2) e^{-x_2} y[3]\Big(1 +\frac{1}{4}\Big(\frac{1}{2} +1\Big)\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{1}{2} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big)+y[1] = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{2}\Big) ^2\Big) e^{-\frac{1}{2}} y[3]\Big(1 +\frac{3}{8}\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{3}{2}\Big)- \frac{1}{8}\Big)+y[1] = = \Big(\frac{1}{16}\Big) \Big(1-\frac{1}{4}\Big) e^{-\frac{1}{2}} \frac{11}{8} y[3] - \frac{21}{8} y[2]+y[1] = \frac{1}{16} \Big(\frac{3}{4}\Big) e^{-\frac{1}{2}}

multiplicando ambos lados por 8 y reordenando,

8y[1] - 21 y[2] + 11 y[3] = \frac{3}{8} e^{-\frac{1}{2}}

i = 3

y[4](1 +h(x_3 +1)) + y[3](-2 - h(x_3 +1)-2 h^2)+y[2] = = h^2 (1-x_3 ^2) e^{-x_3} (0) \Big(1 +h\Big(x_3 +1\Big)\Big) + y[3]\Big(-2 - \frac{1}{4}\Big(\frac{3}{4} +1\Big)-2 \frac{1}{16}\Big)+y[2] = = \frac{1}{16}\Big (1-\Big(\frac{3}{4}\Big) ^2 \Big) e^{-3/4} y[3]\Big(-2 - \frac{7}{16} - \frac{2}{16})+y[2] = \frac{1}{16}\Big (1-\frac{9}{16}\Big) e^{-3/4}

multiplicando todo por 16 y reordenando:

16 y[2] - 41 y[3] = \frac{7}{16} e^{-3/4}

Con lo que se puede crear la forma matricial del sistema de ecuaciones Ax=B

\begin{bmatrix} -39 && 21 && 0\\8 && 21 && 11 \\ 0 && 16 && 41 \end{bmatrix}\begin{bmatrix} y[1]\\ y[2] \\y[3] \end{bmatrix} =\begin{bmatrix} 16+ \frac{15}{16} e^{-1/4} \\ \frac{3}{8} e^{-\frac{1}{2}} \\ \frac{7}{16} e^{-3/4} \end{bmatrix}

con lo que resolviendo la matriz se obtienen los valroes para y[1], y[2], y[3]

solución de matriz: 
[-0.59029143 -0.29958287 -0.12195088]

con lo que se completan los puntos de la tabla,

Solución de ecuación
x[i] =  [ 0.   0.25        0.5         0.75        1.  ]
y[i] =  [-1.  -0.59029143 -0.29958287 -0.12195088  0.  ]

con la siguiente gráfica:


Algoritmo para solución de matriz con Python

tarea: modificar para cambiar el valor del tamaño de paso.

# Problema de Frontera
import numpy as np
import matplotlib.pyplot as plt

# INGRESO

h = 1/4
y0 = -1
y1 = 0

xi = np.arange(0,1+h,h)
n = len(xi)
yi = np.zeros(n,dtype=float)
yi[0] = y0
yi[n-1] = y1
    
A = np.array([[-39.0, 21,  0],
              [  8.0, -21, 11],
              [  0.0, 16, -41]])
B = np.array([16+(15/16)*np.exp(-1/4),
              (3/8)*np.exp(-1/2),
              (7/16)*np.exp(-3/4)])

# PROCEDIMIENTO
x = np.linalg.solve(A,B)

for i in range(1,n-1,1):
    yi[i] = x[i-1]

# SALIDA
print('solución de matriz: ')
print(x)
print('Solución de ecuación')
print('x[i] = ',xi)
print('y[i] = ',yi)

# Grafica
plt.plot(xi,yi,'ro')
plt.plot(xi,yi)
plt.show()