Solving Simultaneous Equations with Python

I own a very old fashion scientific calculator and it can’t solve any simultaneous equations like those new calculators (not even 2×2!). The situation goes worst when I try to do my Circuit Theory tutorial, in which I need to solve many simultaneous equations. Can’t I just concentrate on forming the equations and let someone to solve them for me?

Well! Let the Python do it!

All you have to do is forming the proper equations. Then, fire up the following script. Key in the proper coefficients of each equation when Python asks.

Please feel free to copy and use it anywhere you like.


import numpy as np

print 'How many variables ?'
print '>> ',
total_variable = int(raw_input())

# User key in all the variable numbers

equation_list = []

for i in range(total_variable):
    print 'Please enter the coefficients of the #%d equation.' % (i+1)
    print 'For example: (1)*x0 + (2)*x1 = (3) ---> "1 2 3"'
    print '>> ',
    user_input = raw_input()
    user_token_str = user_input.split()
    assert (len(user_token_str)== (total_variable +1))
    user_token = [float(u) for u in user_token_str]

equation_arr = np.array(equation_list)

# A * X = Y
# Given A and Y, we need to find X.

Y_arr = equation_arr[:, -1:]
A_arr = equation_arr[:, :-1]

print 'Y_arr:' , Y_arr
print 'A_arr:' , A_arr

X_arr = np.linalg.solve(A_arr, Y_arr)

print 'Answer:'
for i, x in enumerate (X_arr):
    print 'x%d value: %f' % (i, x)

Solving Math’s Assignment with Sage

Someone sent me this question:

Solve for the currents in the circult of Figure 2, if E(t)=5H(t-2) and the initial currents are zero. [Hint : Use Lapalce transform to solve this problem.]

So, to solve it, form mesh analysis of two loops. Then, convert them from time domain to complex domain with Laplace transform. Next, solve I1 and I2 with normal algebra. Then only inverse I1 and I2 back to time domain.

Of cause, if you familiar with Sage, you can solve it within 30min (or lesser?).

t = var('t')
s = var('s')
I1 = var('I1')
I2 = var('I2')

E(t) = 5*unit_step(t-2)

E(s) = E(t).laplace(t, s); E(s)
# >> 5*e^(-2*s)/s

equation = [
	-E(s) + I1*20*s + 10*(I1-I2) == 0,
	10*(I2-I1) + I2*30*s + I2*10 == 0 ]

solution = solve(equation, I1, I2); solution
# >> [[I1 == 1/2*(3*s + 2)*e^(-2*s)/(6*s^3 + 7*s^2 + s), I2 == 1/2*e^(-2*s)/(6*s^3 + 7*s^2 + s)]]

# Note that Sage cannot inverse-Laplace time-delay function. So, taking out e^(-2*s)
I1(s) = 1/2*(3*s + 2)/(6*s^3 + 7*s^2 + s)
I2(s) = 1/2/(6*s^3 + 7*s^2 + s)

i1_temp(t) = I1(s).inverse_laplace(s, t); i1_temp
# >>t |--> -1/10*e^(-t) - 9/10*e^(-1/6*t) + 1

i2_temp(t) = I2(s).inverse_laplace(s, t); i2_temp
# >> t |--> 1/10*e^(-t) - 3/5*e^(-1/6*t) + 1/2

# Referring to Table. For G(s)= e^(as)F(s), the inverse is g(t) = f(t-a).
u(t) = unit_step(t)
i1(t) = u(t-2) * ( -1/10*e^(-(t-2)) - 9/10*e^(-1/6*(t-2)) + 1 ) # Answer for i1
i2(t) = u(t-2) * ( 1/10*e^(-(t-2)) - 3/5*e^(-1/6*(t-2)) + 1/2 ) # Answer for i2

p1 = plot(i1(t), 0, 10, color='blue', legend_label='i1(t)')
p2 = plot(i2(t), 0, 10, color='red', legend_label='i2(t)')
show(p1 + p2)

And, the final answers are:

i_1(t) =u(t)\left( \frac{-1}{10}e^{-(t-2)} -\frac{9}{10}e^{-(t-2)/6} +1 \right)
i_2(t) =u(t)\left( \frac{1}{10}e^{-(t-2)} -\frac{3}{5}e^{-(t-2)/6} +\frac{1}{2} \right)

Current I1 and I2 vs time.

Bye SAINT2012

I happened to be the last batch of speakers in SAINT2012 conference. Next year it will be merged with COMPSAC as “New COMPSAC”, which will be held in Kyoto in 2013. The name of “SAINT” will not be used anymore.

Badge in COMPSAC/SAINT2012

Happened to know that I have to teeeaaaach in the coming short-semester. It seems like I am getting far far away from my PhD completion. I blame noone but myself. When will I have the gut to tender resignation letter? When will I have the time to do nothing but my own research?

Ephesus, Turkey.

Visit some old cities of Turkey during the SAINT2012 conference. Tortured by the sunshine, and was amazed at the legacy of Ephesus.