Desarrollo analítico del método explícito de una ecuación diferencial parcial con Sympy
Ejemplos:
\frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}El algoritmo desarrolla la parte analítica de una EDP modelo con el siguiente resultado:
EDP=0 :
2
d d
- K*--(u(x, t)) + ---(u(x, t)) = 0
dt 2
dx
Lambda L :
dt
-----
2
K*dx
discreta = 0 :
K*(-U(i, j) + U(i, j + 1)) -2*U(i, j) + U(i - 1, j) + U(i + 1, j)
- -------------------------- + --------------------------------------
dt 2
dx
discreta_L = 0 :
L*U(i - 1, j) + L*U(i + 1, j) + (1 - 2*L)*U(i, j) - U(i, j + 1) = 0
U(i, j + 1) :
L*U(i - 1, j) + L*U(i + 1, j) + (1 - 2*L)*U(i, j)
Instrucciones en Phyton
# Ecuaciones Diferenciales Parciales Parabólicas con Sympy import sympy as sym # variables continuas t = sym.Symbol('t',real=True) x = sym.Symbol('x',real=True) u = sym.Function('u')(x,t) K = sym.Symbol('K',real=True) # variables discretas i = sym.Symbol('i',integer=True,positive=True) j = sym.Symbol('j',integer=True,positive=True) dx = sym.Symbol('dx',real=True,positive=True) dt = sym.Symbol('dt',real=True,positive=True) L = sym.Symbol('L',real=True) U = sym.Function('U')(i,j) # INGRESO # ecuacion: LHS=RHS LHS = sym.diff(u,x,2) RHS = K*sym.diff(u,t,1) ecuacion = sym.Eq(LHS-RHS,0) dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(dx**2), sym.diff(u,t,1): (U.subs(j,j+1)-U)/dt} buscar = U.subs(j,j+1) # PROCEDIMIENTO desarrollo={} desarrollo['EDP=0'] = ecuacion lamb = dt/K/(dx**2) desarrollo['Lambda L'] = lamb discreta = ecuacion.lhs # forma discreta EDP for factor_k in dif_dividida: discreta = discreta.subs(factor_k,dif_dividida[factor_k]) desarrollo['discreta = 0'] = discreta # simplifica con Lambda L discreta_L = sym.expand(discreta*dt/K) discreta_L = discreta_L.subs(dt/K/(dx**2),L) discreta_L = sym.collect(discreta_L,U) discreta_L = sym.Eq(discreta_L,0) desarrollo['discreta_L = 0'] = discreta_L # busca el término no conocido u_explicito = sym.solve(discreta_L,buscar) u_explicito = sym.collect(u_explicito[0],U) desarrollo[buscar] = u_explicito # SALIDA for k in desarrollo: print('\n',k,' :') sym.pprint(desarrollo[k])