Difference between revisions of "Maple/Differential Equations"

From PrattWiki
Jump to navigation Jump to search
(→‎Not Great Version 2)
 
(19 intermediate revisions by the same user not shown)
Line 10: Line 10:
 
=== Initialization ===
 
=== Initialization ===
 
As always, it is a good idea to include the <code>restart</code> command in your worksheet:
 
As always, it is a good idea to include the <code>restart</code> command in your worksheet:
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
 
restart
 
restart
 
</syntaxhighlight>
 
</syntaxhighlight>
  
 
=== Define Equation ===
 
=== Define Equation ===
−
In the same way that you were able to assign a linear algebra expression to a variable, you can do the same with a differential equation.  The key is to note that the Maple <code>diff</code> command can be used to calculate or represent a derivative.  You will need to explicitly let Maple know that your variable is a function (in our case, a function of $$t$$) by including that parameter with the variable.  Given that, you can store the differential equation in a variable called <code>eqn</code> with:
+
In the same way that you were able to assign a linear algebra expression to a variable, you can do the same with a differential equation.  The key is to note that the Maple <code>diff</code> command can be used to calculate or represent a derivative.  You will need to explicitly let Maple know that your variable is a function (in our case, a function of $$t$$) by including that parameter with the variable.  Given that, you can store the differential equation in a variable called <code>deqn1</code> with:
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
eqn1:=2*diff(y(t), t)+3*y(t)=4
+
deqn1:=2*diff(y(t), t)+3*y(t)=4
 
</syntaxhighlight>
 
</syntaxhighlight>
  
 
=== Solve Equation ===
 
=== Solve Equation ===
 
Solving a system of differential equations is also similar to solving a system of linear algebra equations - the main differences are that you will use <code>dsolve</code> instead of <code>solve</code>, you must continue to use $$y(t)$$ instead of just $$y$$, and you may end up needing to add some initial conditions.  The code  
 
Solving a system of differential equations is also similar to solving a system of linear algebra equations - the main differences are that you will use <code>dsolve</code> instead of <code>solve</code>, you must continue to use $$y(t)$$ instead of just $$y$$, and you may end up needing to add some initial conditions.  The code  
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
soln1:=dsolve([eqn1], [y(t)])
+
dsoln1:=dsolve([deqn1], [y(t)])
 
</syntaxhighlight>
 
</syntaxhighlight>
−
will result in a solution of:<center><math>\mathit{soln1}:=\left\{y\! \left(t\right)=\frac{4}{3}+{\mathrm e}^{-\frac{3 t}{2}} \textit{_}\mathit{C1}\right\}</math></center>
+
will result in a solution of:<center><math>\mathit{dsoln1}:=\left\{y\! \left(t\right)=\frac{4}{3}+{\mathrm e}^{-\frac{3 t}{2}} \textit{_}\mathit{C1}\right\}</math></center>
  
 
=== Initial Condition ===
 
=== Initial Condition ===
 
To incorporate initial conditions, you will give the <code>dsolve</code> command information about the value of the variable (or, for higher order differential equations, the value and values of the derivatives of the variable).  For example, to solve our sample equation with $$y(0)=5$$, you will include the initial condition by adding <code>y(0)=5</code> to the equations:
 
To incorporate initial conditions, you will give the <code>dsolve</code> command information about the value of the variable (or, for higher order differential equations, the value and values of the derivatives of the variable).  For example, to solve our sample equation with $$y(0)=5$$, you will include the initial condition by adding <code>y(0)=5</code> to the equations:
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
soln2:=dsolve([eqn1, y(0)=5], [y(t)])
+
dsoln2:=dsolve([deqn1, y(0)=5], [y(t)])
 
</syntaxhighlight>
 
</syntaxhighlight>
 
will produce <center><math>
 
will produce <center><math>
−
\mathit{soln2}:=y\! \left(t\right)=\frac{4}{3}+\frac{11 {\mathrm e}^{-\frac{3 t}{2}}}{3}
+
\mathit{dsoln2}:=y\! \left(t\right)=\frac{4}{3}+\frac{11 {\mathrm e}^{-\frac{3 t}{2}}}{3}
 
</math></center>
 
</math></center>
  
 
Note that the initial condition does not have to be at time 0; if you know that $$y(6)=7$$, you can use that as well:
 
Note that the initial condition does not have to be at time 0; if you know that $$y(6)=7$$, you can use that as well:
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
soln3:=dsolve([eqn1, y(6)=7], [y(t)])
+
dsoln3:=dsolve([deqn1, y(6)=7], [y(t)])
 
