I have an ordinary differential equation like this:
DiffEq = Eq(-ℏ*ℏ*diff(Ψ,x,2)/(2*m) + m*w*w*(x*x)*Ψ/2 - E*Ψ , 0)
I want to perform a variable change :
sp.Eq(u , x*sqrt(m*w/ℏ)) sp.Eq(Ψ, H*exp(-u*u/2))
How can I do this with sympy?
Advertisement
Answer
Use the following function:
def variable_change(ODE,dependent_var, independent_var, new_dependent_var = None, new_independent_var= None, dependent_var_relation = None, independent_var_relation = None, order = 2): if new_dependent_var == None: new_dependent_var = dependent_var if new_independent_var == None: new_independent_var = independent_var # dependent variable change if new_independent_var != independent_var: for i in range(order, -1, -1): # remplace derivate a = D(dependent_var , independent_var, i ) ξ = Function("ξ")(independent_var) b = D( dependent_var.subs(independent_var, ξ), independent_var ,i) rel = solve(independent_var_relation, new_independent_var)[0] for j in range(order, 0, -1): b = b.subs( D(ξ,independent_var,j), D(rel,independent_var,j)) b = b.subs(ξ, new_independent_var) rel = solve(independent_var_relation, independent_var)[0] b = b.subs(independent_var, rel) ODE = ODE.subs(a,b) ODE = ODE.subs(independent_var, rel) # change of variables of indpendent variable if new_dependent_var != dependent_var: ODE = (ODE.subs(dependent_var.subs(independent_var,new_independent_var) , (solve(dependent_var_relation, dependent_var)[0]))) ODE = ODE.doit().expand() return ODE.simplify()
For the example posted:
from sympy import * from sympy import diff as D E, ℏ ,w,m,x,u = symbols("E, ℏ , w,m,x,u") Ψ ,H = map(Function, ["Ψ ","H"]) Ψ ,H = Ψ(x), H(u) DiffEq = Eq(-ℏ*ℏ*D(Ψ,x,2)/(2*m) + m*w*w*(x*x)*Ψ/2 - E*Ψ,0) display(DiffEq) display(Eq(u , x*sqrt(m*w/ℏ))) display(Eq(Ψ, H*exp(-u*u/2))) newODE = variable_change(ODE = DiffEq, independent_var = x, new_independent_var= u, independent_var_relation = Eq(u , x*sqrt(m*w/ℏ)), dependent_var = Ψ, new_dependent_var = H, dependent_var_relation = Eq(Ψ, H*exp(-u*u/2)), order = 2) display(newODE)
Under this substitution the differential equation outputted is then:
Eq((-E*H + u*w*ℏ*D(H, u) + w*ℏ*H/2 - w*ℏ*D(H, (u, 2))/2)*exp(-u**2/2), 0)