Ejercicio : 3Eva_IIT2007_T1 EDP Eliptica, problema de frontera
Con los datos del ejercicio se plantean de la siguiente forma en los bordes:
La ecuación a resolver es:
δ 2 u δ x 2 + δ 2 u δ y 2 = 4 \frac{\delta^2u}{\delta x^2} +\frac{\delta^2 u}{\delta y^2} = 4 δ x 2 δ 2 u + δ y 2 δ 2 u = 4
Que en su forma discreta, con diferencias divididas centradas es:
u [ i − 1 , j ] − 2 u [ i , j ] + u [ i + 1 , j ] ( Δ x ) 2 + \frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{(\Delta x)^2} + ( Δ x ) 2 u [ i − 1 , j ] − 2 u [ i , j ] + u [ i + 1 , j ] +
u [ i , j − 1 ] − 2 u [ i , j ] + u [ i , j + 1 ] ( Δ y ) 2 = 4 \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{(\Delta y)^2} = 4 ( Δ y ) 2 u [ i , j − 1 ] − 2 u [ i , j ] + u [ i , j + 1 ] = 4
Al agrupar constantes se convierte en:
u [ i − 1 , j ] − 2 u [ i , j ] + u [ i + 1 , j ] + u[i-1,j]-2u[i,j]+u[i+1,j] + u [ i − 1 , j ] − 2 u [ i , j ] + u [ i + 1 , j ] +
+ ( Δ x ) 2 ( Δ y ) 2 ( u [ i , j − 1 ] − 2 u [ i , j ] + u [ i , j + 1 ] ) = 4 ( Δ x ) 2 + \frac{(\Delta x)^2}{(\Delta y)^2} \Big(u[i,j-1]-2u[i,j]+u[i,j+1] \Big)= 4(\Delta x)^2 + ( Δ y ) 2 ( Δ x ) 2 ( u [ i , j − 1 ] − 2 u [ i , j ] + u [ i , j + 1 ] ) = 4 ( Δ x ) 2
Siendo Δx= 1/3 y Δy =2/3, se mantendrá la relación λ = (Δx/Δy)2
u [ i − 1 , j ] − 2 u [ i , j ] + u [ i + 1 , j ] + u[i-1,j]-2u[i,j]+u[i+1,j] + u [ i − 1 , j ] − 2 u [ i , j ] + u [ i + 1 , j ] +
+ λ u [ i , j − 1 ] − 2 λ u [ i , j ] + λ u [ i , j + 1 ] = + \lambda u[i,j-1]-2\lambda u[i,j]+\lambda u[i,j+1] = + λ u [ i , j − 1 ] − 2 λ u [ i , j ] + λ u [ i , j + 1 ] =
= 4 ( Δ x ) 2 = 4(\Delta x)^2 = 4 ( Δ x ) 2
u [ i − 1 , j ] − 2 ( 1 + λ ) u [ i , j ] + u [ i + 1 , j ] + u[i-1,j]-2(1+\lambda)u[i,j]+u[i+1,j] + u [ i − 1 , j ] − 2 ( 1 + λ ) u [ i , j ] + u [ i + 1 , j ] +
+ λ u [ i , j − 1 ] + λ u [ i , j + 1 ] = 4 ( Δ x ) 2 + \lambda u[i,j-1]+\lambda u[i,j+1] = 4(\Delta x)^2 + λ u [ i , j − 1 ] + λ u [ i , j + 1 ] = 4 ( Δ x ) 2
Se pueden realizar las iteraciones para los nodos y reemplaza los valores de frontera:
i=1 y j=1
u [ 0 , 1 ] − 2 ( 1 + λ ) u [ 1 , 1 ] + u [ 2 , 1 ] + u[0,1]-2(1+\lambda)u[1,1]+u[2,1] + u [ 0 , 1 ] − 2 ( 1 + λ ) u [ 1 , 1 ] + u [ 2 , 1 ] +
+ λ u [ 1 , 0 ] + λ u [ 1 , 2 ] = 4 ( Δ x ) 2 + \lambda u[1,0]+\lambda u[1,2] = 4(\Delta x)^2 + λ u [ 1 , 0 ] + λ u [ 1 , 2 ] = 4 ( Δ x ) 2
Δ y 2 − 2 ( 1 + λ ) u [ 1 , 1 ] + u [ 2 , 1 ] + \Delta y^2-2(1+\lambda)u[1,1]+u[2,1] + Δ y 2 − 2 ( 1 + λ ) u [ 1 , 1 ] + u [ 2 , 1 ] +
+ Δ x 2 + λ u [ 1 , 2 ] = 4 ( Δ x ) 2 + \Delta x^2+\lambda u[1,2] = 4(\Delta x)^2 + Δ x 2 + λ u [ 1 , 2 ] = 4 ( Δ x ) 2
− 2 ( 1 + λ ) u [ 1 , 1 ] + u [ 2 , 1 ] + λ u [ 1 , 2 ] -2(1+\lambda)u[1,1]+u[2,1] +\lambda u[1,2] − 2 ( 1 + λ ) u [ 1 , 1 ] + u [ 2 , 1 ] + λ u [ 1 , 2 ]
= 4 ( Δ x ) 2 − Δ y 2 − Δ x 2 = 4(\Delta x)^2 - \Delta y^2 - \Delta x^2 = 4 ( Δ x ) 2 − Δ y 2 − Δ x 2
− 2 ( 1 + 0 . 2 5 ) u [ 1 , 1 ] + u [ 2 , 1 ] + 0 . 2 5 u [ 1 , 2 ] -2(1+0.25)u[1,1]+u[2,1] +0.25 u[1,2] − 2 ( 1 + 0 . 2 5 ) u [ 1 , 1 ] + u [ 2 , 1 ] + 0 . 2 5 u [ 1 , 2 ]
= 4 ( 1 / 3 ) 2 − ( 2 / 3 ) 2 − ( 1 / 3 ) 2 = 4(1/3)^2 - (2/3)^2- (1/3)^2 = 4 ( 1 / 3 ) 2 − ( 2 / 3 ) 2 − ( 1 / 3 ) 2
− 2 . 5 u [ 1 , 1 ] + u [ 2 , 1 ] + 0 . 2 5 u [ 1 , 2 ] = − 0 . 1 1 1 1 -2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111 − 2 . 5 u [ 1 , 1 ] + u [ 2 , 1 ] + 0 . 2 5 u [ 1 , 2 ] = − 0 . 1 1 1 1
i=2 , j=1
u [ 1 , 1 ] − 2 ( 1 + λ ) u [ 2 , 1 ] + u [ 3 , 1 ] + u[1,1]-2(1+\lambda)u[2,1]+u[3,1] + u [ 1 , 1 ] − 2 ( 1 + λ ) u [ 2 , 1 ] + u [ 3 , 1 ] +
+ λ u [ 2 , 0 ] + λ u [ 2 , 2 ] = 4 ( Δ x ) 2 + \lambda u[2,0]+\lambda u[2,2] = 4(\Delta x)^2 + λ u [ 2 , 0 ] + λ u [ 2 , 2 ] = 4 ( Δ x ) 2
u [ 1 , 1 ] − 2 ( 1 + λ ) u [ 2 , 1 ] + ( δ y − 1 ) 2 + u[1,1]-2(1+\lambda)u[2,1]+(\delta y-1)^2 + u [ 1 , 1 ] − 2 ( 1 + λ ) u [ 2 , 1 ] + ( δ y − 1 ) 2 +
+ Δ x 2 + λ u [ 2 , 2 ] = 4 ( Δ x ) 2 + \Delta x^2 +\lambda u[2,2] = 4(\Delta x)^2 + Δ x 2 + λ u [ 2 , 2 ] = 4 ( Δ x ) 2
u [ 1 , 1 ] − 2 ( 1 + λ ) u [ 2 , 1 ] + λ u [ 2 , 2 ] u[1,1]-2(1+\lambda)u[2,1] + \lambda u[2,2] u [ 1 , 1 ] − 2 ( 1 + λ ) u [ 2 , 1 ] + λ u [ 2 , 2 ]
= 4 ( Δ x ) 2 − ( Δ y − 1 ) 2 − Δ x 2 = 4(\Delta x)^2 - (\Delta y-1)^2 - \Delta x^2 = 4 ( Δ x ) 2 − ( Δ y − 1 ) 2 − Δ x 2
u [ 1 , 1 ] − 2 ( 1 + 0 . 2 5 ) u [ 2 , 1 ] + 0 . 2 5 u [ 2 , 2 ] u[1,1]-2(1+0.25)u[2,1] + 0.25 u[2,2] u [ 1 , 1 ] − 2 ( 1 + 0 . 2 5 ) u [ 2 , 1 ] + 0 . 2 5 u [ 2 , 2 ]
= 4 ( 1 / 3 ) 2 − ( 2 / 3 − 1 ) 2 − ( 1 / 3 ) 2 = 4(1/3)^2 - (2/3-1)^2 - (1/3)^2 = 4 ( 1 / 3 ) 2 − ( 2 / 3 − 1 ) 2 − ( 1 / 3 ) 2
u [ 1 , 1 ] − 2 . 5 u [ 2 , 1 ] + 0 . 2 5 u [ 2 , 2 ] = 0 . 2 2 2 2 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222 u [ 1 , 1 ] − 2 . 5 u [ 2 , 1 ] + 0 . 2 5 u [ 2 , 2 ] = 0 . 2 2 2 2
i=1 , j=2
u [ 0 , 2 ] − 2 ( 1 + λ ) u [ 1 , 2 ] + u [ 2 , 2 ] + u[0,2]-2(1+\lambda)u[1,2]+u[2,2] + u [ 0 , 2 ] − 2 ( 1 + λ ) u [ 1 , 2 ] + u [ 2 , 2 ] +
+ λ u [ 1 , 1 ] + λ u [ 1 , 3 ] = 4 ( Δ x ) 2 + \lambda u[1,1]+\lambda u[1,3] = 4(\Delta x)^2 + λ u [ 1 , 1 ] + λ u [ 1 , 3 ] = 4 ( Δ x ) 2
( 2 Δ y 2 ) − 2 ( 1 + λ ) u [ 1 , 2 ] + u [ 2 , 2 ] + (2\Delta y^2)-2(1+\lambda)u[1,2]+u[2,2] + ( 2 Δ y 2 ) − 2 ( 1 + λ ) u [ 1 , 2 ] + u [ 2 , 2 ] +
+ λ u [ 1 , 1 ] + λ ( Δ x − 1 ) 2 = 4 ( Δ x ) 2 + \lambda u[1,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2 + λ u [ 1 , 1 ] + λ ( Δ x − 1 ) 2 = 4 ( Δ x ) 2
− 2 ( 1 + λ ) u [ 1 , 2 ] + u [ 2 , 2 ] + λ u [ 1 , 1 ] -2(1+\lambda)u[1,2]+u[2,2] + \lambda u[1,1] − 2 ( 1 + λ ) u [ 1 , 2 ] + u [ 2 , 2 ] + λ u [ 1 , 1 ]
= 4 ( Δ x ) 2 − ( 2 Δ y 2 ) − λ ( Δ x − 1 ) 2 = 4(\Delta x)^2 - (2\Delta y^2) -\lambda (\Delta x-1)^2 = 4 ( Δ x ) 2 − ( 2 Δ y 2 ) − λ ( Δ x − 1 ) 2
− 2 ( 1 + 0 . 2 5 ) u [ 1 , 2 ] + u [ 2 , 2 ] + 0 . 2 5 u [ 1 , 1 ] -2(1+0.25)u[1,2]+u[2,2] + 0.25 u[1,1] − 2 ( 1 + 0 . 2 5 ) u [ 1 , 2 ] + u [ 2 , 2 ] + 0 . 2 5 u [ 1 , 1 ]
= 4 ( 1 / 3 ) 2 − ( 2 ( 2 / 3 ) 2 ] − 0 . 2 5 ( 1 / 3 − 1 ) 2 = 4(1/3)^2 - (2(2/3)^2] -0.25 (1/3-1)^2 = 4 ( 1 / 3 ) 2 − ( 2 ( 2 / 3 ) 2 ] − 0 . 2 5 ( 1 / 3 − 1 ) 2
− 2 . 5 u [ 1 , 2 ] + u [ 2 , 2 ] + 0 . 2 5 u [ 1 , 1 ] = − 0 . 5 5 5 5 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555 − 2 . 5 u [ 1 , 2 ] + u [ 2 , 2 ] + 0 . 2 5 u [ 1 , 1 ] = − 0 . 5 5 5 5
i=2, j=2
u [ 1 , 2 ] − 2 ( 1 + λ ) u [ 2 , 2 ] + u [ 3 , 2 ] + u[1,2]-2(1+\lambda)u[2,2]+u[3,2] + u [ 1 , 2 ] − 2 ( 1 + λ ) u [ 2 , 2 ] + u [ 3 , 2 ] +
+ λ u [ 2 , 1 ] + λ u [ 2 , 3 ] = 4 ( Δ x ) 2 + \lambda u[2,1]+\lambda u[2,3] = 4(\Delta x)^2 + λ u [ 2 , 1 ] + λ u [ 2 , 3 ] = 4 ( Δ x ) 2
u [ 1 , 2 ] − 2 ( 1 + λ ) u [ 2 , 2 ] + ( Δ y − 1 ) 2 + u[1,2]-2(1+\lambda)u[2,2]+(\Delta y-1)^2 + u [ 1 , 2 ] − 2 ( 1 + λ ) u [ 2 , 2 ] + ( Δ y − 1 ) 2 +
+ λ u [ 2 , 1 ] + λ ( Δ x − 1 ) 2 = 4 ( Δ x ) 2 + \lambda u[2,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2 + λ u [ 2 , 1 ] + λ ( Δ x − 1 ) 2 = 4 ( Δ x ) 2
u [ 1 , 2 ] − 2 ( 1 + λ ) u [ 2 , 2 ] + λ u [ 2 , 1 ] u[1,2]-2(1+\lambda)u[2,2] + \lambda u[2,1] u [ 1 , 2 ] − 2 ( 1 + λ ) u [ 2 , 2 ] + λ u [ 2 , 1 ]
= 4 ( Δ x ) 2 − ( Δ y − 1 ) 2 − λ ( Δ x − 1 ) 2 =4(\Delta x)^2- (\Delta y-1)^2 -\lambda (\Delta x-1)^2 = 4 ( Δ x ) 2 − ( Δ y − 1 ) 2 − λ ( Δ x − 1 ) 2
u [ 1 , 2 ] − 2 ( 1 + 0 . 2 5 ) u [ 2 , 2 ] + 0 . 2 5 u [ 2 , 1 ] u[1,2]-2(1+0.25)u[2,2] + 0.25 u[2,1] u [ 1 , 2 ] − 2 ( 1 + 0 . 2 5 ) u [ 2 , 2 ] + 0 . 2 5 u [ 2 , 1 ]
= 4 ( 1 / 3 ) 2 − ( 2 / 3 − 1 ) 2 − 0 . 2 5 ( 1 / 3 − 1 ) 2 =4(1/3)^2- (2/3-1)^2 -0.25(1/3-1)^2 = 4 ( 1 / 3 ) 2 − ( 2 / 3 − 1 ) 2 − 0 . 2 5 ( 1 / 3 − 1 ) 2
u [ 1 , 2 ] − 2 . 5 u [ 2 , 2 ] + 0 . 2 5 u [ 2 , 1 ] = − 0 . 2 2 2 2 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222 u [ 1 , 2 ] − 2 . 5 u [ 2 , 2 ] + 0 . 2 5 u [ 2 , 1 ] = − 0 . 2 2 2 2
Obtiene el sistema de ecuaciones a resolver:
− 2 . 5 u [ 1 , 1 ] + u [ 2 , 1 ] + 0 . 2 5 u [ 1 , 2 ] = − 0 . 1 1 1 1 -2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111 − 2 . 5 u [ 1 , 1 ] + u [ 2 , 1 ] + 0 . 2 5 u [ 1 , 2 ] = − 0 . 1 1 1 1
u [ 1 , 1 ] − 2 . 5 u [ 2 , 1 ] + 0 . 2 5 u [ 2 , 2 ] = 0 . 2 2 2 2 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222 u [ 1 , 1 ] − 2 . 5 u [ 2 , 1 ] + 0 . 2 5 u [ 2 , 2 ] = 0 . 2 2 2 2
− 2 . 5 u [ 1 , 2 ] + u [ 2 , 2 ] + 0 . 2 5 u [ 1 , 1 ] = − 0 . 5 5 5 5 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555 − 2 . 5 u [ 1 , 2 ] + u [ 2 , 2 ] + 0 . 2 5 u [ 1 , 1 ] = − 0 . 5 5 5 5
u [ 1 , 2 ] − 2 . 5 u [ 2 , 2 ] + 0 . 2 5 u [ 2 , 1 ] = − 0 . 2 2 2 2 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222 u [ 1 , 2 ] − 2 . 5 u [ 2 , 2 ] + 0 . 2 5 u [ 2 , 1 ] = − 0 . 2 2 2 2
Se escribe en forma matricial Ax=B
[ − 2 . 5 1 0 . 2 5 0 1 − 2 . 5 0 0 . 2 5 0 . 2 5 0 − 2 . 5 1 0 0 . 2 5 1 − 2 . 5 ] [ u [ 1 , 1 ] u [ 2 , 1 ] u [ 1 , 2 ] u [ 2 , 2 ] ] = [ − 0 . 1 1 1 1 0 . 2 2 2 2 − 0 . 5 5 5 5 − 0 . 2 2 2 2 ] \begin {bmatrix} -2.5 && 1&&0.25&&0 \\1&&-2.5&&0&&0.25\\0.25&&0&&-2.5&&1\\0&&0.25&&1&&-2.5\end{bmatrix} \begin {bmatrix} u[1,1] \\ u[2,1]\\ u[1,2]\\u[2,2] \end{bmatrix} = \begin {bmatrix}-0.1111\\0.2222\\-0.5555\\-0.2222 \end{bmatrix} ⎣ ⎢ ⎢ ⎡ − 2 . 5 1 0 . 2 5 0 1 − 2 . 5 0 0 . 2 5 0 . 2 5 0 − 2 . 5 1 0 0 . 2 5 1 − 2 . 5 ⎦ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎡ u [ 1 , 1 ] u [ 2 , 1 ] u [ 1 , 2 ] u [ 2 , 2 ] ⎦ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎡ − 0 . 1 1 1 1 0 . 2 2 2 2 − 0 . 5 5 5 5 − 0 . 2 2 2 2 ⎦ ⎥ ⎥ ⎤
que se resuelve con un algoritmo:
import numpy as np
A = np.array([[-2.5,1,0.25,0],
[1,-2.5,0,0.25],
[0.25,0,-2.5,1],
[0,0.25,1,-2.5]])
B = np.array([-0.1111,0.2222,-0.5555,-0.2222])
x = np.linalg.solve(A,B)
print(x)
y los resultados son:
[ 0.05762549 -0.04492835 0.31156835 0.20901451]