#!/usr/bin/env python
# coding: utf-8

# <style>
# @import url(https://www.numfys.net/static/css/nbstyle.css);
# </style>
# <a href="https://www.numfys.net"><img class="logo" /></a>
# 
# # Resistor Network
# 
# ### Examples - Electromagnetism
# <section class="post-meta">
# By Magnus A. Gjennestad, Vegard Hagen, Aksel Kvaal, Morten Vassvik, Trygve B. Wiig, Peter Berg and Magnus Dahle
# </section>
# Last edited: March 22nd 2018 
# ___
# 
# 
# ### Problem Formulation and Notation
# 
# Within this example, we'll consider a resistor network consisting of $N$ repetitive units. Each unit has two resistors of magnitude $R$ and $12R$, respectively. The units are connected to a battery of voltage $V_0$, as shown in the figure below.
# 
# 

# ![circuit1](images/4ec_image_1.png)

# The goal of this example is to calculate the toal current propagating through the cicuit, provided by the battery. The following image explains the applied notations and variables.

# ![circuit2](images/4ec_image_2.png)

# ### The Special Case: $N=1$

# When solving a problem like this, it is often useful to look at special scenarios, if available, before one tries to solve the general problem. We will now look at two such simplified problems.
# 
# First, consider the special case when $N=1$. The circuit then looks like the following,

# ![circuit2](images/4ec_image_3.png)

# This reduced problem is easily solved by defining $R_{eff}$, the effective resistance of the circuit connected to the battery. This implies that we may represent the circuit by the following diagram

# ![circuit2](images/4ec_image_4.png)

# where $R_{eff} = R+12R=13R$. Then, by Ohm's law,
# 
# $$
# I_{1,1} = \frac{V_0}{13R} = \frac{1}{13} \frac{\textrm{V}}{\Omega} \approx 0.0769 \textrm{A}
# $$

# ### The Special Case: $N\to\infty$

# As another special scenario, we will consider the case when the number of units $N$ goes to infinity. This scenario is not so trivial. Obviously, we are providing more options for the current to flow as $N\to\infty$. Hence, we expect the resistance to decrease in this limit. Take a minute to think about how you would solve it before you read on.

# Again, we can consider an effective resistance $R_{eff}$ and represent the whole circuit as before (last image). Now, take the circuit from the previuos page with infinitely many repetitive units and add one more unit to it, so that it can be represented by the following circuit diagram

# ![circuit2](images/4ec_image_5.png)

# Since $N$ is infinitely large, adding one more unit should not change the effective resistance $R_{eff}$ of the whole circuit. In other words, $12R$ and $R_{eff}$ are in parallel in the above sketch and the resistance of the entire circuit $R_{total}$ is still $R_{eff}$. By this argument, we can relate $R_{eff}$ to itself by
# 
# $$
# R_{total} = R_{eff} = R + \frac{1}{\frac{1}{12R} + \frac{1}{R_{eff}}}
# $$

# By solving this quadric equation for $R_{eff}$ and demanding that $R_{eff}$ be positive, we obtain
# 
# $$
# R_{eff} = 4R = 4\Omega
# $$

# Such that by Ohm's law, we find
# $$
# I_{1,1} = \frac{V_0}{4R} = \frac{1}{4} \frac{\textrm{V}}{\Omega} \approx 0.25 \textrm{A}
# $$

# ### Setting Up the System of Equations

# Now, we turn to the more general case when $1<N<\infty$. To solve it, we will set up a system of $N$ equations and $N$ unknowns that we will formulate as a matrix problem. It will then be solven using Python. The unknowns will be the $N$ voltages $V_i$, $i=1,\ldots,N;$.

# To obtain the $N$ equations and $N$ unknowns, we first apply Ohm's law to all resistors in the circuit. This yields
# 
# $$
# I_{i,1} = \frac{V_{i-1} - V_i}{R}, \quad i=1,\ldots,N;
# $$
# 
# for the $N$ resistors of resistance $R$ and
# 
# $$
# I_{i,2} = \frac{V_i}{12R}, \quad i=1,\ldots,N;
# $$
# 
# for the $N$ resistors with resistance 12R.

# The next step is to eliminate the currents $I_{i,1}$ and $I_{i,2}$ for $i=1,\ldots,N;$ in the equations above. To do so, we turn to the principle of conservation of current, namely that the sum of all currents flowing into a node in a circuit diagram is equal to the sum of all currents flowing out of it. This statement equates to saying that no charges are created or destroyed in a node and is often reffered to as Krichoff's current law.

# For the nodes labelled $V_i$, where $i=1,\ldots,N-1;$ we get
# 
# $$
# I_{i,1} = I_{i,2} + I_{i+1,1}
# $$
# 
# For the last node, labelled $V_N$, we get
# 
# $$
# I_{N,1} = I_{N,2}
# $$
# 

# Substituting the earliest expressions for $I_{i,1}$ and $I_{i,2}$ into the each of these two last equations separately (the upper one first), we find for the first node $i=1$,
# 
# $$
# \frac{25}{12R} V_1 - \frac{1}{R}V_2 = V_0
# $$
# 
# and for the nodes labelled $V_i$, with $i=2,\ldots,N-1;$
# 
# $$
# -\frac{1}{R}V_{i-1} + \frac{25}{12R}V_i - \frac{1}{R} V_{i+1} = 0
# $$
# 
# and then for the lower equation,
# 
# $$
# -\frac{1}{R}V_{N-1} + \frac{13}{12R}V_N = 0
# $$
# 
# Counting the number of equations, we see that these three last expressions contain $N$ equations in total. This is <i>exactly</i> the amount we need to determine all $N$ voltages $V_i$ for $i=2,\ldots,N-1;$ uniquely.

