Solving differential equations numerically¶
Analytic functions are certainly very nice ... but sometimes you just want to know the answer, eh?
References:
- wikipedia: Numerical Methods for ordinary differential equations
- wikipedia: Nondimensionalization
- Jupyter notebooks - a python environment for numerical
Jim Mahoney | Jan 2018
In [29]:
# I'm using python 2.
# Install some numerical and plotting packages. ("numpy" is "numerical python")
from numpy import *
from matplotlib.pyplot import plot
import matplotlib.pyplot as plt
import numpy as np
% matplotlib inline
In [30]:
# Define t as a list of 100 numbers from 0.0 to 10.0.
t = np.linspace(0, 50, num=512)
print "Here are the first few values of t :"
print t[0:10]
# Then to see for example sin(t), we can just take the sin of each of those.
# (There is an implied loop here - we're finding [sin(t[0]), sin(t[1]), ...]
print "And here are the first few values of sin(t) :"
print sin(t)[0:10]
print "And here's a picture ..."
plot( t, sin(t) , 'ro')
plt.show()
We can use these arrays of numbers to solve differentiatl equations.
Well, difference equations actually, but for small $\Delta t$ it's nearly the same. :)
As an example, let's solve
$$ \frac{dx}{dt} = \frac{1}{3+sin(\sqrt{t}))} = f(t) $$where x=1 at t=0. (I'm just making this up - it's an example, eh?)
We put t on a grid , something like
$$ t = [ 0.01, 0.02, 0.03, ... ] $$and approximate the derivative with a difference equation
$$ \frac{dx}{dt} = \frac{x[n+1] - x[n]}{\Delta t} $$This lets us solve for x[n+1] in terms of t and x[n] :
$$ x[n+1] = x[n] + (\Delta t) f(t) $$In [31]:
# Here we go.
x = zeros(len(t))
In [32]:
x
Out[32]:
In [33]:
dt = t[1] - t[0]
dt
Out[33]:
In [34]:
def f(t):
return 1.0/(3.0 + sin(sqrt(t)))
In [35]:
f(0.1234)
Out[35]:
In [36]:
# Start here :
x[0] = 1.0
# And now use the difference equation to get the next value ... over and over.
for i in range(1, len(t)):
_t = (i - 0.5) * dt # midway between x[i-1] and x[i]
x[i] = x[i-1] + dt * f(_t)
plot(t, x)
plt.show()
In [37]:
# All right, not a very exciting plot ... but that's the idea.
# The point is that you can do this sort of thing to get
# numerical solutions to equations that can't be done analytically
# ... which is most of 'em.
TODO
- more examples
- 2nd order equations
- some physics equations : forced damped harmonic oscillator, schrodinder equation, ...
- nondimensionalization (i.e. "get rid of the physical units")
In [ ]: