Difference between revisions of "Python:Ordinary Differential Equations/Examples"

From PrattWiki
Jump to navigation Jump to search
 
 
(7 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
The following examples show different ways of setting up and solving initial value problems in Python.  It is part of the page on [[Python:Ordinary Differential Equations| Ordinary Differential Equations in Python]] and is very much based on [[MATLAB:Ordinary Differential Equations/Examples]].
 
The following examples show different ways of setting up and solving initial value problems in Python.  It is part of the page on [[Python:Ordinary Differential Equations| Ordinary Differential Equations in Python]] and is very much based on [[MATLAB:Ordinary Differential Equations/Examples]].
 +
 +
== Preamble ==
 +
The examples below assume a file called <code>ode_helpers.py</code> that contains the code below is in the same folder as the example codes; for the moment, this code contains a function that makes it easier to plot all the different dependent variables from a solver.
 +
<syntaxhighlight lang='python'>
 +
import numpy as np
 +
import matplotlib.pyplot as plt
 +
 +
def state_plotter(times, states, fig_num):
 +
    num_states = np.shape(states)[0]
 +
    num_cols = int(np.ceil(np.sqrt(num_states)))
 +
    num_rows = int(np.ceil(num_states / num_cols))
 +
    plt.figure(fig_num)
 +
    plt.clf()
 +
    fig, ax = plt.subplots(num_rows, num_cols, num=fig_num, clear=True,
 +
                        squeeze=False)
 +
    for n in range(num_states):
 +
        row = n // num_cols
 +
        col = n % num_cols
 +
        ax[row][col].plot(times, states[n], 'k.:')
 +
        ax[row][col].set(xlabel='Time',
 +
                        ylabel='$y_{:0.0f}(t)$'.format(n),
 +
                        title='$y_{:0.0f}(t)$ vs. Time'.format(n))
 +
       
 +
    for n in range(num_states, num_rows * num_cols):
 +
        fig.delaxes(ax[n // num_cols][n % num_cols])
 +
 +
    fig.tight_layout()
 +
 +
    return fig, ax
 +
</syntaxhighlight>
  
 
== Examples ==
 
== Examples ==
Note - each example began with the [[MATLAB:Ordinary Differential Equations/Templates|Templates]] provided at this web site.  Some comments have been removed from the templates to conserve space while some comments may have been added to provide a clearer explanation of the process for a particular example.
+
Note - each example began with the [[Python:Ordinary Differential Equations/Templates|Templates]] provided at this web site.  Some comments may have been removed from the templates to conserve space while some comments may have been added to provide a clearer explanation of the process for a particular example. The highlighted lines are the only lines that change between examples!
 +
 
 
=== Constant Rate of Change ===
 
=== Constant Rate of Change ===
[[File:ODEConstDiffPlot.png|thumb|Result using constant rate of change.]]
+
[[File:ODEConstDiffPlot_p.png|thumb|Result using constant rate of change.]]
 
If the dependent variable has a constant rate of change:
 
If the dependent variable has a constant rate of change:
 
<center><math>
 
<center><math>
Line 11: Line 42:
 
</math></center>
 
</math></center>
 
where <math>C</math> is some constant, you can provide the differential equation  
 
where <math>C</math> is some constant, you can provide the differential equation  
with a function called <code>ConstDiff.m</code> that contains the code:
+
in the <code>f</code> function and then calculate answers using this model with the code below. 
<source lang="matlab">
+
The code assumes there are 100 evenly spaced times between 0 and 10, the
function dydt = ConstDiff(t, y, C)
 
% Differential equation for constant growth
 
% t is time
 
% y is the state vector
 
% C contains any required constants
 
% dydt must be a column vector
 
dydt = C(1); % or just C since there is only one
 
</source>
 
 
 
You could calculate answers using this model with the following code
 
called <code>RunConstDiff.m</code>,
 
which assumes there are 100 evenly spaced times between 0 and 10, the
 
 
initial value of <math>y</math> is 6, and the rate of change is 1.2:
 
initial value of <math>y</math> is 6, and the rate of change is 1.2:
<source lang="matlab">
+
<syntaxhighlight lang="python" line highlight="9,13-15">
% Initialize workspace and graph
+
# %% Imports
clear; format short e; figure(1); clf
+
import numpy as np
 
+
import matplotlib.pyplot as plt
% Set name of file containing derivatives
+
from scipy.integrate import solve_ivp
DiffFileName = 'ConstDiff';
+
from ode_helpers import state_plotter
  
% Set up time span, initial value(s), and constant(s)
+
# %% Define independent function and derivative function
% Note: Variables should be in columns
+
def f(t, y, c):
tspan = linspace(0, 10);
+
    dydt = [c[0]]
yinit = 6;
+
    return dydt
C = 1.2;
 
  
% Determine if states should be plotted
+
# %% Define time spans, initial values, and constants
PlotStates = 1;
+
tspan = np.linspace(0, 10, 100)
 +
yinit = [6]
 +
c = [1.2]
  
%% Under the hood
+
# %% Solve differential equation
% Use ODE function of choice to get output times and states
+
sol = solve_ivp(lambda t, y: f(t, y, c),  
DE = eval(sprintf('@(t, y, C) %s(t,y,C)', DiffFileName))
+
                [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
[tout, yout] = ode45(@(t,y) DE(t,y,C), tspan, yinit);
 
  
% Plot results
+
# %% Plot states
if PlotStates
+
state_plotter(sol.t, sol.y, 1)
    figure(1); clf
+
</syntaxhighlight >
    StatePlotter(tout, yout)
 
end
 
</source>
 
  
 
===Time-dependent Rate of Change===
 
===Time-dependent Rate of Change===
[[File:ODETimeDiffPlot.png|thumb|Result using time-varying rate of change]]
+
[[File:ODETimeDiffPlot_p.png|thumb|Result using time-varying rate of change]]
 
If the dependent variable's rate of change is some function of time,
 
If the dependent variable's rate of change is some function of time,
this can be easily written using MATLAB.  For example, if the
+
this can be easily coded.  For example, if the
 
differential equation is some quadratic function given as:
 
differential equation is some quadratic function given as:
 
<center><math>
 
<center><math>
Line 65: Line 81:
 
</math></center>
 
</math></center>
 
then the function providing the values of the derivative may be
 
then the function providing the values of the derivative may be
written in a file called <code>TimeDiff.m</code>
+
written using <code>np.polyval</code>.
<source lang="matlab">
+
You could calculate answers using this model with the following code;
function dydt = TimeDiff(t, y, C)
+
it assumes there are 20 evenly spaced times between 0 and 4, the
% Differential equation for time-based polynomial derivative
 
% t is time
 
% y is the state vector
 
% C contains any required constants
 
% dydt must be a column vector
 
dydt = polyval(C, t);
 
</source>
 
 
 
You could calculate answers using this model with the following code
 
called <code>RunTimeDiff.m</code>,
 
which assumes there are 20 evenly spaced times between 0 and 4, the
 
 
initial value of <math>y</math> is 6, and the polynomial is defined by the vector
 
initial value of <math>y</math> is 6, and the polynomial is defined by the vector
[2 -6 3]:
+
[2, -6, 3]:
<source lang="matlab">
+
<syntaxhighlight lang="python" line highlight="9,13-15">
% Initialize workspace and graph
+
# %% Imports
clear; format short e; figure(1); clf
+
import numpy as np
 
+
import matplotlib.pyplot as plt
% Set name of file containing derivatives
+
from scipy.integrate import solve_ivp
DiffFileName = 'TimeDiff';
+
from ode_helpers import state_plotter
 
 
% Set up time span, initial value(s), and constant(s)
 
% Note: Variables should be in columns
 
tspan = linspace(0, 4, 20);
 
yinit = 6;
 
C    = [2 -6 3];
 
  
% Determine if states should be plotted
+
# %% Define derivative function
PlotStates = 1;
+
def f(t, y, c):
 +
    dydt = np.polyval(c, t)
 +
    return dydt
 +
   
 +
# %% Define time spans, initial values, and constants
 +
tspan = np.linspace(0, 4, 20)
 +
yinit = [6]
 +
c = [2, -6, 3]
  
%% Under the hood
+
# %% Solve differential equation
% Use ODE function of choice to get output times and states
+
sol = solve_ivp(lambda t, y: f(t, y, c),  
DE = eval(sprintf('@(t, y, C) %s(t,y,C)', DiffFileName))
+
                [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
[tout, yout] = ode45(@(t,y) DE(t,y,C), tspan, yinit);
 
  
% Plot results
+
# %% Plot states
if PlotStates
+
state_plotter(sol.t, sol.y, 1)
    StatePlotter(tout, yout)
+
</syntaxhighlight >
end
 
</source>
 
  
 
===Population Growth===
 
===Population Growth===
[[File:ODEPopDiffPlot.png|thumb|Result using rate of change proportional to measurement]]
+
[[File:ODEPopDiffPlot_p.png|thumb|Result using rate of change proportional to measurement]]
 
For population growth, the rate of change of population is dependent
 
For population growth, the rate of change of population is dependent
 
upon the number of people as well as some constant of
 
upon the number of people as well as some constant of
Line 119: Line 122:
 
</math></center>
 
</math></center>
 
where <math>C</math> is again some constant.
 
where <math>C</math> is again some constant.
In that case, the function may be written in a file called <code>PopDiff.m</code> as follows:
+
The following code will calculate the population for a span of 3 seconds with 25
<source lang="matlab">
 
function dydt = PopDiff(t, y, C)
 
% Differential equation for population growth
 
% t is time
 
% y is the state vector
 
% C contains any required constants
 
% dydt must be a column vector
 
dydt = C(1)*y(1); % or just C*y since both are 1x1
 
</source>
 
 
 
The following code, <code>RunPopDiff.m</code>, will calculate the population for
 
a span of 3 seconds with 25
 
 
points for the population model above with an initial population of 10
 
points for the population model above with an initial population of 10
 
and a constant of proportionality of 1.02:
 
and a constant of proportionality of 1.02:
<source lang="matlab">
+
<syntaxhighlight lang="python" line highlight="9,13-15">
% Initialize workspace and graph
+
# %% Imports
clear; format short e; figure(1); clf
+
import numpy as np
 +
import matplotlib.pyplot as plt
 +
from scipy.integrate import solve_ivp
 +
from ode_helpers import state_plotter
  
% Set name of file containing derivatives
+
# %% Define derivative function
DiffFileName = 'PopDiff';
+
def f(t, y, c):
 +
    dydt = [c[0] * y[0]]
 +
    return dydt
 +
   
 +
# %% Define time spans, initial values, and constants
 +
tspan = np.linspace(0, 3, 25)
 +
yinit = [10]
 +
c = [1.02]
  
% Set up time span, initial value(s), and constant(s)
+
# %% Solve differential equation
% Note: Variables should be in columns
+
sol = solve_ivp(lambda t, y: f(t, y, c),  
tspan = linspace(0, 3, 25);
+
                [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
yinit = 10;
 
C    = 1.02;
 
  
% Determine if states should be plotted
+
# %% Plot states
PlotStates = 1;
+
state_plotter(sol.t, sol.y, 1)
 
+
</syntaxhighlight>
%% Under the hood
 
% Use ODE function of choice to get output times and states
 
DE = eval(sprintf('@(t, y, C) %s(t,y,C)', DiffFileName))
 
[tout, yout] = ode45(@(t,y) DE(t,y,C), tspan, yinit);
 
 
 
% Plot results
 
if PlotStates
 
    StatePlotter(tout, yout)
 
end
 
</source>
 
  
  
 
===Multiple Variable Models===
 
===Multiple Variable Models===
[[File:ODETwoDiffPlot.png|thumb|Result for system with two variables]]
+
[[File:ODETwoDiffPlot_p.png|thumb|Result for system with two variables]]
  
 
It is possible to solve multiple-variable systems by making sure the
 
It is possible to solve multiple-variable systems by making sure the
Line 172: Line 161:
 
<center><math>
 
<center><math>
 
\begin{align}
 
\begin{align}
\frac{dy_1}{dt}&=\alpha\cos(\beta t) &
+
\frac{dy_0}{dt}&=\alpha\cos(\beta t) &
\frac{dy_2}{dt}&=\gamma y_1+\delta t
+
\frac{dy_1}{dt}&=\gamma y_0+\delta t
 
\end{align}
 
\end{align}
 
</math></center>
 
</math></center>
The differential file <code>TwoDiff.m</code> for this system will have a 2x1
+
The differential function <code>f</code> for this system will have a 2 element list as the output.
matrix as the output:
+
Also, if you have systems with multiple dependent variables, just
<source lang="matlab">
+
be sure to put the initial conditions in a list.  For
function dydt = TwoDiff(t, y, C)
 
% Differential equations for two variables
 
% t is time
 
% y is the state vector
 
% C contains any required constants
 
% dydt must be a column vector
 
dydt = [...
 
    C(1)*cos(C(2)*t);...
 
    C(3)*y(1)+C(4)*t];
 
</source>
 
 
 
Finally, if you have systems with multiple dependent variables, just
 
be sure to put the initial conditions in a column vector.  For
 
 
example, with the system defined as:
 
example, with the system defined as:
 
<center><math>
 
<center><math>
 
\begin{align}
 
\begin{align}
\frac{dy_1}{dt}&=4\cos(3t) &
+
\frac{dy_0}{dt}&=4\cos(3t) &
\frac{dy_2}{dt}&=-2y_1+0.5t
+
\frac{dy_1}{dt}&=-2y_0+0.5t
 
\end{align}
 
\end{align}
 
</math></center>
 
</math></center>
you could use the following script called <code>RunTwoDiff</code> to solve for both
+
you could use the following script to solve for both
<math>y_1</math> and <math>y_2</math> assuming <math>y_1</math> starts as 0 and <math>y_2</math> starts at -3:
+
<math>y_0</math> and <math>y_1</math>; the code assumes <math>y_0</math> starts as 0 and <math>y_1</math> starts at -3:
<source lang="matlab">
+
<syntaxhighlight lang="python" line highlight="9,13-15">
% Initialize workspace and graph
+
# %% Imports
clear; format short e; figure(1); clf
+
import numpy as np
 +
import matplotlib.pyplot as plt
 +
from scipy.integrate import solve_ivp
 +
from ode_helpers import state_plotter
  
% Set name of file containing derivatives
+
# %% Define derivative function
DiffFileName = 'TwoDiff';
+
def f(t, y, c):
 +
    dydt = [c[0]*np.cos(c[1]*t), c[2]*y[0]+c[3]*t]
 +
    return dydt
  
% Set up time span, initial value(s), and constant(s)
+
# %% Define time spans, initial values, and constants
% Note: Variables should be in columns
+
tspan = np.linspace(0, 5, 100)
tspan = linspace(0, 5);
+
yinit = [0, -3]
yinit = [0 -3]';
+
c = [4, 3, -2, 0.5]
C    = [4 3 -2 0.5];
 
  
% Determine if states should be plotted
+
# %% Solve differential equation
PlotStates = 1;
+
sol = solve_ivp(lambda t, y: f(t, y, c),
 +
                [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
  
%% Under the hood
+
# %% Plot states
% Use ODE function of choice to get output times and states
+
state_plotter(sol.t, sol.y, 1)
DE = eval(sprintf('@(t, y, C) %s(t,y,C)', DiffFileName))
+
</syntaxhighlight>
[tout, yout] = ode45(@(t,y) DE(t,y,C), tspan, yinit);
 
 
 
% Plot results
 
if PlotStates
 
    StatePlotter(tout, yout)
 
end
 
 
 
</source>
 
  
 
=== Higher Order Differential Equations ===
 
=== Higher Order Differential Equations ===
[[File:ODEJerkDiffPlot.png|thumb|Result using constant third derivative.]]
+
[[File:ODEJerkDiffPlot_p.png|thumb|Result using constant third derivative.]]
 
The system must be written in terms of first-order  
 
The system must be written in terms of first-order  
 
differential equations only.  To solve a system with
 
differential equations only.  To solve a system with
 
higher-order derivatives, you will first write a cascading system of simple
 
higher-order derivatives, you will first write a cascading system of simple
first-order equations then use them in your differential file.  For
+
first-order equations then use them in your differential function.  For
 
example, assume you have a system characterized by constant jerk:
 
example, assume you have a system characterized by constant jerk:
 
<center><math>
 
<center><math>
Line 245: Line 218:
 
<center><math>
 
<center><math>
 
\begin{align}
 
\begin{align}
y_1 &=y &~& &~\\
+
y_0 &=y &~& &~\\
\frac{dy_1}{dt}&=\frac{dy}{dt}=y_2 & &\longrightarrow & \frac{dy_1}{dt}&=y_2\\
+
\frac{dy_0}{dt}&=\frac{dy}{dt}=y_1 & &\longrightarrow & \frac{dy_0}{dt}&=y_1\\
\frac{dy_2}{dt}&=\frac{d^2y_1}{dt^2}=\frac{d^2y}{dt^2}=y_3 & &\longrightarrow &   
+
\frac{dy_1}{dt}&=\frac{d^2y_0}{dt^2}=\frac{d^2y}{dt^2}=y_2 & &\longrightarrow &   
\frac{dy_2}{dt}&=y_3\\
+
\frac{dy_1}{dt}&=y_2\\
\frac{dy_3}{dt}&=\frac{d^2y_2}{dt^2}=\frac{d^3y_1}{dt^3}=\frac{d^3y}{dt^3}=j=C &
+
\frac{dy_2}{dt}&=\frac{d^2y_1}{dt^2}=\frac{d^3y_0}{dt^3}=\frac{d^3y}{dt^3}=j=C &
 
&\longrightarrow &   
 
&\longrightarrow &   
\frac{dy_3}{dt}&=C
+
\frac{dy_2}{dt}&=C
 
\end{align}</math></center>
 
\end{align}</math></center>
 
Notice how the derivatives cascade so that the constant jerk equation
 
Notice how the derivatives cascade so that the constant jerk equation
can now be written as a set of three first-order equations.  The
+
can now be written as a set of three first-order equations.  Note that in this system,
differential file <code>JerkDiff.m</code> would thus be:
+
<code>y[0]</code> represents the position, <code>y[1]</code> represents the velocity, and
<source lang=matlab>
+
<code>y[2]</code> represents the acceleration. This type of cascading system will
function dydt = JerkDiff(t, y, C)
 
% Differential equations for constant jerk
 
% t is time
 
% y is the state vector
 
% C contains any required constants
 
% dydt must be a column vector
 
dydt = [...
 
    y(2); ...
 
    y(3);...
 
    C(1)]; % or just C since there is only one
 
</source>
 
to represent the three equations given above. Note that in this system,
 
<math>y_1</math> represents the position, <math>y_2</math> represents the velocity, and
 
<math>y_3</math> represents the acceleration. This type of cascading system will
 
 
show up often when modeling equations of motion.
 
show up often when modeling equations of motion.
  
Line 277: Line 236:
 
position of 6, and initial velocity of 2, an initial acceleration of
 
position of 6, and initial velocity of 2, an initial acceleration of
 
-4, and a constant jerk of 1.3:
 
-4, and a constant jerk of 1.3:
<source lang=matlab>
+
<syntaxhighlight lang='python' line highlight="9,13-15">
% Set name of file containing derivatives
+
# %% Imports
DiffFileName = 'JerkDiff';
+
import numpy as np
 
+
import matplotlib.pyplot as plt
% Set up time span, initial value(s), and constant(s)
+
from scipy.integrate import solve_ivp
% Note: Variables should be in columns
+
from ode_helpers import state_plotter
tspan = linspace(0, 8, 50);
 
yinit = [6 2 -4];
 
C    = 1.3;
 
 
 
% Determine if states should be plotted
 
PlotStates = 1;
 
  
%% Under the hood
+
# %% Define derivative function
% Use ODE function of choice to get output times and states
+
def f(t, y, c):
DE = eval(sprintf('@(t, y, C) %s(t,y,C)', DiffFileName))
+
    dydt = [y[1], y[2], c[0]]
[tout, yout] = ode45(@(t,y) DE(t,y,C), tspan, yinit);
+
    return dydt
  
% Plot results
+
# %% Define time spans, initial values, and constants
if PlotStates
+
tspan = np.linspace(0, 8, 50)
    StatePlotter(tout, yout)
+
yinit = [6, 2, -4]
end
+
c = [1.3]
</source>
 
Note in the above that the <code>StatePlotter</code> function is called to plot the states. That code is:
 
<source lang=matlab>
 
function StatePlotter(Time, States)
 
  
StateCount = size(States, 2);
+
# %% Solve differential equation
 +
sol = solve_ivp(lambda t, y: f(t, y, c),
 +
                [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
  
NumCols = ceil(sqrt(StateCount));
+
# %% Plot states
NumRows = ceil(StateCount / NumCols);
+
state_plotter(sol.t, sol.y, 1)
clf;
+
</syntaxhighlight>
for PlotNumber = 1:StateCount
 
        subplot(NumRows, NumCols, PlotNumber);
 
        plot(Time, States(:,PlotNumber), 'ko:');
 
        xlabel('Time');
 
        ylabel(sprintf('y_{%0.0f}(t)', PlotNumber))
 
        title(sprintf('y_{%0.0f}(t) vs. Time', PlotNumber));
 
end
 
</source>
 
  
 
== Questions ==
 
== Questions ==

Latest revision as of 14:38, 14 April 2021

The following examples show different ways of setting up and solving initial value problems in Python. It is part of the page on Ordinary Differential Equations in Python and is very much based on MATLAB:Ordinary Differential Equations/Examples.

Preamble

The examples below assume a file called ode_helpers.py that contains the code below is in the same folder as the example codes; for the moment, this code contains a function that makes it easier to plot all the different dependent variables from a solver.

import numpy as np
import matplotlib.pyplot as plt

def state_plotter(times, states, fig_num):
    num_states = np.shape(states)[0]
    num_cols = int(np.ceil(np.sqrt(num_states)))
    num_rows = int(np.ceil(num_states / num_cols))
    plt.figure(fig_num)
    plt.clf()
    fig, ax = plt.subplots(num_rows, num_cols, num=fig_num, clear=True,
                         squeeze=False)
    for n in range(num_states):
        row = n // num_cols
        col = n % num_cols
        ax[row][col].plot(times, states[n], 'k.:')
        ax[row][col].set(xlabel='Time',
                         ylabel='$y_{:0.0f}(t)$'.format(n),
                         title='$y_{:0.0f}(t)$ vs. Time'.format(n))
        
    for n in range(num_states, num_rows * num_cols):
        fig.delaxes(ax[n // num_cols][n % num_cols])

    fig.tight_layout()

    return fig, ax

Examples

Note - each example began with the Templates provided at this web site. Some comments may have been removed from the templates to conserve space while some comments may have been added to provide a clearer explanation of the process for a particular example. The highlighted lines are the only lines that change between examples!

Constant Rate of Change

Result using constant rate of change.

If the dependent variable has a constant rate of change:

\( \begin{align} \frac{dy}{dt}=C\end{align} \)

where \(C\) is some constant, you can provide the differential equation in the f function and then calculate answers using this model with the code below. The code assumes there are 100 evenly spaced times between 0 and 10, the initial value of \(y\) is 6, and the rate of change is 1.2:

 1 # %% Imports
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 from scipy.integrate import solve_ivp
 5 from ode_helpers import state_plotter
 6 
 7 # %% Define independent function and derivative function
 8 def f(t, y, c):
 9     dydt = [c[0]]
10     return dydt
11 
12 # %% Define time spans, initial values, and constants
13 tspan = np.linspace(0, 10, 100)
14 yinit = [6]
15 c = [1.2]
16 
17 # %% Solve differential equation
18 sol = solve_ivp(lambda t, y: f(t, y, c), 
19                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
20 
21 # %% Plot states
22 state_plotter(sol.t, sol.y, 1)

Time-dependent Rate of Change

Result using time-varying rate of change

If the dependent variable's rate of change is some function of time, this can be easily coded. For example, if the differential equation is some quadratic function given as:

\( \begin{align} \frac{dy}{dt}&=\alpha t^2+\beta t+\gamma \end{align} \)

then the function providing the values of the derivative may be written using np.polyval. You could calculate answers using this model with the following code; it assumes there are 20 evenly spaced times between 0 and 4, the initial value of \(y\) is 6, and the polynomial is defined by the vector [2, -6, 3]:

 1 # %% Imports
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 from scipy.integrate import solve_ivp
 5 from ode_helpers import state_plotter
 6 
 7 # %% Define derivative function
 8 def f(t, y, c):
 9     dydt = np.polyval(c, t)
10     return dydt
11     
12 # %% Define time spans, initial values, and constants
13 tspan = np.linspace(0, 4, 20)
14 yinit = [6]
15 c = [2, -6, 3]
16 
17 # %% Solve differential equation
18 sol = solve_ivp(lambda t, y: f(t, y, c), 
19                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
20 
21 # %% Plot states
22 state_plotter(sol.t, sol.y, 1)

Population Growth

Result using rate of change proportional to measurement

For population growth, the rate of change of population is dependent upon the number of people as well as some constant of proportionality:

\( \begin{align} \frac{dy}{dt}=C\cdot y \end{align} \)

where \(C\) is again some constant. The following code will calculate the population for a span of 3 seconds with 25 points for the population model above with an initial population of 10 and a constant of proportionality of 1.02:

 1 # %% Imports
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 from scipy.integrate import solve_ivp
 5 from ode_helpers import state_plotter
 6 
 7 # %% Define derivative function
 8 def f(t, y, c):
 9     dydt = [c[0] * y[0]]
10     return dydt
11     
12 # %% Define time spans, initial values, and constants
13 tspan = np.linspace(0, 3, 25)
14 yinit = [10]
15 c = [1.02]
16 
17 # %% Solve differential equation
18 sol = solve_ivp(lambda t, y: f(t, y, c), 
19                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
20 
21 # %% Plot states
22 state_plotter(sol.t, sol.y, 1)


Multiple Variable Models

Result for system with two variables

It is possible to solve multiple-variable systems by making sure the differential function returns values for each of the variables. For instance, in the following system the first variable's rate of change depends only on time while the second is dependent upon both time and the first variable:

\( \begin{align} \frac{dy_0}{dt}&=\alpha\cos(\beta t) & \frac{dy_1}{dt}&=\gamma y_0+\delta t \end{align} \)

The differential function f for this system will have a 2 element list as the output. Also, if you have systems with multiple dependent variables, just be sure to put the initial conditions in a list. For example, with the system defined as:

\( \begin{align} \frac{dy_0}{dt}&=4\cos(3t) & \frac{dy_1}{dt}&=-2y_0+0.5t \end{align} \)

you could use the following script to solve for both \(y_0\) and \(y_1\); the code assumes \(y_0\) starts as 0 and \(y_1\) starts at -3:

 1 # %% Imports
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 from scipy.integrate import solve_ivp
 5 from ode_helpers import state_plotter
 6 
 7 # %% Define derivative function
 8 def f(t, y, c):
 9     dydt = [c[0]*np.cos(c[1]*t), c[2]*y[0]+c[3]*t]
10     return dydt
11 
12 # %% Define time spans, initial values, and constants
13 tspan = np.linspace(0, 5, 100)
14 yinit = [0, -3]
15 c = [4, 3, -2, 0.5]
16 
17 # %% Solve differential equation
18 sol = solve_ivp(lambda t, y: f(t, y, c), 
19                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
20 
21 # %% Plot states
22 state_plotter(sol.t, sol.y, 1)

Higher Order Differential Equations

Result using constant third derivative.

The system must be written in terms of first-order differential equations only. To solve a system with higher-order derivatives, you will first write a cascading system of simple first-order equations then use them in your differential function. For example, assume you have a system characterized by constant jerk:

\( \begin{align} j&=\frac{d^3y}{dt^3}=C \end{align} \)

The first thing to do is write three first-order differential equations to represent the third-order equation:

\( \begin{align} y_0 &=y &~& &~\\ \frac{dy_0}{dt}&=\frac{dy}{dt}=y_1 & &\longrightarrow & \frac{dy_0}{dt}&=y_1\\ \frac{dy_1}{dt}&=\frac{d^2y_0}{dt^2}=\frac{d^2y}{dt^2}=y_2 & &\longrightarrow & \frac{dy_1}{dt}&=y_2\\ \frac{dy_2}{dt}&=\frac{d^2y_1}{dt^2}=\frac{d^3y_0}{dt^3}=\frac{d^3y}{dt^3}=j=C & &\longrightarrow & \frac{dy_2}{dt}&=C \end{align}\)

Notice how the derivatives cascade so that the constant jerk equation can now be written as a set of three first-order equations. Note that in this system, y[0] represents the position, y[1] represents the velocity, and y[2] represents the acceleration. This type of cascading system will show up often when modeling equations of motion.

The following script, RunJerkDiff.m, calculates the position, velocity, and speed over a period of 8 seconds assuming an initial position of 6, and initial velocity of 2, an initial acceleration of -4, and a constant jerk of 1.3:

 1 # %% Imports
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 from scipy.integrate import solve_ivp
 5 from ode_helpers import state_plotter
 6 
 7 # %% Define derivative function
 8 def f(t, y, c):
 9     dydt = [y[1], y[2], c[0]]
10     return dydt
11 
12 # %% Define time spans, initial values, and constants
13 tspan = np.linspace(0, 8, 50)
14 yinit = [6, 2, -4]
15 c = [1.3]
16 
17 # %% Solve differential equation
18 sol = solve_ivp(lambda t, y: f(t, y, c), 
19                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
20 
21 # %% Plot states
22 state_plotter(sol.t, sol.y, 1)

Questions

Post your questions by editing the discussion page of this article. Edit the page, then scroll to the bottom and add a question by putting in the characters *{{Q}}, followed by your question and finally your signature (with four tildes, i.e. ~~~~). Using the {{Q}} will automatically put the page in the category of pages with questions - other editors hoping to help out can then go to that category page to see where the questions are. See the page for Template:Q for details and examples.

External Links

References