# Moreover, these equations can be formulated as a matrix problem $\mathcal{A}\boldsymbol{V}=\boldsymbol{b}$
# 
# $$
# \begin{bmatrix} 25/12R & -1/R  & 0    & \dots & 0 \\ 
#                 -1/R  & 25/12R & -1/R & \dots & 0 \\
#                 \vdots& \ddots &\ddots& \ddots& \vdots \\
#                 0     & \dots &-1/R  & 25/12R & -1/R  \\
#                 0     & \dots & 0    & -1/R   & 25/12R
# \end{bmatrix} \cdot 
# \begin{bmatrix} V_1   \\ V_2 \\ \vdots \\ V_{N-1} \\ V_N \end{bmatrix} =
# \begin{bmatrix} V_0/R \\ 0   \\ \vdots \\ 0       \\ 0   \end{bmatrix} 
# $$

# Now, finding the unknown voltages amounts to solving this matrix equation $\mathcal{A}\boldsymbol{V}=\boldsymbol{b}$, for the unknown voltage vector $\boldsymbol{V}$. The matrix $\mathcal{A}$ and the vector $\boldsymbol{b}$ are given by the resistances within the circuit and the voltage $V_0$. Subsequently, we can calculate the total current by Ohm's law,
# 
# $$
# I_{1,1} = \frac{V_0 - V_1}{R}
# $$
# 
# In Python, the first thing we need to do is defineing $V_0$, $R$ and $N$.

# In[1]:


R  = 1.0 # Resistance [Ohm]
V0 = 1.0 # Applied voltage [V]
N  = 10  # Number of repitative units [dimless]


# Then, we need to set up the matrix $\mathcal{A}$ and the vector $\boldsymbol{b}$ of the matrix equation. We start by initializing them both with zeros before filling them and also define two variables, $a=25/12R$ and $c=-1/R$ to simplify their filling.

# In[2]:


# We use numpy arrays for increased efficiency
# of matrix operations and functionality.
import numpy as np 

A = np.zeros((N,N)) # Matrix of dimension NxN
b = np.zeros(N)     # Vector of length N
a = 25.0/(12*R)     # Scalar (constant)
c = -1.0/R          # Scalar (constant)


# Next, we set up $\boldsymbol{b}$, and $\mathcal{A}$ row by row.

# In[3]:


# the b-vector is all zeros except for first entry. Thus,
b[0] = V0/R

# Set up first row
A[0,0] = a
A[0,1] = c 

# Set up last row
A[N-1,N-1] = a
A[N-1,N-2] = c

# Set up all other rows
# OBS: if you're running Python 2.7 (or lower) you might want to use
#      the 'xrange()' instead of 'range()' depending on the size of 'N'.
for row in range(1,N-1):
    A[row,row-1] = c
    A[row,row  ] = a
    A[row,row+1] = c

# You may want to print A,b to see if it was initialized correctly:
# print(A,b)


# Then we can solve the system of equations by using the built-in <i>Numerical Python Linear Algebra Solver</i>

# In[4]:


Voltages = np.linalg.solve(A,b)
print("Voltages = ", Voltages)

I11 = (V0 - Voltages[0])/R
print("\nI_11     = ", I11)


# We see that when $N$ becomes large (or actually already at $N\geq15$), $I_{1,1}$ approaches the limit we found analytically as $N\to\infty$. That is $I_{1,1}\to 1/4$.
# 
# You should note that, even though the built-in solver we use in this example is implemented to be effective, it is a general solver which do <strong>not</strong> take advantage of the fact that $\mathcal{A}$ is a sparse matrix for this problem.
# 
# If $\mathcal{A}$ becomes really large, then this solver would eventually become very slow as it unnessecarily iterates over all the zeros. As an alternative, one can use the built-in functionality of the <i>Scientific Python Sparse Linear Algebra</i> package. However, this module requires $\mathcal{A}$ to be stored in a specific (sparse) manner.

# In[5]:


from scipy import sparse
import scipy.sparse.linalg as ssl
from scipy.sparse.linalg import spsolve


# Since our matrix only have non-zero values on 
# THE diagonal, the UPPER(sup) diagonal and the LOWER(sub) diagonal,
# this is in fact all we need to tell Python.

#Create create a sparse matrix

sup_diag = np.ones(N)*c
sub_diag = np.ones(N)*c
the_diag = np.ones(N)*a
           
# Offset:           -1     0        1
all_diags  = [sub_diag, the_diag, sup_diag]
offsets    = np.array([-1, 0, 1])


csc="csc" # Computer format of which the matrix is stored

# Define the SPARSE matrix
A_sparse = sparse.spdiags(all_diags, offsets, N,N, format=csc)


# print(A_sparse.todense()) # prints the SPARSE matrix in a DENSE (normal NxN) format.


Voltages = spsolve(A_sparse,b)
print("Voltages = ", Voltages)

I11 = (V0 - Voltages[0])/R
print("\nI_11     = ", I11)



# Now, try to compare running times of the sparse solver, with the standard solver as you increase $N$ to $10, 100, 1000, \ldots$.