From Symbolics to Numerics

An engineering analysis typically requires that a symbolic solution be applied via the substitution of numbers into a symbolic expression. In subsection 4.2.7, we considered how to subsitute numerical values into expressions using SymPy’s evalf() method. This is fine for a single value, but frequently an expression is to be evaluated at an array of numerical values. Looping through the array and applying evalf() is cumbersome and computationally slow. An easier and computationally efficient technique using the sp.lambdify() function is introduced in this section. The function sp.lambdify() creates an efficient, numerically evaluable function from a SymPy expression. The basic usage of the function is as follows:

x = sp.symbols("x", real=True)
expr = x**2 + 7
f = sp.lambdify(x, expr)
f(2)

By default, if NumPy is present, sp.lambdify() vectorizes the function such that the function can be provided with NumPy array arguments and return NumPy array values. However, it is best to avoid relying on the function’s implicit behavior, which can change when different modules are present, it is best to provide the numerical module explicitly, as follows:

f = sp.lambdify(x, expr, modules="numpy")
f(np.array([1, 2, 3.5]))

array([ 8.  , 11.  , 19.25])

Multiple arguments are supported, as in the following example:

x, y = sp.symbols("x, y", real=True)
expr = sp.cos(x) * sp.sin(y)
f = sp.lambdify([x, y], expr, modules="numpy")
f(3, 4)

0.7492287917633428

All the usual NumPy broadcasting rules will apply for the function. For instance,

X = np.array([[1], [2]])  # 2x1 matrix
Y = np.array([[1, 2, 3]])  # 1x3 matrix
f(X, Y)

array([[ 0.45464871,  0.4912955 ,  0.07624747],
       [-0.35017549, -0.37840125, -0.05872664]])

Example 4.2

You are designing the circuit shown in [@fig:ur-circuit]. Treat the source voltage V_S, the source resistance R_S, and the overall circuit topology as known constants. The circuit design requires the selection of resistances R₁, R₂, and R₃ such that the voltage across R₃, v_R₃ = V_R₃, and the current through R₁, i_R₁ = I_R₁, where V_R₃ and I_R₁ are known constants (i.e., design requirements). Proceed through the following steps:

Solve for all the resistor voltages v_{R_k} and currents i_{R_k} in terms of known constants and R₁, R₂, and R₃ using circuit laws
Apply the constraints v_R₃ = V_R₃ and i_R₁ = I_R₁ to obtain two equations relating R₁, R₂, and R₃
Solve for R₂ and R₃ as functions of R₁ and known constants
Create a design graph for the selection of R₁, R₂, and R₃ given the following design parameters: V_S = 10 V, R_S = 50 Ω, V_R₃ = 1 V, and I_R₁ = 20 mA.

Figure 4.4: A resistor circuit design for example.

Solve for the Resistor Voltages and Currents

Each resistor has an unknown voltage and current. We will develop and solve a system of equations using circuit laws. Begin by defining symbolic variables as follows:

v_RS, i_RS, v_R1, i_R1, v_R2, i_R2, v_R3, i_R3 = sp.symbols(
	"v_RS, i_RS, v_R1, i_R1, v_R2, i_R2, v_R3, i_R3", real=True
)
viR_vars = [v_RS, i_RS, v_R1, i_R1, v_R2, i_R2, v_R3, i_R3]
R1, R2, R3 = sp.symbols("R1, R2, R3", positive=True)
V_S, R_S, V_R3, I_R1 = sp.symbols("V_S, R_S, V_R3, I_R1", real=True)

There are 4 resistors, so there are 2 ⋅ 4 = 8 unknown voltages and currents; therefore, we need 8 independent equations. The first circuit law we apply is Ohm’s law, which states that the ratio of voltage over current for a resistor is approximately constant. Applying this to each resistor, we obtain the following 4 equations:

Ohms_law = [
	v_RS - R_S*i_RS,  # == 0
	v_R1 - R1*i_R1,  # == 0
	v_R2 - R2*i_R2,  # == 0
	v_R3 - R3*i_R3,  # == 0
]

The second circuit law we apply is Kirchhoff’s current law (KCL), which states that the sum of the current into a node must equal 0. Applying this to the upper-middle and upper-right nodes, we obtain the following 2 equations:

KCL = [
	i_RS - i_R1 - i_R2,  # == 0
	i_R2 - i_R3,  # == 0
]

The third circuit law we apply is Kirchhoff’s voltage law (KVL), which states that the sum of the voltage around a closed loop must equal 0. Applying this to the left and right inner loops, we obtain the following 2 equations:

KVL = [
	V_S - v_R1 - v_RS,  # == 0
	v_R1 - v_R3 - v_R2, # == 0
]

Our collection of 8 equations are independent because none can be derived from another. They make a linear system of equations, which can be solved simultaneously as follows:

viR_sol = sp.solve(Ohms_law + KCL + KVL, viR_vars, dict=True)[0]
print(viR_sol)

Apply the Requirement Constraints

The requirements that v_R₃ = V_R₃ and i_R₁ = I_R₁ can be encoded symbolically as two equations as follows:

constraints = {v_R3: V_R3, i_R1: I_R1}  # Design constraints
constraint_equations = [
	sp.Eq(v_R3.subs(constraints), v_R3.subs(viR_sol)),
	sp.Eq(i_R1.subs(constraints), i_R1.subs(viR_sol)),
]
print(constraint_equations)

Solve for Resistances

The system of 2 constraint equations and 3 unkowns (R₁, R₂, and R₃) is underdetermined, which means there are infinite solutions. The two equations can be solved for R₁ and R₂ in terms of R₃ and parameters as follows:

constraints_sol = sp.solve(
	constraint_equations, [R1, R2], dict=True
)[0]
print(constraints_sol)

Create a Design Graph

Applying the design parameters and defining numerically evaluable functions for R₁ and R₂ as functions of R₃,

design_params = {V_S: 10, R_S: 50, V_R3: 1, I_R1: 0.02}
R1_fun = sp.lambdify(
	[R3],
	R1.subs(constraints_sol).subs(design_params),
	modules="numpy",
)
R2_fun = sp.lambdify(
	[R3],
	R2.subs(constraints_sol).subs(design_params),
	modules="numpy",
)

And now we are ready to create the design graph, as follows:

R3_ = np.linspace(10, 100, 101)  # Values of $R_3$
fig, ax = plt.subplots()
ax.plot(R3_, R1_fun(R3_), label="$R_1$ ($\\Omega$)")
ax.plot(R3_, R2_fun(R3_), label="$R_2$ ($\\Omega$)")
ax.set_xlabel("$R_3$ ($\\Omega$)")
ax.legend()
ax.grid()
plt.show()

Figure 4.5: A design graph for resistors R₁, R₂, and R₃.

Online Resources for Section 4.4

No online resources.