Ordinary Differential Equations

A good thing about Sympy that it has support for solving several kinds of ordinary differential equations via its dsolve command. We need to set up the ODE and pass it as the first argument, eq. The second argument is the function f(x) to solve for. An optimal third argument, hint, influences the method that dsolve uses: some methods are better-suited to certain classes of ODE’s, or will express the solution more simply, than others.

In order to set up ODE solver, we need a way to refer to the unknown function for which we a re solving as well as its derivatives. The Function and Derivative classes facilitate this:

>>> from sympy import Symbol, dsolve, Function, Derivative, Eq
>>> y = Function('y')
>>> x = Symbol('x')
>>> y_ = Derivative(y(x), x)
>>> print dsolve(y_ + 5*y(x), y(x))
y(x) == C1*exp(-5*x)

It is worth mentioning that dsolve has introduced a constant of integration C1. It will introduce as many constants as are required, and they will all be named Cn, whene n is an integer. Also the first argument to dsolve is taken to be equal to zero unless we use Eq() function to specify otherwise:

>>> print dsolve(y_ + 5*y(x), y(x))
y(x) == C1*exp(-5*x)
>>> print dsolve(Eq(y_ + 5*y(x), 0), y(x))
y(x) == C1*exp(-5*x)
>>> print dsolve(Eq(y_ + 5*y(x), 12), y(x))
y(x) == C1*exp(-5*x)/5 + 12/5

The results from dsolve are an instance of the Equality class. This has consequences when we wish to numerically evaluate the function and use the result elsewhere, because even after using subs() and evalf(), we still have an Equalitym not any sort of scalar. The way to evaluate the function to a number is via the rhs attribute of the Equality.

Note that here, we use z to store Equality returned by dsolve, even though it is an expression for a function called y(x), to emphasis the distinction between the Equality itself and the data it contains.

>>> z=dsolve(y_ + 5*y(x), y(x))
>>> print z
y(x) == C1*exp(-5*x)
>>> type(z)

>>> print z.rhs
C1*exp(-5*x)
>>> C1 = Symbol('C1')
>>> y3 = z.subs({C1:2, x:3})
>>> print y3
y(3) == 2*exp(-15)
>>> y3.evalf(10)
y(3) == 6.11804641e-7
>>> z.rhs.subs({C1:2, x:4}).evalf(10)
4.122307245e-9
>>> type(z.rhs.subs({C1:2, x:4}).evalf(10))

Sometimes the dsolve my return too general solution. One example is when there is possibility that some coefficients may be complex. If we know that, for example, they are always real and positive, we can provide dsolve this information to avoid the solution becoming unnecessarily complicated.

>> from sympy import *
>>> a, x = symbols('a, x')
>>> f = Function('f')
>>> print dsolve(Derivative(f(x), x, 2) - a**4*f(x), f(x))
f(x) == C1*exp(-x*sqrt(a**4)) + C2*exp(x*sqrt(a**4))
>>> print dsolve(Derivative(f(x), x, 4) - a**4*f(x), f(x))
f(x) == C1*exp(-x*(a**4)**(1/4)) + C2*exp(x*(a**4)**(1/4)) + C3*exp(-I*x*(a**4)**(1/4)) + C4*exp(I*x*(a**4)**(1/4))