# Solving multiple linear ordinary differential equations in Sympy

Presume we wish to solve the coupled linear ordinary differential equations given by

\begin{align}
3\dot{x} &= y\\
\dot{y}&=-3y-15x+4u(t)
\end{align}
where $u(t)$ is the step function and $x(0)=5$ and $y(0) = 10$. 




```python
%matplotlib inline

# import symbolic capability to Python- namespace is a better idea in a more general code. 
from sympy import *

# print things all pretty
from sympy.abc import *
init_printing()
```

Need to define variables as symbolic for sympy to use them. 


```python
t, C1, C2= symbols("t C1 C2")
x, y =  symbols("x y", cls = Function, Function = True)
```

First we must define the governing equations as equalities or expressions. We have the option, I will do one each way for the sake of illustration, as opposed to simplicity. 


```python
eq1 = Eq(3 * diff(x(t), t), y(t))
eq1
```

I'd like to use the Heaviside function. It doesn't work, but I'm showing how I tried anyway. 


```python
import sympy.functions.special.delta_functions as sfuncs
```


```python
eq2 = Eq(diff(y(t),t), - 3 * y(t) - 15 * x(t) + 4 * sfuncs.Heaviside(t))
# Note: Heaviside doesn't work for me in the solver, so I've left it here as "proper" but dokn't use it. 
eq2 = Eq(diff(y(t),t), - 3 * y(t) - 15 * x(t) + 4 * 1)
eq2
```

Solving the differential equations. The `ics = `... should apply the initial conditions. Doesn't work. Perhaps some day. 


```python
soln = dsolve((eq1, eq2), ics = {x: 5, y: 0})
soln
```

Solving for the constants. I'm substituting `t = 0`, then the initial values for `x` and `y`. You'll note that the first equation doesn't need `y(0)` substituted, while the second doesn't need `x(0)` substituted. I found this out after and did that for brevity.   


```python
constants = solve((soln[0].subs(t,0).subs(x(0),5), soln[1].subs(t,0).subs(y(0),10)),{C1,C2})
constants
```

Let's put in our constants and see what we get. I'm using `.rhs` to pull out the right side of the solution. You can look at `soln[0]` to see what I mean, or try `help(soln[0])` and read the results.  


```python
xsoln = expand(soln[0].rhs.subs(constants))
xsoln
```


```python
ysoln = soln[1].rhs.subs(constants)
ysoln
```


```python
eq1.subs(x(t),xsoln).subs(y(t),ysoln)
```

Is equation 1 true with this solution set?


```python
Eq(simplify(3*diff(xsoln,t)),simplify(ysoln))
```

Is equation 2 true with this solution set?


```python
Eq(simplify(diff(ysoln,t)),simplify(4-3*ysoln-15*xsoln))
```

So that's it. Out answers are indeed given by


```python
Eq(x(t),xsoln)
```

and


```python
Eq(y(t),ysoln)
```


```python

```


```python

```

d8e460c4d7bd23b541d8cbc41e00bfcba1179ff0
ipynb
Jupyter Notebook
Coupled_Differential_Equations.ipynb
josephcslater/JupyterExamples
4f2af75b9fda80dcab6cac0b9210713eeb703839
BSD-3-Clause
Coupled_Differential_Equations.ipynb
josephcslater/JupyterExamples
4f2af75b9fda80dcab6cac0b9210713eeb703839
BSD-3-Clause
Coupled_Differential_Equations.ipynb
josephcslater/JupyterExamples
4f2af75b9fda80dcab6cac0b9210713eeb703839
BSD-3-Clause
2020-02-07T20:28:39.000Z
2020-02-07T20:28:39.000Z
Qwen/Qwen-72B
1. YES
2. YES