</syntaxhighlight>
 
</syntaxhighlight>
 
will produce<center><math>
 
will produce<center><math>
−
\mathit{soln3}:=y\! \left(t\right)=\frac{4}{3}+\frac{17 {\mathrm e}^{-\frac{3 t}{2}}}{3}
+
\mathit{dsoln3}:=y\! \left(t\right)=\frac{4}{3}+\frac{17 {\mathrm e}^{-\frac{3 t}{2}}e^9}{3}
 
</math></center>
 
</math></center>
  
 
=== Making a Plot ===
 
=== Making a Plot ===
−
Once you have a solution or set of solutions, you can plot them using subs.  For instance, to plot $$y(t)$$ in ''soln2'' above (where we set $$y(0)=5$$) you can use:
+
Once you have a solution or set of solutions, you can plot them using subs.  For instance, to plot $$y(t)$$ in ''dsoln2'' above (where we set $$y(0)=5$$) you can use:
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
plot(subs(soln2, y(t)), t = 0 .. 5)
+
plot(subs(dsoln2, y(t)), t = 0 .. 5)
 
</syntaxhighlight>
 
</syntaxhighlight>
−
and can of course add other plotting options as needed.  Note that ''soln2'' is a single equation, not a collection.  If you end up getting a collection of results, you may need to de-bracket them.
+
and can of course add other plotting options as needed.  Note that ''dsoln2'' is a single equation, not a collection.  If you end up getting a collection of results, you may need to de-bracket them.
  
 
== Second Order Example ==
 
== Second Order Example ==
 
Now imagine that you want to solve for $$y(t)$$ in<center><math>
 
Now imagine that you want to solve for $$y(t)$$ in<center><math>
−
\frac{d^2y(t)}{dt^2}=-g</math></center>. We can solve this using all symbols or numbers, and we can solve with or without initial conditions.
+
\frac{d^2y(t)}{dt^2}=-g</math></center>   
 +
We can solve this using all symbols or numbers, and we can solve with or without initial conditions.
  
 
=== Initialization ===
 
=== Initialization ===
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
 
restart
 
restart
 
</syntaxhighlight>
 
</syntaxhighlight>
  
 
=== Define Equation ===
 
=== Define Equation ===
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
eqn1 := diff(y(t), t, t) = -g
+
deqn1 := diff(y(t), t, t) = -g
 
</syntaxhighlight>
 
</syntaxhighlight>
  
 
=== Solve Equation / Initial Conditions ===
 
=== Solve Equation / Initial Conditions ===
 
Without initial conditions, you can use
 
Without initial conditions, you can use
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
soln1 := dsolve([eqn1], [y(t)])
+
dsoln1 := dsolve([deqn1], [y(t)])
 
</syntaxhighlight>
 
</syntaxhighlight>
−
to get:<center><math>\mathit{soln1}:=\left\{y\! \left(t\right)=-\frac{1}{2} g \,t^{2}+\textit{_}\mathit{C1} t+\textit{_}\mathit{C2}\right\}</math></center>
+
to get:<center><math>\mathit{dsoln1}:=\left\{y\! \left(t\right)=-\frac{1}{2} g \,t^{2}+\textit{_}\mathit{C1} t+\textit{_}\mathit{C2}\right\}</math></center>
  
 
You can also put in symbolic initial conditions.  To put in a derivative condition, use the format <code>D^n(var)(time)=val</code> to establish a condition for the ''n''th derivative of <code>var</code> at time <code>time</code>.  For instance, if you know some initial velocity and position, we can use:
 
You can also put in symbolic initial conditions.  To put in a derivative condition, use the format <code>D^n(var)(time)=val</code> to establish a condition for the ''n''th derivative of <code>var</code> at time <code>time</code>.  For instance, if you know some initial velocity and position, we can use:
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
soln2 := dsolve([eqn1, y(0) = y0, D(y)(0) = v0], [y(t)])
+
dsoln2 := dsolve([deqn1, y(0) = y0, D(y)(0) = v0], [y(t)])
−
</syntaxhighlight>to get:<center><math>\mathit{soln2}:=y\! \left(t\right)=-\frac{1}{2} g \,t^{2}+\mathit{v0} t+\mathit{y0}</math></center>
+
</syntaxhighlight>to get:<center><math>\mathit{dsoln2}:=y\! \left(t\right)=-\frac{1}{2} g \,t^{2}+\mathit{v0} t+\mathit{y0}</math></center>
  
  
 
