2

This code only works for solving the differential equation v_equation if v(t) isn't squared. When I squared it it returned the error PolynomialDivisionFailed. Is there another way of doing this with Sympy or should I find a different python package for doing these sorts of calculations.

from sympy import *
from matplotlib import pyplot as plt
import numpy as np

m = float(raw_input('Mass:\n> '))
g = 9.8
k = float(raw_input('Drag Coefficient:\n> '))
f1 = g * m
t = Symbol('t')
v = Function('v')
v_equation = dsolve(f1 - k * (v(t) ** 2) - m * Derivative(v(t)), 0)
C1 = Symbol('C1')
C1_ic = solve(v_equation.rhs.subs({t:0}),C1)[0]
v_equation = v_equation.subs({C1:C1_ic})

func = lambdify(t, v_equation.rhs,'numpy')
Kklj8
  • 137
  • 2
  • 11

1 Answers1

4

From my experience with symbolic math packages, I would not recommend performing (symbolic) calculations using floating point constants. It is better to define equations using symbolic constants, perform calculations as far as possible, and then substitute with numerical values.

With this approach, Sympy can provide a solution for this D.E.

First, define symbolic constants. To aid calculations, note that we can provide additional information about these constants (e.g., real, positive, e.t.c)

import sympy as sp

t = sp.symbols('t', real = True)
g, k, m = sp.symbols('g, k, m', real = True, positive = True)
v = sp.Function('v')

The symbolic solution for the DE can be obtained as follows

f1 = g * m
eq = f1 - k * (v(t) ** 2) - m * sp.Derivative(v(t))
sol = sp.dsolve(eq,v(t)).simplify()

The solution sol will be a function of k, m, g, and a constant C1. In general, there will be two, complex C1 values corresponding to the initial condition. However, both values of C1 result in the same (real-valued) solution when substituted in sol.

Note that if you don't need a symbolic solution, a numerical ODE solver, such as Scipy's odeint, may be used. The code would be the following (for an initial condition 0):

from scipy.integrate import odeint

def fun(v, t, m, k, g):
    return (g*m - k*v**2)/m

tn = np.linspace(0, 10, 101)
soln = odeint(fun, 0, tn, args=(1000, 0.2, 9.8))

soln is an array of samples v(t) corresponding to the tn elements

Stelios
  • 5,271
  • 1
  • 18
  • 32
  • Thanks, this helped a lot. Also, would I able to have the solution to the differential equation be returned as a function using scipy? – Kklj8 Aug 22 '16 at 14:28
  • @KKlj8 In principle, the answer is no, scipy (and any other *numerical* solver), can only provide samples of a function corresponding to a finite set of values of the independent variable. However, you can always use this set of samples as a basis for setting up a function that provides values for any other (arbitrary) value of the independent variable by means of [interpolation](http://docs.scipy.org/doc/scipy/reference/interpolate.html). – Stelios Aug 22 '16 at 14:40