I have a homogeneous second-order ODE, whose general solution is ur(r) = C1*u1(r) + C2*u2(r).
How can I get u1(r) and u2(r) as separate functions by setting (C1=1, C2=0) and (C1=0, C2=1), respectively?
This is what I tried (and other variants as well)
import IPython.display as disp
import sympy as sym
sym.init_printing(latex_mode='equation')
r = sym.symbols('r', real=True, positive=True)
ur = sym.Function('ur', real=True)
def Naxp(ur):
return r**2 * sym.diff(ur(r), r, 2) + r * sym.diff(ur(r), r) - ur(r)
Naxp_ur = Naxp(ur)
EqNaxpur = sym.Eq(Naxp_ur)
#print('Radial ODE:')
#display(Naxp_ur)
ur_sol = sym.dsolve(Naxp_ur, ur(r))
ur_sol1 = ur_sol.rhs
#display(sym.Eq(ur(r), ur_sol1))
# I still have to find how to extract the two basis functions
#u1 = sym.Function('u1', real=True)
C1, C2 = sym.symbols('C1, C2', real=True)
u1 = ur_sol1.subs([(C1, 1), (C2, 0)])
print(u1)
and I always get a result as if subs were ignored
C1/r + C2*r
Sympy isn't recognizing the equality between the constants it generated and the ones you generated. This seems to be due to setting your constants with real=True. Remove that and the code works.
C1, C2 = sym.symbols("C1 C2")
print(ur_sol1.subs({C1:1, C2:0})) # 1/r
print(ur_sol1.subs({C1:0, C2:1})) # r
There's probably a way to extract the constants and avoid defining them yourself, but I'm not sure the best way to do that. It could involve ur_sol1.free_symbols, but that also returns r. Edit: The sympy documentation seems to recommend creating the variables yourself, but there is probably still a better way.
Also, you might want to consider setting the initial conditions for your ODE to get a particular solution.