Live Coding Session: Modeling and Simulation Basics

Objectives

• Work with vectors and matrices in Python.

• Learn how to iterate through arrays and perform matrix operations.

• Explore symbolic math using SymPy.

• Create simple plots with Matplotlib.

• Define functions

---

Section 1: Arrays and Matrices

Step 1.1: Create a Row and Column Vector

We start by creating a row vector and a column vector using NumPy.

import numpy as np

# Make a row vector using Numpy
v_flat = np.array([1, 2, 3])
print("v_flat: ", v_flat)

# Make a column vector using Numpy
v_tall = np.array([[1], [2], [3]])
print("v_tall:\n ", v_tall)

v_flat:  [1 2 3]
v_tall:
  [[1]
 [2]
 [3]]

Step 1.2: Create a Matrix and perform matrix-vector multiplications

We can also use Numpy to defind matrices and perform various mathematical operations on matrices and vectors

# Create a matrix
M = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print("Matrix: \n", M)

# Multipy the matrix by a vector using Nympy's dot method
new_vector = np.dot(M, v_tall)
print("new vector:\n ", new_vector)

# Make a scalar by dotting the column vector, matrix and row vector
scalar = np.dot(np.dot(v_flat, M), v_tall)

print("scalar: ", scalar)

Matrix:
 [[1 2 3]
 [4 5 6]
 [7 8 9]]
new vector:
  [[14]
 [32]
 [50]]
scalar:  [228]

Section 2: Arrays

Step 2.1: Generating an array

Our task is to generate an array of time instances from 0 to 10 seconds with 0.5 second time steps.

The linspace method takes the following arguments:

• The starting value of the sequence

• The end value of hte sequence

• The number of samples to generate

# Generate an array using Numpy's linspace method
step_size = 0.5
n_steps = int(10/0.5)+1

time_steps = np.linspace(0, 10, n_steps)
print(time_steps)

[ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4.   4.5  5.   5.5  6.   6.5
  7.   7.5  8.   8.5  9.   9.5 10. ]

Step 2.2: Iterate through the array

We will use a for-loop to iterate through the array we generated to print each time step

# Iterate through the array and print something

for t in time_steps:
    print("time is: ", t)

time is:  0.0
time is:  0.5
time is:  1.0
time is:  1.5
time is:  2.0
time is:  2.5
time is:  3.0
time is:  3.5
time is:  4.0
time is:  4.5
time is:  5.0
time is:  5.5
time is:  6.0
time is:  6.5
time is:  7.0
time is:  7.5
time is:  8.0
time is:  8.5
time is:  9.0
time is:  9.5
time is:  10.0

Section 3: Symbolic math

Step 3.1: Defining variables and functions

With symbolic math tools, we can define variables and do matematical operations on them symbolically. For example we may define functions and differentiate them.

# Import diff, symbols and Matrix from sympy
from sympy import symbols, diff, Matrix, sin

# Define variables x1 and x2
x1, x2 = symbols("x1 x2")

# Define two functions f1 and f2
f1 = x1 ** 3 + 2*x2
f2 = x2 * sin(x1)

print(f1)

x1**3 + 2*x2

Step 3.2: Find Partial derivatives

Finding partial derivatives is something we have to deal with often in this course, and often the expressions are not nice to work with by hand.

# Find the partial derivatives of f1 and f2 wrt x1 and x2
df1_dx1 = diff(f1, x1)
df1_dx2 = diff(f1, x1)
df2_dx1 = diff(f1, x1)
df2_dx2 = diff(f1, x1)

print(df1_dx1)

3*x1**2

Step 3.3: Use the chain rule in sympy

Given the function \(f(x_1(t), x_2(t))\) as an example, we have that \(\frac{df}{dt} = \frac{\partial f}{\partial x_1} \frac{dx_1}{dt} + \frac{\partial f}{\partial x_2} \frac{dx_2}{dt}\)

Let’s pretend like x1 and x2, defined above, are functions of time, and then try to find the time derivative of our funcitons f1 and f2

# Define variables for time derivatives of x1 and x2
dx1, dx2 = symbols("dx1 dx2")


# Find the time derivative of f1 and f2
df1_dt = df1_dx1 * dx1 + df1_dx2 * dx2
print(df1_dt)

3*dx1*x1**2 + 3*dx2*x1**2

Step 3.4: Find the Jacobian for a Vector-Valued Function

Given a function \(\boldsymbol{f}(\boldsymbol{x})\), where \(\boldsymbol{f}\) is an n-dimensional vector of functions and \(\boldsymbol{x}\) is an m-dimensional vector, the Jacobian is defined as:

(1)\[\begin{split}\mathbf{J} = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \dots & \frac{\partial f_1}{\partial x_m} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \dots & \frac{\partial f_n}{\partial x_m} \end{bmatrix}\end{split}\]

We will find the Jacobian for \(f_1\) and \(f_2\).

# Define the Jacobian for f1 and f2
functions = Matrix([f1, f2])
variables = Matrix([x1, x2])

jacobian = functions.jacobian(variables)
print(jacobian)

Matrix([[3*x1**2, 2], [x2*cos(x1), sin(x1)]])

Section 4: Plotting

We will make a time array and plot teh function \(\sin t\) using the Python module Matplotlib

# Make a time array from 0 to 2*pi

# Make an array for the value of sin(t) for each value of the time array

# Set up the plot

Section 5: Defining functions

In Python (and MATLAB) you can define functions (methods). We will make an example here

# Define a function to calculate the square root of an input, including handling of negative numbers

# Call the function and show the result