So my previous post wasn't the last post of 2012. The streets are full of snow and the cold is keeping me indoors with nothing else to do but code. A couple of changes to the previous solver program:
The changes are a couple of bug fixes. In the function "mat_equation"
matrix_x.data[c1][0]=0.0
has been added as the result vector x has to be initialized to zero if the input is a vector. Similarly, in function "mat_ode"
ddt_mat.data[c1][0]=0.0
has been added to initialize value of dx/dt if input is a vector.
A major change has been to be able to solve the ODE when matrix E is singular. In the case, for an nXn matrix E, the rank will be k<n. This means that after the row operations, the lowest (n-k) rows will be zero rows. But if the diagonal elements of matrix A of these rows are not zero, the variable could still be calculated as that what this code does:
The concept here is that through the matrix A, these variables are still linked to the dynamic variables and so even if they are not directly affected by an input signal, they should be non-zero.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #! /usr/bin/env python | |
| import sys, math, matrix | |
| def row_op(matrix_a, matrix_b): | |
| """Function reduces a set of a set of equations Ax=b to upper triangular form.""" | |
| # Check if it is a genuine set of equations | |
| if not (matrix_a.rows==matrix_b.rows): | |
| print "Matrices do not have the same number of rows." | |
| else: | |
| for c1 in range(matrix_a.rows): | |
| #Check if the diagonal element is zero. | |
| if not matrix_a.data[c1][c1]: | |
| #If diagonal element is zero try to interchange it with a lower row. | |
| for c2 in range(c1+1, matrix_a.rows): | |
| if matrix_a.data[c2][c1]: | |
| for c3 in range(matrix_a.columns): | |
| matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3] | |
| for c3 in range(matrix_b.columns): | |
| matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3] | |
| #Make the diagonal element = 1. | |
| div_row=matrix_a.data[c1][c1] | |
| #Even after all possible row manipulations, | |
| #a diagonal element is zero, it means rank<size of matrix. | |
| if not div_row: | |
| print "Infinite solutions possible for this set of equations." | |
| else: | |
| for c2 in range(matrix_a.columns): | |
| matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row | |
| for c2 in range(matrix_b.columns): | |
| matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row | |
| #Row operations to remove row entries below the diagonal unity element. | |
| for c2 in range(c1+1, matrix_a.rows): | |
| mul_column=matrix_a.data[c2][c1] | |
| for c3 in range(matrix_a.columns): | |
| matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column | |
| for c3 in range(matrix_b.columns): | |
| matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column | |
| return | |
| def mat_equation(matrix_a, matrix_b, inp_signal, matrix_x): | |
| """ Function returns result x for equations Ax=bu """ | |
| row_op(matrix_a, matrix_b) | |
| #Calculate result matrix in reverse. | |
| for c1 in range(matrix_a.rows-1, -1, -1): | |
| if matrix_a.data[c1][c1]: | |
| try: | |
| if input_vector.rows: | |
| pass | |
| except: | |
| if not (matrix_b.columns==1): | |
| print "Input signal has to be a real number and not a vector." | |
| else: | |
| matrix_x.data[c1][0]=matrix_b.data[c1][0]*input_vector | |
| else: | |
| if not (matrix_b.columns==input_signal.rows): | |
| print "Dimension of input vector incorrect." | |
| else: | |
| matrix_x.data[c1][0]=0.0 # Added on Dec. 29, 2012. | |
| for c2 in range(matrix_b.columns): | |
| matrix_x.data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0] | |
| for c3 in range(c1+1, matrix_a.columns): | |
| matrix_x.data[c1][0]-=matrix_a.data[c1][c3]*matrix_x.data[c3][0] | |
| else: | |
| #If rank<size, set the (size-rank) solutions to zero. | |
| matrix_x.data[c1][0]=0.0 | |
| return | |
| def mat_ode_reduce(matrix_e, matrix_a, matrix_b): | |
| """ Function to modify the matrices for E d/dt(x)=Ax+bu. | |
| Upper triangularization is used to calculate the value of d/dt(x).""" | |
| # Check if it is a genuine set of equations | |
| if not ((matrix_e.rows==matrix_a.rows) and (matrix_a.rows==matrix_b.rows)): | |
| print "Matrices do not have the same number of rows." | |
| else: | |
| for c1 in range(matrix_e.rows): | |
| #Check if the diagonal element is zero. | |
| if not matrix_e.data[c1][c1]: | |
| #If diagonal element is zero try to interchange it with a lower row. | |
| for c2 in range(c1+1, matrix_e.rows): | |
| if matrix_e.data[c2][c1]: | |
| for c3 in range(matrix_e.columns): | |
| matrix_e.data[c1][c3], matrix_e.data[c2][c3] = matrix_e.data[c2][c3], matrix_e.data[c1][c3] | |
| for c3 in range(matrix_a.columns): | |
| matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3] | |
| for c3 in range(matrix_b.columns): | |
| matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3] | |
| #Make the diagonal element = 1. | |
| div_row=matrix_e.data[c1][c1] | |
| #Even after all possible row manipulations, | |
| #a diagonal element is zero, it means rank<size of matrix. | |
| if not div_row: | |
| # print "Check the circuit. Some variables have no dynamics." | |
| pass | |
| else: | |
| for c2 in range(matrix_e.columns): | |
| matrix_e.data[c1][c2]=matrix_e.data[c1][c2]/div_row | |
| for c2 in range(matrix_a.columns): | |
| matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row | |
| for c2 in range(matrix_b.columns): | |
| matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row | |
| #Row operations to remove row entries below the diagonal unity element. | |
| for c2 in range(c1+1, matrix_e.rows): | |
| mul_column=matrix_e.data[c2][c1] | |
| for c3 in range(matrix_e.columns): | |
| matrix_e.data[c2][c3]=matrix_e.data[c2][c3]-matrix_e.data[c1][c3]*mul_column | |
| for c3 in range(matrix_a.columns): | |
| matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column | |
| for c3 in range(matrix_b.columns): | |
| matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column | |
| return | |
| def mat_ode(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, runge_vectors): | |
| """ Solves the ODE using Runge Kutta 4th order method. """ | |
| mat_ode_reduce(matrix_e, matrix_a, matrix_b) | |
| def runge_function(x_plus_dxdt, ddt_mat): | |
| """ Defines the function dx/dt=f(x) """ | |
| #Calculate d/dt vector in reverse. | |
| for c1 in range(matrix_e.rows-1, -1, -1): | |
| if matrix_e.data[c1][c1]: | |
| try: | |
| if input_vector.rows: | |
| pass | |
| except: | |
| if not (matrix_b.columns==1): | |
| print "Input signal has to be a real number and not a vector." | |
| else: | |
| ddt_mat.data[c1][0]=matrix_b.data[c1][0]*input_vector | |
| else: | |
| if not (matrix_b.columns==input_vector.rows): | |
| print "Dimension of input vector incorrect." | |
| else: | |
| ddt_mat.data[c1][0]=0.0 # Added on Dec. 29, 2012. | |
| for c2 in range(matrix_b.columns): | |
| ddt_mat.data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0] | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat.data[c1][0]-=matrix_a.data[c1][c3]*x_plus_dxdt.data[c3][0] | |
| for c3 in range(c1+1, matrix_e.columns): | |
| ddt_mat.data[c1][0]-=matrix_e.data[c1][c3]*ddt_mat.data[c3][0] | |
| else: | |
| if not (matrix_a.data[c1][c1]): | |
| #The variable has no dynamics and can't even be calculated statically! | |
| print "This result will be wrong!" | |
| else: | |
| ddt_mat.data[c1][0]=0.0 | |
| state_vector[0].data[c1][0]=0.0 | |
| state_vector[1].data[c1][0]=0.0 | |
| for c2 in range(matrix_a.columns): | |
| if not (c1==c2): | |
| state_vector[0].data[c1][0]+=matrix_a.data[c1][c2]*state_vector[0].data[c2][0] | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] | |
| try: | |
| if input_vector.rows: | |
| pass | |
| except: | |
| if not (matrix_b.columns==1): | |
| print "Input signal has to be a real number and not a vector." | |
| else: | |
| state_vector[0].data[c1][0]+=matrix_b.data[c1][0]*input_vector | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] | |
| else: | |
| if not (matrix_b.columns==input_vector.rows): | |
| print "Dimension of input vector incorrect." | |
| else: | |
| for c2 in range(matrix_b.columns): | |
| state_vector[0].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0] | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] | |
| state_vector[0].data[c1][0]=-state_vector[0].data[c1][0]/matrix_a.data[c1][c1] | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] | |
| return ddt_mat | |
| # Runge-Kutta 4th order method | |
| k1=runge_vectors[0] | |
| k1(runge_function(state_vector[0], runge_vectors[4])*dt) #k1=dt*(dx/dt) | |
| k2=runge_vectors[1] | |
| k2(runge_function(state_vector[0]+0.5*k1, runge_vectors[4])*dt) #k2=dt*d/dt(x+k1/2) | |
| k3=runge_vectors[2] | |
| k3(runge_function(state_vector[0]+0.5*k2, runge_vectors[4])*dt) #k3=dt*d/dt(x+k2/2) | |
| k4=runge_vectors[3] | |
| k4(runge_function(state_vector[0]+k3, runge_vectors[4])*dt) #k4=dt*d/dt(x+k3) | |
| for c1 in range(state_vector[1].rows): | |
| state_vector[1].data[c1][0]=(state_vector[0].data[c1][0]+(k1.data[c1][0]+2.0*k2.data[c1][0]+2.0*k3.data[c1][0]+k4.data[c1][0])*(1/6.0)) | |
| return | |
| R1=10.0 | |
| R2=10.0 | |
| R3=10.0 | |
| L1=20.0e-3 | |
| L2=0.0e-3 | |
| L3=0.0e-3 | |
| a=matrix.Matrix(2,2) | |
| b=matrix.Matrix(2) | |
| b.data=[[1.0],[0.0]] | |
| a.data=[[R1+R2, -R2], [-R2, R2+R3]] | |
| e=matrix.Matrix(2,2) | |
| e.data=[[L1+L2, -L2], [-L2, L2+L3]] | |
| curr_state_vec = matrix.Matrix(2) | |
| curr_state_vec.zeros(2) | |
| next_state_vec = matrix.Matrix(2) | |
| next_state_vec.zeros(2) | |
| dt=20.0e-6 | |
| f=open("data_file.dat","w") | |
| input_volt=matrix.Matrix(2) | |
| ode_k1=matrix.Matrix(2) | |
| ode_k2=matrix.Matrix(2) | |
| ode_k3=matrix.Matrix(2) | |
| ode_k4=matrix.Matrix(2) | |
| ode_dbydt=matrix.Matrix(2) | |
| ode_solver=[ode_k1, ode_k2, ode_k3, ode_k4, ode_dbydt] | |
| for time_step in range(0,10000): | |
| t=time_step*dt | |
| u=100.0*math.sin(2*60*math.pi*t) | |
| mat_ode(e, a, b, [curr_state_vec, next_state_vec], u, 20.0e-6, ode_solver) | |
| curr_state_vec=next_state_vec | |
| f.write("%s \t %s \t %s \t %s \n" %(str(t), str(u), str(next_state_vec.data[0][0]), str(next_state_vec.data[1][0]))) | |
| print a | |
| print e | |
| print b | |
| f.close() |
The changes are a couple of bug fixes. In the function "mat_equation"
matrix_x.data[c1][0]=0.0
has been added as the result vector x has to be initialized to zero if the input is a vector. Similarly, in function "mat_ode"
ddt_mat.data[c1][0]=0.0
has been added to initialize value of dx/dt if input is a vector.
A major change has been to be able to solve the ODE when matrix E is singular. In the case, for an nXn matrix E, the rank will be k<n. This means that after the row operations, the lowest (n-k) rows will be zero rows. But if the diagonal elements of matrix A of these rows are not zero, the variable could still be calculated as that what this code does:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| if not (matrix_a.data[c1][c1]): | |
| #The variable has no dynamics and can't even be calculated statically! | |
| print "This result will be wrong!" | |
| else: | |
| ddt_mat.data[c1][0]=0.0 | |
| state_vector[0].data[c1][0]=0.0 | |
| state_vector[1].data[c1][0]=0.0 | |
| for c2 in range(matrix_a.columns): | |
| if not (c1==c2): | |
| state_vector[0].data[c1][0]+=matrix_a.data[c1][c2]*state_vector[0].data[c2][0] | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] | |
| try: | |
| if input_vector.rows: | |
| pass | |
| except: | |
| if not (matrix_b.columns==1): | |
| print "Input signal has to be a real number and not a vector." | |
| else: | |
| state_vector[0].data[c1][0]+=matrix_b.data[c1][0]*input_vector | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] | |
| else: | |
| if not (matrix_b.columns==input_vector.rows): | |
| print "Dimension of input vector incorrect." | |
| else: | |
| for c2 in range(matrix_b.columns): | |
| state_vector[0].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0] | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] | |
| state_vector[0].data[c1][0]=-state_vector[0].data[c1][0]/matrix_a.data[c1][c1] | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0] |
The concept here is that through the matrix A, these variables are still linked to the dynamic variables and so even if they are not directly affected by an input signal, they should be non-zero.
No comments:
Post a Comment