CODEWITHSUNDEEP
pythonsympy
CrepeGoat
CrepeGoat

Reputation: 2515

SymPy dsolve with a parameter gives a wrong answer when the parameter is zero

My Code

I have a program that calculates the solutions to a 2nd-order differential equation, like in the code snippet below:

import sympy as sp
print('sympy version:', sp.__version__)

t = sp.symbols('t', real=True, nonnegative=True)
n = sp.symbols('n', integer=True, nonnegative=True)
f = sp.symbols('f', cls=sp.Function)

diff_eq = sp.Eq(f(t).diff(t, 2) + n**2*f(t), 0)

print('general solution:', sp.dsolve(diff_eq, f(t)))
print('solution at n=0 (pre-subs):', sp.dsolve(diff_eq.subs(n, 0), f(t)))
print('solution at n=0 (post-subs):', sp.dsolve(diff_eq, f(t)).subs(n, 0))

Results:

sympy version: 1.3

general solution: Eq(f(t), C1*sin(n*t) + C2*cos(n*t))

solution at n=0 (pre-subs): Eq(f(t), C1 + C2*t)

solution at n=0 (post-subs): Eq(f(t), C2)

My Problem

The solution form for a general n does not seem to accurately describe the specific solution form for n=0. Specifically, using dsolve first and subs(n, 0) second produces different results than using subs(n, 0) first and dsolve second, even though the two should be logically equivalent.

Can somebody explain the reason for the discrepancy in my results? Am I doing something wrong, or is this a bug?

Upvotes: 1

Views: 178

Answers (1)

user6655984
user6655984

Reputation:

It can be considered a bug in dsolve logic: it finds two eigenvalues n and -n and treats them as different without considering the special case n=0 when they are equal. Ideally it would output a Piecewise like the following code does. 

sol_nonzero = sp.dsolve(diff_eq, f(t)).rhs
sol_zero = sp.dsolve(diff_eq.subs(n, 0), f(t)).rhs
sol_complete = sp.Piecewise((sol_nonzero, sp.Ne(n, 0)), (sol_zero, True))
print('general solution:', sol_complete)
print('solution at n=0:', sol_complete.subs(n, 0))

This prints

general solution: Piecewise((C1*sin(n*t) + C2*cos(n*t), Ne(n, 0)), (C1 + C2*t, True))
solution at n=0: C1 + C2*t

A more familiar mathematical form is provided by sp.pprint(sol_complete). 

⎧C₁⋅sin(n⋅t) + C₂⋅cos(n⋅t)  for n ≠ 0
⎨
⎩        C₁ + C₂⋅t          otherwise

Upvotes: 1

Related Questions