View source for SymPy/Differential Equations/RC Example

== Description ==
The following page is adapted from [[Maple/Differential Equations/RC Example]] and will go through an example of using SymPy's ability to symbolically solve a differential equation - specifically to analyze a circuit that undergoes a change in source values.  In this particular case, the independent source is given as a constant for all times before 0 sec, at which time it changes to a non-constant source.  The source is connected to a network containing both resistors and a capacitor. While there are several ways to put all this work together in Python, the following will provide a flexible framework for solving the equations and using substitution to determine a numerical list of initial conditions, substituting element values into the differential equations and putting them into a list, solving the differential equations, and finally plotting the results.  

'''''Note:''''' at the moment, SymPy's <code>dsolve</code> command works pretty well for up to two coupled ODEs but starts having issues once a third equation is in the mix.  Unlike [[Maple]], which has the option to use Laplace techniques with its <code>dsolve</code> command, SymPy does not.  Roundoff and other errors creep in, causing the solution to fail to match the initial conditions.  For higher-order systems, at the moment, a numerical solver such as <code>solve_ivp</code> might be required (which also requires breaking the system down into a series of first-order differential equations).

== Circuit ==
For this
demonstration code, the following circuit is used:
<center>
[[File:RCD DemoCircuit.png]]
</center>
where <math>R_1</math>=20 k<math>\Omega</math>,  
<math>R_2</math>=30 k<math>\Omega</math>, 
<math>C</math>=50 <math>\mu</math>F, 
and <math>v_s(t)</math> changes from 5 V to 10 cos(8<math>t</math>) V when <math>t=0</math> s. 

== DC Steady-State Analysis ==
Assuming the circuit has been in place for a "long time" before <math>t</math>=0 sec, and given the topology of the circuit and the fact that the independent source is a constant for all times before 0 sec, you can use the DC steady-state equivalent for the circuit at <math>t=0^-</math> sec:
<center>
[[File:RCD DemoCircuitSS.png]]
</center>
Using KCL at the top right node gives:
<center><math>
\frac{v_C(0^-)-v_S}{R_1}+
\frac{v_C(0^-)}{R_2}=0\!
</math></center>
which can be solved to find the capacitor voltage at the time just before the source <math>v_S(t)</math> changes.  

== Model Equations for t>0 ==
In general, after <math>t</math>=0 sec you can label the circuit as:
<center>
[[File:RCD DemoCircuitDE.png]]
</center>
Using KCL at the top right node again gives:
<center><math>
\frac{v_C(t)-v_S(t)}{R_1}+
\frac{v_C(t)}{R_2}+C\frac{dv_C(t)}{dt}=0\!
</math></center>
 

== Code ==
Now that the equations for DC steady state and for the differential model are known, you can write a Python script to solve for everything.

=== Initialization ===
You will need to import <code>sympy</code> to use symbolic calculations.  You will also need <code>numpy</code> and <code>matplotlib.pyplot</code> to make a graph.
<syntaxhighlight lang=python>
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
</syntaxhighlight>

