s3Eva_IIT2007_T1 EDP Eliptica, problema de frontera

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:

δ2uδx2+δ2uδy2=4 \frac{\delta^2u}{\delta x^2} +\frac{\delta^2 u}{\delta y^2} = 4

Que en su forma discreta, con diferencias divididas centradas es:

u[i1,j]2u[i,j]+u[i+1,j](Δx)2+ \frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{(\Delta x)^2} + u[i,j1]2u[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

Al agrupar constantes se convierte en:

u[i1,j]2u[i,j]+u[i+1,j]+ u[i-1,j]-2u[i,j]+u[i+1,j] + +(Δx)2(Δy)2(u[i,j1]2u[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

Siendo Δx= 1/3 y Δy =2/3, se mantendrá la relación λ = (Δx/Δy)2

u[i1,j]2u[i,j]+u[i+1,j]+ u[i-1,j]-2u[i,j]+u[i+1,j] + +λu[i,j1]2λu[i,j]+λu[i,j+1]= + \lambda u[i,j-1]-2\lambda u[i,j]+\lambda u[i,j+1] = =4(Δx)2 = 4(\Delta x)^2
u[i1,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,j1]+λu[i,j+1]=4(Δx)2 + \lambda u[i,j-1]+\lambda u[i,j+1] = 4(\Delta 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[1,0]+λu[1,2]=4(Δx)2 + \lambda u[1,0]+\lambda u[1,2] = 4(\Delta x)^2 Δy22(1+λ)u[1,1]+u[2,1]+ \Delta y^2-2(1+\lambda)u[1,1]+u[2,1] + +Δx2+λu[1,2]=4(Δx)2 + \Delta x^2+\lambda u[1,2] = 4(\Delta 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] =4(Δx)2Δy2Δx2 = 4(\Delta x)^2 - \Delta y^2 - \Delta x^2 2(1+0.25)u[1,1]+u[2,1]+0.25u[1,2] -2(1+0.25)u[1,1]+u[2,1] +0.25 u[1,2] =4(1/3)2(2/3)2(1/3)2 = 4(1/3)^2 - (2/3)^2- (1/3)^2 2.5u[1,1]+u[2,1]+0.25u[1,2]=0.1111 -2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111

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[2,0]+λu[2,2]=4(Δx)2 + \lambda u[2,0]+\lambda u[2,2] = 4(\Delta x)^2 u[1,1]2(1+λ)u[2,1]+(δy1)2+ u[1,1]-2(1+\lambda)u[2,1]+(\delta y-1)^2 + +Δx2+λu[2,2]=4(Δx)2 + \Delta x^2 +\lambda u[2,2] = 4(\Delta 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] =4(Δx)2(Δy1)2Δx2 = 4(\Delta x)^2 - (\Delta y-1)^2 - \Delta x^2 u[1,1]2(1+0.25)u[2,1]+0.25u[2,2] u[1,1]-2(1+0.25)u[2,1] + 0.25 u[2,2] =4(1/3)2(2/31)2(1/3)2 = 4(1/3)^2 - (2/3-1)^2 - (1/3)^2 u[1,1]2.5u[2,1]+0.25u[2,2]=0.2222 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222

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[1,1]+λu[1,3]=4(Δx)2 + \lambda u[1,1]+\lambda u[1,3] = 4(\Delta x)^2 (2Δy2)2(1+λ)u[1,2]+u[2,2]+ (2\Delta y^2)-2(1+\lambda)u[1,2]+u[2,2] + +λu[1,1]+λ(Δx1)2=4(Δx)2 + \lambda u[1,1]+\lambda (\Delta x-1)^2 = 4(\Delta 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] =4(Δx)2(2Δy2)λ(Δx1)2 = 4(\Delta x)^2 - (2\Delta y^2) -\lambda (\Delta x-1)^2 2(1+0.25)u[1,2]+u[2,2]+0.25u[1,1] -2(1+0.25)u[1,2]+u[2,2] + 0.25 u[1,1] =4(1/3)2(2(2/3)2]0.25(1/31)2 = 4(1/3)^2 - (2(2/3)^2] -0.25 (1/3-1)^2 2.5u[1,2]+u[2,2]+0.25u[1,1]=0.5555 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555

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[2,1]+λu[2,3]=4(Δx)2 + \lambda u[2,1]+\lambda u[2,3] = 4(\Delta x)^2 u[1,2]2(1+λ)u[2,2]+(Δy1)2+ u[1,2]-2(1+\lambda)u[2,2]+(\Delta y-1)^2 + +λu[2,1]+λ(Δx1)2=4(Δx)2 + \lambda u[2,1]+\lambda (\Delta x-1)^2 = 4(\Delta 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] =4(Δx)2(Δy1)2λ(Δx1)2 =4(\Delta x)^2- (\Delta y-1)^2 -\lambda (\Delta x-1)^2 u[1,2]2(1+0.25)u[2,2]+0.25u[2,1] u[1,2]-2(1+0.25)u[2,2] + 0.25 u[2,1] =4(1/3)2(2/31)20.25(1/31)2 =4(1/3)^2- (2/3-1)^2 -0.25(1/3-1)^2 u[1,2]2.5u[2,2]+0.25u[2,1]=0.2222 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222

Obtiene el sistema de ecuaciones a resolver:

2.5u[1,1]+u[2,1]+0.25u[1,2]=0.1111 -2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111 u[1,1]2.5u[2,1]+0.25u[2,2]=0.2222 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222 2.5u[1,2]+u[2,2]+0.25u[1,1]=0.5555 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555 u[1,2]2.5u[2,2]+0.25u[2,1]=0.2222 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222

Se escribe en forma matricial Ax=B

[2.510.25012.500.250.2502.5100.2512.5][u[1,1]u[2,1]u[1,2]u[2,2]]=[0.11110.22220.55550.2222]\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}

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]

s3Eva_IT2009_T3 Integrar Simpson compuesta

Ejercicio: 3Eva_IT2009_T3 Integrar Simpson compuesta

La fórmula a integrar es:

01cos(2x)x1/3δx \int_0^1 \frac{\cos (2x)}{x^{1/3}} \delta x

Que tiene la forma:

da como resultado:

el área bajo la curva es:  0.879822622256

Algoritmo en Python

# 3Eva_IT2009_T3 Integrar Simpson compuesta
import numpy as np
import matplotlib.pyplot as plt

funcionx = lambda x: np.cos(2*x)/(x**(1/3))

# INGRESO
a = 0.00001
b = 1
tramos = 10000

# PROCEDIMIENTO
h = (b-a)/tramos
x = a
area = 0
for i in range(0,tramos,2):
    deltaA = (h/3)*(funcionx(x)+4*funcionx(x+h)+funcionx(x+2*h))
    area = area + deltaA
    x = x + 2*h
    
# para la gráfica
xi = np.linspace(a,b,tramos+1)
fxi = funcionx(xi)

print('el área bajo la curva es: ', area)

# Gráfica
plt.plot(xi,fxi)
plt.title('Funcion a integrar')
plt.grid()
plt.xlabel('x')
plt.show()

Tarea: convertir a función el cálculo de Simpson