# Ordinary differential equations


## Introduction 

Engineering systems are often represented as differential equations. ``sympy`` has robust support for the following types of ordinary differential equations (ODEs):

- 1st order separable differential equations
- 1st order differential equations whose coefficients or dx and dy are functions homogeneous of the same order.
- 1st order exact differential equations.
- 1st order linear differential equations
- 1st order Bernoulli differential equations.
- 2nd order Liouville differential equations.
- nth order linear homogeneous differential equation with constant coefficients.
- nth order linear inhomogeneous differential equation with constant coefficients using the method of undetermined coefficients.
- nth order linear inhomogeneous differential equation with constant coefficients using the method of variation of parameters.

In addition to ODE support, sympy can even solve separable (either multiplicative or additive) partial differential equations (PDEs).

The main functionality for ODE solving in sympy is the ``sympy.dsolve`` function. The call signature for this function is ``sympy.dsolve(eq, goal, **kwargs)``, where ``eq`` is the differential equation, ``goal`` is the function you want to end up with, and ``**kwargs`` is a placeholder for other arguments that could be passed to the function to help it out a bit. For more detail on what ``**kwargs`` could be, see the documentation.


```python
from sympy import *
```


```python
init_session()
plt.style.use(u"seaborn-notebook")
```

    IPython console for SymPy 1.0 (Python 2.7.13-64-bit) (ground types: python)
    
    These commands were executed:
    >>> from __future__ import division
    >>> from sympy import *
    >>> x, y, z, t = symbols('x y z t')
    >>> k, m, n = symbols('k m n', integer=True)
    >>> f, g, h = symbols('f g h', cls=Function)
    >>> init_printing()
    
    Documentation can be found at http://docs.sympy.org/1.0/


To solve differential equations, use dsolve. First, create an undefined function by passing cls=Function to the symbols function.


```python
f, g = symbols('f g', cls=Function)
```

Derivatives of f(x) are unevaluated.


```python
f(x).diff(x)
```

To represent the differential equation $f′′(x)−2f′(x)+f(x)=sin(x)$, we would use


```python
ode = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))
```


```python
ode
```

To solve the ODE, pass it and the function to solve for to ``dsolve``.


```python
sol = dsolve(ode, f(x))
display(sol)
```

``dsolve`` returns an instance of ``Eq``. This is because in general, solutions to differential equations cannot be solved explicitly for the function.


```python
dsolve(f(x).diff(x)*(1 - sin(f(x))), f(x))
```

The arbitrary constants in the solutions from ``dsolve`` are symbols of the form $C_1, C_2, C_3$, and so on.

## Applying Boundary conditions

This isn't implemented yet in ``dsolve`` stable version, but it will be available in the next release (probably!).

For now, solve for contants on your own. For example, if

$$f(0)=1\quad \left.\frac{d f}{d x}\right\vert_{x=0}=0,$$

solve the following equations:


```python
constants = solve([sol.rhs.subs(x,0) - 1, sol.rhs.diff(x,1).subs(x,0)- 0])
constants
```


```python
C1, C2 = symbols('C1,C2')
sol = sol.subs(constants)
sol
```

## Systems of equations


```python
g1, g2 = symbols('g1 g2', cls=Function)
```


```python
eq1 = g1(x).diff(x) - 2*g1(x) - g2(x)
eq2 = g2(x).diff(x) - 3*g1(x) - 4*g2(x)
```


```python
dsolve((eq1, eq2), fun=(g1(x), g2(x)))
```


```python

```


## Vector field


```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```

The '''Lotka–Volterra equations''', also known as the '''predator&ndash;prey equations''', are a pair of first-order, nonlinear, differential equations frequently used to describe the dynamics of biological systems in which two species interact, one as a predator and the other as prey. The populations change through time according to the pair of equations:

\begin{align}
\frac{dx}{dt} & = \alpha x - \beta x y \\[6pt]
\frac{dy}{dt} & = \delta x y  - \gamma y
\end{align}

where

* $x$ is the number of prey (for example, rabbits);
* $y$ is the number of some predator (for example, foxes);
* $\tfrac{dy}{dt}$ and $\tfrac{dx}{dt}$ represent the growth rates of the two populations over time;
* $t$  represents time; and
* $α, β, γ, δ$ are positive real parameters describing the interaction of the two species.


```python
y, x = np.mgrid[0:2.5:100j, 0:4:100j]
```


```python
alpha = 2/3
beta = 4/3
gamma = 1
delta = 1
dxdt = alpha*x - beta*x*y
dydt = delta*x*y - gamma*y
speed = np.sqrt(dxdt**2 + dydt**2)
```

In the first case we just plot the velocity field.


```python
fig0 = plt.figure(figsize=(8,5))
strm = plt.streamplot(x, y, dxdt, dydt, linewidth=2)
plt.xlabel("Preys")
plt.ylabel("Predators")
plt.xlim(0,4)
plt.ylim(0,2.5)
```

In this example we change the color of the line according to the magnitude of the
_speed_ at each point.


```python
fig1 = plt.figure(figsize=(8,5))
strm = plt.streamplot(x, y, dxdt, dydt, color=speed, linewidth=2, cmap="summer")
plt.colorbar(strm.lines)
plt.xlabel("Preys")
plt.ylabel("Predators")
plt.xlim(0,4)
plt.ylim(0,2.5)
```

In this case the densities in $x$ and $y$ directions are different.


```python
fig2 = plt.figure(figsize=(8,5))
plt.streamplot(x, y, dxdt, dydt, density=[0.5, 1])
plt.xlabel("Preys")
plt.ylabel("Predators")
plt.xlim(0,4)
plt.ylim(0,2.5)
```

In this last case we vary the width of the lines according to the magnitude
of the _speed_ at each point.


```python
fig3 = plt.figure(figsize=(8,5))
lw = 5*speed / speed.max()
plt.streamplot(x, y, dxdt, dydt, density=0.6, color='k', linewidth=lw)
plt.xlabel("Preys")
plt.ylabel("Predators")
plt.xlim(0,4)
plt.ylim(0,2.5)
```

## Power series approximation


```python
def picard_iteration(f, t, y0, n):
    y = y0
    for cont in range(n):
        y = y0 + f(t, y, n).integrate((t, 0, t))

    return y
```


```python
def fun(t, y, n):
    return Matrix([y[1],  t*y[0]])
```


```python
x0, x1 = symbols("x0 x1")
y = Matrix([airyai(0), airyaiprime(0)])
y_approx = picard_iteration(fun, t, y, 10)
```


```python
p0 = plot(airyai(t), (t, -4, 4), line_color="black");
p1 = plot(N(y_approx[0]), (t, -4, 4));
p0.extend(p1);
p0.show()
```

## References

1. SymPy Development Team (2016). [Sympy Tutorial: Matrices](http://docs.sympy.org/latest/tutorial/matrices.html)
2. Ivan Savov (2016). [Taming math and physics using SymPy](https://minireference.com/static/tutorials/sympy_tutorial.pdf)

The following cell change the style of the notebook.


```python
from IPython.core.display import HTML
def css_styling():
    styles = open('./styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>








```python

```

32bfd7887de96ff41ad6086e8175ac94c1db807c
ipynb
Jupyter Notebook
notebooks/sympy/ode.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
2017-06-29T17:45:20.000Z
2022-02-06T20:14:29.000Z
notebooks/sympy/ode.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
notebooks/sympy/ode.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
2019-04-22T08:08:56.000Z
2022-01-27T08:15:53.000Z
Qwen/Qwen-72B
1. YES
2. YES