=== Making a Plot ===
 
=== Making a Plot ===
 
Since your answer has symbols, you will need to replace them with numbers before plotting:
 
Since your answer has symbols, you will need to replace them with numbers before plotting:
−
<syntaxhighlight lang=Maple>
+
<syntaxhighlight lang=text>
−
plot(subs(soln2, g = 9.81, y0 = 5, v0 = 10, y(t)), t = 0 .. 2)
+
plot(subs(dsoln2, g = 9.81, y0 = 5, v0 = 10, y(t)), t = 0 .. 2)
 
</syntaxhighlight>
 
</syntaxhighlight>
−
and can of course add other plotting options as needed.  Note that ''soln2'' is a single equation, not a collection.  If you end up getting a collection of results, you may need to de-bracket them.
+
and can of course add other plotting options as needed.  Once again, note that ''dsoln2'' is a single equation, not a collection, so there is no need to de-bracket.
  
  
−
<!--
+
== More Complicated Coupled Example ==
−
== Basic Symbolic Example ==
+
For a series RLC circuit with an applied voltage $$v_S(t)$$ across the series network, you can find a coupled set of differential equations using KVL and the model equation for the capacitor, respectively, as:<center><math>
−
For some simple systems of differential equations, SymPy can actually solve with symbolic coefficients and undefined forcing functions.  The result will be an unevaluated integral, but still - pretty impressive!  Here is an example that solves $$a\,\frac{dy(t)}{dt} + b\,y(t)=x(t)$$ assuming $$y(t_0)=y_0$$.  If you want to follow along, the code is in the <code>dsolve_2</code> notebook in the [https://drive.google.com/drive/folders/10Y1kzpu-DTTRwm1mYtdgTwpXF2zKpmoL?usp=sharing Basic Dsolve] Google Drive folder.
+
\begin{align*}
−
<syntaxhighlight lang=python>
+
L\frac{i_L(t)}{st}+R\,i_L(t)+v_C(t)&=0\\
−
import sympy as sym
+
i_L(t)&=C\frac{dv_C(t)}{dt}
−
sym.init_printing(use_latex=True)
+
\end{align*}</math></center>
−
a, b, t, y0, t0 = sym.var('a, b, t, y0, t0')
+
We can '''try''' to solve these symbolically, but the results will not be particularly helpful.
−
x = sym.Function('x')
+
 
−
y = sym.Function('y')
+
=== Initialization ===
−
eqn1 = sym.Eq(a*sym.Derivative(y(t), t) + b * y(t), x(t))
+
<syntaxhighlight lang=text>
−
eqns = [eqn1]
+
restart
−
soln = sym.dsolve(eqns, [y(t)], ics = {y(t0): y0})
 
−
soln
 
 
</syntaxhighlight>
 
</syntaxhighlight>
−
The solution is:
 
−
<center>
 
−
$$
 
−
\left[ y{\left(t \right)} = \left(y_{0} e^{\frac{b t_{0}}{a}} - \frac{\int x{\left(t_{0} \right)} e^{\frac{b t_{0}}{a}}\, dt_{0}}{a}\right) e^{- \frac{b t}{a}} + \frac{e^{- \frac{b t}{a}} \int x{\left(t \right)} e^{\frac{b t}{a}}\, dt}{a}\right]
 
−
$$
 
−
</center>
 
−
=== Making Substitutions ===
 
−
Making substitutions into lists of equations is a little complicated - you need to go into each element in the list and make the substitutions there.  Given that, there is a helper function in <code>sym_helper</code> called <code>makesubseqnlist(EQN LIST, VARS TO SUB, VALS TO SUB)</code> that will take a list of equations, a list of variables (or functions), and a list of values (or functions).  The code for <code>makesubseqnlist(EQN LIST, VARS TO SUB, VALS TO SUB)</code> is
 
  
−
<syntaxhighlight lang=python>
+
=== Define Equations ===
−
def makesubseqnlist(equations, vartosub, valtosub):
+
<syntaxhighlight lang=text>
−
    sublist = list(zip(vartosub, valtosub))
+
deqn1 := -vs(t) + L*diff(iL(t), t) + R*iL(t) + vC(t) = 0;
−
    esub = []   
+
deqn2 := iL(t) = C*diff(vC(t), t);
−
    for eqn in equations:
+
deqns := [deqn1, deqn2]
−
      esub.append(eqn.subs(sublist))
 
−
    return esub
 
 
</syntaxhighlight>
 
</syntaxhighlight>
  
−
and you can get by by adding the following to near the top of your notebook or script:
+
=== Solve Equation / Initial Conditions ===
−
<syntaxhighlight lang=python>
+
==== Not Great Version 1 ====
−
try:
+
If you try
−
  from sym_helper import *
+
<syntaxhighlight lang=text>
−
except: # Get sym_helper.py from GitHub
+
dsoln1 := dsolve(deqns, [iL(t), vC(t)])
−
  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>
 
</syntaxhighlight>
 +
<div class="mw-collapsible mw-collapsed">
 +
you will get something huge - click "Expand" at right to see it->
 +
<div class="mw-collapsible-content">
 +
<center>
 +
<math>
 +
\mathit{dsoln1}:= \left\{\mathit{iL}\! \left(t\right) = \left(\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right)+\left(-\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right)+\left(-\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}+\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \textit{_}\mathit{C1} {\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}+\left(\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}-\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \textit{_}\mathit{C2} {\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}},
 +
\\
 +
\mathit{vC}\! \left(t\right)=-\frac{-{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \textit{_}\mathit{C2} \sqrt{C \left(R^{2} C-4 L\right)}-{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \textit{_}\mathit{C1} \sqrt{C \left(R^{2} C-4 L\right)}+\left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}-\left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{\sqrt{C \left(R^{2} C-4 L\right)}}\mathrm{
 +
\\}\right\}
 +
</math>
 +
</center></div>
 +
</div>
  
−
If we specifically want to set $$a=2$$, $$b=3$$, and $$x(t)=4\,\cos(5t)$$, and then solve that where $$y(0)=6$$, we can:
+
==== Not Great Version 2 ====
−
<syntaxhighlight lang=python>
+
Even if you put in initial conditions with
−
eqns2 = makesubseqnlist(eqns, [a, b, x(t)], [2, 3, 4*sym.cos(5*t)])
+
<syntaxhighlight lang=text>
−
soln2 = sym.dsolve(eqns2, [y(t)], ics = {y(0): 6})
+
dsoln2 := dsolve([deqns[], iL(0) = iL0, vC(0) = vC0], [iL(t), vC(t)])
 
</syntaxhighlight>
 
</syntaxhighlight>
−
yields:
+
<div class="mw-collapsible mw-collapsed">
 +
you will still get something huge - click "Expand" at right to see it ->
 +
<div class="mw-collapsible-content">
 
<center>
 
<center>
−
$$
+
<math>
−
\left[ y{\left(t \right)} = \frac{40 \sin{\left(5 t \right)}}{109} + \frac{12 \cos{\left(5 t \right)}}{109} + \frac{642 e^{- \frac{3 t}{2}}}{109}\right]
+
\mathit{dsoln2}:= \left\{\mathit{iL}\, \left(t\right) = \left(\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right)+\left(-\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right)-\frac{\left(-\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}+\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \left(R C \mathit{vC0}+2 L \mathit{iL0}-\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, {\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 C \left(R^{2} C-4 L\right)}+\frac{\left(\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}-\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \left(R C \mathit{vC0}+2 L \mathit{iL0}+\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, {\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 C \left(R^{2} C-4 L\right)},
−
$$
+
\\
−
</center>
+
\mathit{vC}\, \left(t\right)=\frac{-\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \left(R C \mathit{vC0}+2 L \mathit{iL0}-\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, \sqrt{C \left(R^{2} C-4 L\right)}}{2 C \left(R^{2} C-4 L\right)}+\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \left(R C \mathit{vC0}+2 L \mathit{iL0}+\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, \sqrt{C \left(R^{2} C-4 L\right)}}{2 C \left(R^{2} C-4 L\right)}-\left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}+\left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{\sqrt{C \left(R^{2} C-4 L\right)}}\mathrm{
 +
\\}\right\}
 +
</math>
 +
</center></div>
 +
</div>
  
−
== Making a Dictionary ==
+
==== Numerical Version ====
−
For linear algebra, the results came in a dictionary that had each of the unknowns as keys and their solutions as variables.  Since <code>dsolve</code> produces a list, versus a ductionary, it will be useful to convert that list to a dictionary for making substitutions later.  <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:
+
Instead, what you will want to do is substitute in numerical values first, then get the solution.  For second-order and higher systems, the types of responses can include exponentials, sinusoids, exponential sinusoids, and polynomials.  Without knowing the relative values of the symbolic components, Maple cannot easily simplify.  For this example, assume that we have a 1000 $\Omega$ resistor, a 100 nF capacitor, and a 100 mH inductor with a source that turns on to 5 V at time 0.  Further assume that we know that the initial current in the inductor is 1 mA and the initial voltage across the capacitor is -2 V.  We can get numerical versions of the equations and use them with the numerical versions of the initial conditions as follows:
−
<syntaxhighlight lang=python>
+
   
−
def makedict(sollist):
+
<syntaxhighlight lang=text>
−
    d = {}
+
vals := R = 1000, C = 0.100*10^(-6), L = 0.100, vs(t) = 5;
−
    for eqn in sollist:
+
numdeqns := subs(vals, deqns)[];
−
        d[eqn.lhs] = eqn.rhs
+
dsoln3 := dsolve([numdeqns, iL(0) = 0.001, vC(0) = -2], [iL(t), vC(t)])
−
    return d
 
 
</syntaxhighlight>
 
</syntaxhighlight>
−
and so the code you can add to your script is:
+
<div class="mw-collapsible mw-collapsed">
−
<syntaxhighlight lang=python>
+
You will get something less...complicated but still not small - click "Expand" at right to see it ->
−
def makedict(sollist):
+
<div class="mw-collapsible-content">
−
soln = makedict(solnlist)
+
<center>
−
</syntaxhighlight>
+
<math>
−
and, as long as <code>sym_helper</code> has been imported, you can run it by writing code like:
+
\mathit{dsoln3}:= \left\{\mathit{iL}\! \left(t\right) = \frac{{\mathrm e}^{-5000 t} \left(2 \cos\! \left(5000 t \sqrt{3}\right)+\frac{26 \sqrt{3}\, \sin\left(5000 t \sqrt{3}\right)}{3}\right)}{2000},
−
<syntaxhighlight lang=python>
+
\\
−
solndict = makedict(soln2)
+
\mathit{vC}\! \left(t\right)=5+{\mathrm e}^{-5000 t} \left(-7 \cos\! \left(5000 t \sqrt{3}\right)-\frac{5 \sqrt{3}\, \sin\! \left(5000 t \sqrt{3}\right)}{3}\right)\right\}
 +
</math>
 +
</center></div>
 +
</div>
 +
You can look at a version where Maple does floating point evaluation and rounds things; to see that version, you can type:
 +
<syntaxhighlight lang=text>
 +
evalf[4](dsoln3)
 
</syntaxhighlight>
 
</syntaxhighlight>
−
which yields:
+
and will see<center>
−
<center>
+
<math>
−
$$
+
\mathit{iL}\! \left(t\right) =  0.0005000 {\mathrm e}^{- 5000.0 t} \left( 2.0 \cos\! \left( 8660.0 t\right)+ 15.01 \sin\! \left( 8660.0 t\right)\right),
−
\left\{ y{\left(t \right)} : \frac{40 \sin{\left(5 t \right)}}{109} + \frac{12 \cos{\left(5 t \right)}}{109} + \frac{642 e^{- \frac{3 t}{2}}}{109}\right\}
+
\\
−
 
+
\mathit{vC}\! \left(t\right)= 5.0+{\mathrm e}^{- 5000.0 t} \left(- 7.0 \cos\! \left( 8660.0 t\right)- 2.887 \sin\! \left( 8660.0 t\right)\right)
−
$$
+
</math>
 
</center>
 
</center>
  
−
== Plotting Results ==
+
==== Making Plots ====
−
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); note the additional required imports:
+
Once you have the numerical solutions, you can make plots:
−
<syntaxhighlight lang=python>
+
<syntaxhighlight lang=text>
−
import matplotlib.pyplot as plt
+
plot(subs(dsoln3[], iL(t)), t = 0 .. 0.001);
−
import numpy as np
+
plot(subs(dsoln3[], vC(t)), t = 0 .. 0.001)
−
yfun = sym.lambdify(t, solndict[y(t)])
 
−
fig = plt.figure(num=1, clear=True)
 
−
ax = fig.add_subplot(1, 1, 1)
 
−
tn = np.linspace(0, 10, 1000)
 
−
ax.plot(tn, yfun(tn))
 
 
</syntaxhighlight>
 
</syntaxhighlight>
−
 
−
== Lessons Learned ==
 
−
* dsolve is good for up to two coupled ODEs.  Once there are three or more...not so much (even if all the coefficients are ints or floats).  Once there are three or more (linear) coupled differential equations, consider going the Laplace route.
 
−
* If your equations are a combination of differential equations a non-differential equations, dsolve expects all the differential equations to come first.
 
−
* If you are using certain functions from SymPy that aren't exactly the same in NumPy, when you lambdify things you will need to provide the translation
 
−
== Examples ==
 
−
* See [[SymPy/Differential Equations/RC Example]] for a complete example; this will reference IVP_RCDemo in the Google Drive folder.
 
−
* See [[SymPy/CAD]] for a list of SymPy examples from a free textbook.  The ones most relevant to this page are from Chapter 5.
 
−
* There are several other files in the [https://drive.google.com/drive/folders/10Y1kzpu-DTTRwm1mYtdgTwpXF2zKpmoL?usp=sharing Basic Dsolve] Google Drive folder:
 
−
 
−
** Coupled ODEs (just to see how SymPy handles those!)
 
−
*** Coupled2
 
−
*** Coupled3
 
−
*** deq_coupled: symbolically solve two coupled first-order ODEs; substitute numbers later to make graph.  Needs some editing - does not completely follow the general template.
 
−
*** deq_coupled_num: symbolically solve two coupled first-order ODEs with numbers substituted in before solving. Needs some editing - does not completely follow the general template.
 
−
** First-order differential equations:
 
−
*** IVP_demo1st: symbolic solution to 1st order ODE with a constant forcing function; substitutes numbers later to make graph
 
−
*** IVP_demo1stfun: symbolic solution to 1st order ODE with a non-constant forcing function; substitutes numbers later to make graph (including forcing function on graph)
 
−
*** IVP_demo1st_pulse: symbolic solution to 1st order ODE with a non-constant forcing function comprised of step and ramp functions; substitutes numbers later to make graph (including forcing function on graph)
 
−
**** This one shows extra code needed to translate step and ramp functions into a lambdify function
 
−
** Second-order differential equations:
 
−
*** IVP_demo2ndcon: symbolic solution to 2nd order ODE with a constant forcing function; substitutes numbers later to make graph.  Needs some editing - does not completely follow the general template.
 
−
*** IVP_demo2ndfun: symbolic solution to 2nd order ODE with a non-constant forcing function; substitutes numbers later to make graph.  Needs some editing - does not completely follow the general template.
 
−
-->
 

Latest revision as of 21:00, 28 February 2024

Introduction

This page focuses on using Maple to find both the symbolic and the numeric solutions to differential equations obtained from electric circuits.

Note: There is a...decidedly more complicated explanation of these things at Maple/Differential Equations/Old; the goal for Spring 2024 and beyond is to keep things simpler.

Very Basic Example

Imagine you have the equation $$2\frac{dy(t)}{dt} + 3 y(t) = 4$$ with the initial condition $$y(0)=5$$, and you want to solve for $$y(t)$$. You can do this with Maple as follows:

Initialization

As always, it is a good idea to include the restart command in your worksheet: 

restart

Define Equation

In the same way that you were able to assign a linear algebra expression to a variable, you can do the same with a differential equation. The key is to note that the Maple diff command can be used to calculate or represent a derivative. You will need to explicitly let Maple know that your variable is a function (in our case, a function of $$t$$) by including that parameter with the variable. Given that, you can store the differential equation in a variable called deqn1 with: 

deqn1:=2*diff(y(t), t)+3*y(t)=4

Solve Equation

Solving a system of differential equations is also similar to solving a system of linear algebra equations - the main differences are that you will use dsolve instead of solve, you must continue to use $$y(t)$$ instead of just $$y$$, and you may end up needing to add some initial conditions. The code 

dsoln1:=dsolve([deqn1], [y(t)])

will result in a solution of:

\(\mathit{dsoln1}:=\left\{y\! \left(t\right)=\frac{4}{3}+{\mathrm e}^{-\frac{3 t}{2}} \textit{_}\mathit{C1}\right\}\)

Initial Condition

To incorporate initial conditions, you will give the dsolve command information about the value of the variable (or, for higher order differential equations, the value and values of the derivatives of the variable). For example, to solve our sample equation with $$y(0)=5$$, you will include the initial condition by adding y(0)=5 to the equations: 

dsoln2:=dsolve([deqn1, y(0)=5], [y(t)])

will produce

\( \mathit{dsoln2}:=y\! \left(t\right)=\frac{4}{3}+\frac{11 {\mathrm e}^{-\frac{3 t}{2}}}{3} \)

Note that the initial condition does not have to be at time 0; if you know that $$y(6)=7$$, you can use that as well: 

dsoln3:=dsolve([deqn1, y(6)=7], [y(t)])

will produce

\( \mathit{dsoln3}:=y\! \left(t\right)=\frac{4}{3}+\frac{17 {\mathrm e}^{-\frac{3 t}{2}}e^9}{3} \)

Making a Plot

Once you have a solution or set of solutions, you can plot them using subs. For instance, to plot $$y(t)$$ in dsoln2 above (where we set $$y(0)=5$$) you can use: 

plot(subs(dsoln2, y(t)), t = 0 .. 5)

and can of course add other plotting options as needed. Note that dsoln2 is a single equation, not a collection. If you end up getting a collection of results, you may need to de-bracket them.

Second Order Example

Now imagine that you want to solve for $$y(t)$$ in

\( \frac{d^2y(t)}{dt^2}=-g\)

We can solve this using all symbols or numbers, and we can solve with or without initial conditions.

Initialization

restart

Define Equation

deqn1 := diff(y(t), t, t) = -g

Solve Equation / Initial Conditions

Without initial conditions, you can use 

dsoln1 := dsolve([deqn1], [y(t)])

to get:

\(\mathit{dsoln1}:=\left\{y\! \left(t\right)=-\frac{1}{2} g \,t^{2}+\textit{_}\mathit{C1} t+\textit{_}\mathit{C2}\right\}\)

You can also put in symbolic initial conditions. To put in a derivative condition, use the format D^n(var)(time)=val to establish a condition for the nth derivative of var at time time. For instance, if you know some initial velocity and position, we can use: 

dsoln2 := dsolve([deqn1, y(0) = y0, D(y)(0) = v0], [y(t)])

to get:

\(\mathit{dsoln2}:=y\! \left(t\right)=-\frac{1}{2} g \,t^{2}+\mathit{v0} t+\mathit{y0}\)


Making a Plot

Since your answer has symbols, you will need to replace them with numbers before plotting: 

plot(subs(dsoln2, g = 9.81, y0 = 5, v0 = 10, y(t)), t = 0 .. 2)

and can of course add other plotting options as needed. Once again, note that dsoln2 is a single equation, not a collection, so there is no need to de-bracket.


More Complicated Coupled Example

For a series RLC circuit with an applied voltage $$v_S(t)$$ across the series network, you can find a coupled set of differential equations using KVL and the model equation for the capacitor, respectively, as:

\( \begin{align*} L\frac{i_L(t)}{st}+R\,i_L(t)+v_C(t)&=0\\ i_L(t)&=C\frac{dv_C(t)}{dt} \end{align*}\)

We can try to solve these symbolically, but the results will not be particularly helpful.

Initialization

restart

Define Equations

deqn1 := -vs(t) + L*diff(iL(t), t) + R*iL(t) + vC(t) = 0;
deqn2 := iL(t) = C*diff(vC(t), t);
deqns := [deqn1, deqn2]

Solve Equation / Initial Conditions

Not Great Version 1

If you try 

dsoln1 := dsolve(deqns, [iL(t), vC(t)])

you will get something huge - click "Expand" at right to see it->

\( \mathit{dsoln1}:= \left\{\mathit{iL}\! \left(t\right) = \left(\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right)+\left(-\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right)+\left(-\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}+\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \textit{_}\mathit{C1} {\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}+\left(\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}-\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \textit{_}\mathit{C2} {\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}, \\ \mathit{vC}\! \left(t\right)=-\frac{-{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \textit{_}\mathit{C2} \sqrt{C \left(R^{2} C-4 L\right)}-{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \textit{_}\mathit{C1} \sqrt{C \left(R^{2} C-4 L\right)}+\left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}-\left({\int}\mathit{vS}\! \left(t\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}{d}t\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{\sqrt{C \left(R^{2} C-4 L\right)}}\mathrm{ \\}\right\} \)

Not Great Version 2

Even if you put in initial conditions with 

dsoln2 := dsolve([deqns[], iL(0) = iL0, vC(0) = vC0], [iL(t), vC(t)])

you will still get something huge - click "Expand" at right to see it ->

\( \mathit{dsoln2}:= \left\{\mathit{iL}\, \left(t\right) = \left(\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right)+\left(-\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} C R}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}+\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 L}\right) \left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right)-\frac{\left(-\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}+\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \left(R C \mathit{vC0}+2 L \mathit{iL0}-\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, {\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 C \left(R^{2} C-4 L\right)}+\frac{\left(\frac{C^{2} R^{2}}{2 \sqrt{C \left(R^{2} C-4 L\right)}\, L}-\frac{R C}{2 L}-\frac{2 C}{\sqrt{C \left(R^{2} C-4 L\right)}}\right) \left(R C \mathit{vC0}+2 L \mathit{iL0}+\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, {\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{2 C \left(R^{2} C-4 L\right)}, \\ \mathit{vC}\, \left(t\right)=\frac{-\frac{{\mathrm e}^{-\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \left(R C \mathit{vC0}+2 L \mathit{iL0}-\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, \sqrt{C \left(R^{2} C-4 L\right)}}{2 C \left(R^{2} C-4 L\right)}+\frac{{\mathrm e}^{-\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}} \left(R C \mathit{vC0}+2 L \mathit{iL0}+\sqrt{R^{2} C^{2}-4 L C}\, \mathit{vC0}\right) \sqrt{R^{2} C^{2}-4 L C}\, \sqrt{C \left(R^{2} C-4 L\right)}}{2 C \left(R^{2} C-4 L\right)}-\left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}+\left({\int}_{\,\,\,0}^{t}\mathit{vs}\, \left(\textit{\_}\mathit{z1}\right) {\mathrm e}^{\frac{\left(R C-\sqrt{C \left(R^{2} C-4 L\right)}\right) \textit{\_}\mathit{z1}}{2 L C}}{d}\textit{\_}\mathit{z1}\right) {\mathrm e}^{-\frac{R t}{L}+\frac{\left(R C+\sqrt{C \left(R^{2} C-4 L\right)}\right) t}{2 L C}}}{\sqrt{C \left(R^{2} C-4 L\right)}}\mathrm{ \\}\right\} \)

Numerical Version

Instead, what you will want to do is substitute in numerical values first, then get the solution. For second-order and higher systems, the types of responses can include exponentials, sinusoids, exponential sinusoids, and polynomials. Without knowing the relative values of the symbolic components, Maple cannot easily simplify. For this example, assume that we have a 1000 $\Omega$ resistor, a 100 nF capacitor, and a 100 mH inductor with a source that turns on to 5 V at time 0. Further assume that we know that the initial current in the inductor is 1 mA and the initial voltage across the capacitor is -2 V. We can get numerical versions of the equations and use them with the numerical versions of the initial conditions as follows: 

vals := R = 1000, C = 0.100*10^(-6), L = 0.100, vs(t) = 5;
numdeqns := subs(vals, deqns)[];
dsoln3 := dsolve([numdeqns, iL(0) = 0.001, vC(0) = -2], [iL(t), vC(t)])

You will get something less...complicated but still not small - click "Expand" at right to see it ->

\( \mathit{dsoln3}:= \left\{\mathit{iL}\! \left(t\right) = \frac{{\mathrm e}^{-5000 t} \left(2 \cos\! \left(5000 t \sqrt{3}\right)+\frac{26 \sqrt{3}\, \sin\left(5000 t \sqrt{3}\right)}{3}\right)}{2000}, \\ \mathit{vC}\! \left(t\right)=5+{\mathrm e}^{-5000 t} \left(-7 \cos\! \left(5000 t \sqrt{3}\right)-\frac{5 \sqrt{3}\, \sin\! \left(5000 t \sqrt{3}\right)}{3}\right)\right\} \)

You can look at a version where Maple does floating point evaluation and rounds things; to see that version, you can type: 

evalf[4](dsoln3)

and will see

\( \mathit{iL}\! \left(t\right) = 0.0005000 {\mathrm e}^{- 5000.0 t} \left( 2.0 \cos\! \left( 8660.0 t\right)+ 15.01 \sin\! \left( 8660.0 t\right)\right), \\ \mathit{vC}\! \left(t\right)= 5.0+{\mathrm e}^{- 5000.0 t} \left(- 7.0 \cos\! \left( 8660.0 t\right)- 2.887 \sin\! \left( 8660.0 t\right)\right) \)

Making Plots

Once you have the numerical solutions, you can make plots: 

plot(subs(dsoln3[], iL(t)), t = 0 .. 0.001);
plot(subs(dsoln3[], vC(t)), t = 0 .. 0.001)
Retrieved from "https://pundit.pratt.duke.edu/piki/index.php?title=Maple/Differential_Equations&oldid=30846"