Because you are solving systems of differential equations, the same helper functions that were used to solve [[SymPy/Simultaneous_Equations|simultaneous equations]] will be useful.  Given that, your script should import the helper functions from [https://tinyurl.com/mrgdemos/SymPy_Files/sym_helper.py <code>sym_helper.py</code>]; you can do this with:
<syntaxhighlight lang=python>
# Import the helper functions
try:
  from sym_helper import *
except: # Get sym_helper.py from GitHub
  import requests
  url = 'https://tinyurl.com/mrgdemos/SymPy_Files/sym_helper.py'
  open('sym_helper.py', 'w').write(requests.get(url).text)
  from sym_helper import *
</syntaxhighlight>

Finally, you will want the script to use $$\LaTeX$$ to display your equations if possible:
<syntaxhighlight lang=python>
# Use LaTeX if possible
sym.init_printing(use_latex=True)
</syntaxhighlight>

=== Declaring Variables and Functions ===
Next you should tell SymPy what your variables are.  You will also need to declare the ''functions' that you will be using later.  For this example, the known values include the two resistors, the capacitor, the DC value of the source before time 0, and time.  The functions you will want to declare include the (potentially non-constant) voltage of the source after time 0 and the voltage across the capacitor.  You may also want to create a list of the known variables and a list of the values of those variables:
<syntaxhighlight lang=python>
R1, R2, C, vsdc, t = sym.var('R1, R2, C, v_{s\,dc}, t')
vs = sym.Function("v_s")
vC = sym.Function("v_C")
knowns = R1, R2, C, vsdc, vs(t)
knownvals = [20000, 30000, 50e-6, 5, 10*sym.cos(8*t)]
</syntaxhighlight>


=== Initial Conditions From Steady-State ===
While SymPy can solve some differential equations with symbolic initial conditions and coefficients, most of the time this will result in a very unwieldy and unhelpful representation.  For this assignment, you will be providing SymPy with numerical values for the initial conditions.  To do this, we are going to ask SymPy to solve the linear algebra problem that comes from the DC steady state and then we will substitute in numerical values.  First, we will define a variable to be the initial condition of the capacitor voltage and will add it to a list of unknowns:
<syntaxhighlight lang=python>
vC0 = vC(t).subs(t, 0)
unksSS = [vC0]
</syntaxhighlight>
Next, we will write the conservation equation we have at DC steady-state and add it to a list of equations:
<syntaxhighlight lang=python>
eqn1SS = (vC0-vsdc) / R1 + vC0 / R2
eqnsSS = [eqn1SS]
</syntaxhighlight>
Finally, we can solve the equation and then substitute in numerical values:
<syntaxhighlight lang=python>
solnSS = sym.solve(eqnsSS, unksSS)
solnSSnum = makesubs(solnSS, knowns, knownvals)
</syntaxhighlight>
This results in <code>solnSSnum</code> being a dictionary with the initial value of $$v_C$$ defined in it, $$\{v_C(0):3\}$$

=== Differential Equations ===
Next we can set up and solve the differential equations.  First, make a list of all the unknown functions for which you want to solve:
<syntaxhighlight lang=python>
unks = [vC(t)]
</syntaxhighlight>

You can now write the differential equations (using <code>sym.Derivative(FUN, VAR)</code> to represent the derivatives) and put them in a list:
<syntaxhighlight lang=python>
eq1 = (vC(t)-vs(t))/R1 + vC(t)/R2 + C*sym.Derivative(vC(t), t)
eqns = [eq1]
</syntaxhighlight>

While it is possible for SymPy to solve some systems purely symbolically, there are times where it is much easier for it to solve equations where the coefficients are numbers.  To help with that, the <code>sym_helper</code> file has a command called <code>makesubseqnlist(LIST, KNOWNVARLIST, KNOWNVALLIST)</code> that will take a list of equations and substitute the given values in for the given variables.  The code for <code>makesubeqnlist</code> is:
<syntaxhighlight lang=python>
def makesubseqnlist(equations, vartosub, valtosub):
    sublist = list(zip(vartosub, valtosub))
    esub = []    
    for eqn in equations:
      esub.append(eqn.subs(sublist))
    return esub
</syntaxhighlight>
and the code that will go in your script for this is:
<syntaxhighlight lang=python>
neqns = makesubseqnlist(eqns, knowns, knownvals)
</syntaxhighlight>

Now you can solve the differential equations using <code>dsolve</code>:
<syntaxhighlight lang=python>
solnlist = sym.dsolve(neqns, unks, ics=solnSSnum)
</syntaxhighlight>

'''''Note:''''' unlike <code>solve</code>, which creates a dictionary of solutions, <code>dsolve</code> creates a list.  It turns out having a dictionary makes life much easier later, so <code>sym_helper</code> has a <code>make_dict(LIST)</code> command that will turn a list of equations into a dictionary.  The code for <code>make_dict</code> is:
<syntaxhighlight lang=python>
def makedict(sollist):
    d = {}
    for eqn in sollist:
        d[eqn.lhs] = eqn.rhs
    return d 
</syntaxhighlight>
and so the code you can add to your script is:
<syntaxhighlight lang=python>
def makedict(sollist):
soln = makedict(solnlist)
</syntaxhighlight>


=== Plotting ===
Plotting will require using <code>lambdify</code> to turn one or more functions into a lambda function and then setting up your plot.  Here is the code for this example that does all those things (making a plot that runs from $$t$$=0 to 10 sec):
== Description ==
The following page is adapted from [[Maple/Differential Equations/RC Example]] and will go through an example of using SymPy's ability to symbolically solve a differential equation - specifically to analyze a circuit that undergoes a change in source values.  In this particular case, the independent source is given as a constant for all times before 0 sec, at which time it changes to a non-constant source.  The source is connected to a network containing both resistors and a capacitor. While there are several ways to put all this work together in Python, the following will provide a flexible framework for solving the equations and using substitution to determine a numerical list of initial conditions, substituting element values into the differential equations and putting them into a list, solving the differential equations, and finally plotting the results.  

'''''Note:''''' at the moment, SymPy's <code>dsolve</code> command works pretty well for up to two coupled ODEs but starts having issues once a third equation is in the mix.  Unlike [[Maple]], which has the option to use Laplace techniques with its <code>dsolve</code> command, SymPy does not.  Roundoff and other errors creep in, causing the solution to fail to match the initial conditions.  For higher-order systems, at the moment, a numerical solver such as <code>solve_ivp</code> might be required (which also requires breaking the system down into a series of first-order differential equations).

== Circuit ==
For this
demonstration code, the following circuit is used:
<center>
[[File:RCD DemoCircuit.png]]
</center>
where <math>R_1</math>=20 k<math>\Omega</math>,  
<math>R_2</math>=30 k<math>\Omega</math>, 
<math>C</math>=50 <math>\mu</math>F, 
and <math>v_s(t)</math> changes from 5 V to 10 cos(8<math>t</math>) V when <math>t=0</math> s. 

== DC Steady-State Analysis ==
Assuming the circuit has been in place for a "long time" before <math>t</math>=0 sec, and given the topology of the circuit and the fact that the independent source is a constant for all times before 0 sec, you can use the DC steady-state equivalent for the circuit at <math>t=0^-</math> sec:
<center>
[[File:RCD DemoCircuitSS.png]]
</center>
Using KCL at the top right node gives:
<center><math>
\frac{v_C(0^-)-v_S}{R_1}+
\frac{v_C(0^-)}{R_2}=0\!
</math></center>
which can be solved to find the capacitor voltage at the time just before the source <math>v_S(t)</math> changes.  

== Model Equations for t>0 ==
In general, after <math>t</math>=0 sec you can label the circuit as:
<center>
[[File:RCD DemoCircuitDE.png]]
</center>
Using KCL at the top right node again gives:
<center><math>
\frac{v_C(t)-v_S(t)}{R_1}+
\frac{v_C(t)}{R_2}+C\frac{dv_C(t)}{dt}=0\!
</math></center>
 

== Code ==
Now that the equations for DC steady state and for the differential model are known, you can write a Python script to solve for everything.

=== Initialization ===
You will need to import <code>sympy</code> to use symbolic calculations.  You will also need <code>numpy</code> and <code>matplotlib.pyplot</code> to make a graph.
<syntaxhighlight lang=python>
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
</syntaxhighlight>

Because you are solving systems of differential equations, the same helper functions that were used to solve [[SymPy/Simultaneous_Equations|simultaneous equations]] will be useful.  Given that, your script should import the helper functions from [https://tinyurl.com/mrgdemos/SymPy_Files/sym_helper.py <code>sym_helper.py</code>]; you can do this with:
<syntaxhighlight lang=python>
# Import the helper functions
try:
  from sym_helper import *
except: # Get sym_helper.py from GitHub
  import requests
  url = 'https://tinyurl.com/mrgdemos/SymPy_Files/sym_helper.py'
  open('sym_helper.py', 'w').write(requests.get(url).text)
  from sym_helper import *
</syntaxhighlight>

Finally, you will want the script to use $$\LaTeX$$ to display your equations if possible:
<syntaxhighlight lang=python>
# Use LaTeX if possible
sym.init_printing(use_latex=True)
</syntaxhighlight>

=== Declaring Variables and Functions ===
Next you should tell SymPy what your variables are.  You will also need to declare the ''functions' that you will be using later.  For this example, the known values include the two resistors, the capacitor, the DC value of the source before time 0, and time.  The functions you will want to declare include the (potentially non-constant) voltage of the source after time 0 and the voltage across the capacitor.  You may also want to create a list of the known variables and a list of the values of those variables:
<syntaxhighlight lang=python>
R1, R2, C, vsdc, t = sym.var('R1, R2, C, v_{s\,dc}, t')
vs = sym.Function("v_s")
vC = sym.Function("v_C")
knowns = R1, R2, C, vsdc, vs(t)
knownvals = [20000, 30000, 50e-6, 5, 10*sym.cos(8*t)]
</syntaxhighlight>


=== Initial Conditions From Steady-State ===
While SymPy can solve some differential equations with symbolic initial conditions and coefficients, most of the time this will result in a very unwieldy and unhelpful representation.  For this assignment, you will be providing SymPy with numerical values for the initial conditions.  To do this, we are going to ask SymPy to solve the linear algebra problem that comes from the DC steady state and then we will substitute in numerical values.  First, we will define a variable to be the initial condition of the capacitor voltage and will add it to a list of unknowns:
<syntaxhighlight lang=python>
vC0 = vC(t).subs(t, 0)
unksSS = [vC0]
</syntaxhighlight>
Next, we will write the conservation equation we have at DC steady-state and add it to a list of equations:
<syntaxhighlight lang=python>
eqn1SS = (vC0-vsdc) / R1 + vC0 / R2
eqnsSS = [eqn1SS]
</syntaxhighlight>
Finally, we can solve the equation and then substitute in numerical values:
<syntaxhighlight lang=python>
solnSS = sym.solve(eqnsSS, unksSS)
solnSSnum = makesubs(solnSS, knowns, knownvals)
</syntaxhighlight>
This results in <code>solnSSnum</code> being a dictionary with the initial value of $$v_C$$ defined in it, $$\{v_C(0):3\}$$

=== Differential Equations ===
Next we can set up and solve the differential equations.  First, make a list of all the unknown functions for which you want to solve:
<syntaxhighlight lang=python>
unks = [vC(t)]
</syntaxhighlight>

You can now write the differential equations (using <code>sym.Derivative(FUN, VAR)</code> to represent the derivatives) and put them in a list:
<syntaxhighlight lang=python>
eq1 = (vC(t)-vs(t))/R1 + vC(t)/R2 + C*sym.Derivative(vC(t), t)
eqns = [eq1]
</syntaxhighlight>

While it is possible for SymPy to solve some systems purely symbolically, there are times where it is much easier for it to solve equations where the coefficients are numbers.  To help with that, the <code>sym_helper</code> file has a command called <code>makesubseqnlist(LIST, KNOWNVARLIST, KNOWNVALLIST)</code> that will take a list of equations and substitute the given values in for the given variables.  The code for <code>makesubeqnlist</code> is:
<syntaxhighlight lang=python>
def makesubseqnlist(equations, vartosub, valtosub):
    sublist = list(zip(vartosub, valtosub))
    esub = []    
    for eqn in equations:
      esub.append(eqn.subs(sublist))
    return esub
</syntaxhighlight>
and the code that will go in your script for this is:
<syntaxhighlight lang=python>
neqns = makesubseqnlist(eqns, knowns, knownvals)
</syntaxhighlight>

Now you can solve the differential equations using <code>dsolve</code>:
<syntaxhighlight lang=python>
solnlist = sym.dsolve(neqns, unks, ics=solnSSnum)
</syntaxhighlight>

'''''Note:''''' unlike <code>solve</code>, which creates a dictionary of solutions, <code>dsolve</code> creates a list.  It turns out having a dictionary makes life much easier later, so <code>sym_helper</code> has a <code>make_dict(LIST)</code> command that will turn a list of equations into a dictionary.  The code for <code>make_dict</code> is:
<syntaxhighlight lang=python>
vCfun = sym.lambdify(t, soln[vC(t)])
fig = plt.figure(num=1, clear=True)
ax = fig.add_subplot(1, 1, 1)
tn = np.linspace(0, 10, 1000)
ax.plot(tn, vCfun(tn))
</syntaxhighlight>


== Complete Example ==
For the Spring 2023 semester, students in EGR 224 get the code by going to the [https://colab.research.google.com/drive/17Y-DpLVSSn5um4v3w0MoH2UzRD1v90T4?usp=sharing Google Colab] version of it.  This notebook includes several more comments than are listed above and has multiple instances where it displays intermediate results.