7

So far I have managed to find the particular solution to this equation for any given mass and drag coefficient. I have not however found a way to plot the solution or even evaluate the solution for a specific point. I really want to find a way to plot the solution.

from sympy import *

m = float(raw_input('Mass:\n> '))
g = 9.8
k = float(raw_input('Drag Coefficient:\n> '))
f = Function('f')
f1 = g * m
t = Symbol('t')
v = Function('v')
equation = dsolve(f1 - k * v(t) - m * Derivative(v(t)), 0)
C1 = Symbol('C1')
C1_ic = solve(equation.rhs.subs({t:0}),C1)[0]
equation = equation.subs({C1:C1_ic})
benten
  • 1,995
  • 2
  • 23
  • 38
Kklj8
  • 137
  • 2
  • 11

3 Answers3

7

For completeness, you may also use Sympy's plot, which is probably more convenient if you want a "quick and dirty" plot.

plot(equation.rhs,(t,0,10))

enter image description here

Stelios
  • 5,271
  • 1
  • 18
  • 32
5

Import these libraries (seaborn just makes the plots pretty).

from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np

Then tack this onto the end. This will plot time, t, against velocity, v(t). 

# make a numpy-ready function from the sympy results
func = lambdify(t, equation.rhs,'numpy')
xvals = np.arange(0,10,.1)
yvals = func(xvals)

# make figure
fig, ax = plt.subplots(1,1,subplot_kw=dict(aspect='equal'))     
ax.plot(xvals, yvals)
ax.set_xlabel('t')
ax.set_ylabel('v(t)')
plt.show()

I get a plot like this for a mass of 2 and a drag coefficient of 2. enter image description here

benten
  • 1,995
  • 2
  • 23
  • 38
  • btw, [this question](http://stackoverflow.com/questions/10678843/evaluate-sympy-expression-from-an-array-of-values) is a helpful reference – benten Aug 16 '16 at 01:38
0

If I've understood correctly, you want to represent the right hand side of your solution, here's one of the multiple ways to do it:

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

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

t1 = np.arange(0.0, 50.0, 0.1)
y1 = [equation.subs({t: tt}).rhs for tt in t1]

plt.figure(1)
plt.plot(t1, y1)
plt.show()
BPL
  • 9,632
  • 9
  • 59
  • 117
  • 1
    It's better to use lambdify to create a numpy function than to evaluate the function using subs and a list comprehension. – asmeurer Aug 16 '16 at 16:50