With the changes in the Matrix class, another round of changes to the ODE solver. The Runge Kutta function has been primarily changed to remove all object constructors and object returns in the functions. The new solver functions is:
The full code is:
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
| 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, dbydt_order): | |
| """ 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[dbydt_order-1].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[dbydt_order-1].data[c1][0]=0.0 # Added on Dec. 29, 2012. | |
| for c2 in range(matrix_b.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0] | |
| if (dbydt_order==2): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0] | |
| if (dbydt_order==3): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0] | |
| if (dbydt_order==4): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*ddt_mat[dbydt_order-2].data[c3][0] | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].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[dbydt_order-1].data[c1][0]-=matrix_e.data[c1][c3]*runge_vectors[4].data[c3][0]/dt | |
| ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0]*dt | |
| 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[dbydt_order-1].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 | |
| # Runge-Kutta 4th order method | |
| runge_function(state_vector[0], runge_vectors, 1) #k1=dt*(dx/dt) | |
| runge_function(state_vector[0], runge_vectors, 2) #k2=dt*d/dt(x+k1/2) | |
| runge_function(state_vector[0], runge_vectors, 3) #k3=dt*d/dt(x+k2/2) | |
| runge_function(state_vector[0], runge_vectors, 4) #k4=dt*d/dt(x+k3) | |
| for c1 in range(state_vector[1].rows): | |
| runge_vectors[4].data[c1][0]=(runge_vectors[0].data[c1][0]+2.0*runge_vectors[1].data[c1][0]+2.0*runge_vectors[2].data[c1][0]+runge_vectors[3].data[c1][0])*(1/6.0) | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+runge_vectors[4].data[c1][0] | |
| return |
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 | |
| 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, dbydt_order): | |
| """ 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[dbydt_order-1].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[dbydt_order-1].data[c1][0]=0.0 # Added on Dec. 29, 2012. | |
| for c2 in range(matrix_b.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0] | |
| if (dbydt_order==2): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0] | |
| if (dbydt_order==3): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0] | |
| if (dbydt_order==4): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*ddt_mat[dbydt_order-2].data[c3][0] | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].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[dbydt_order-1].data[c1][0]-=matrix_e.data[c1][c3]*runge_vectors[4].data[c3][0]/dt | |
| ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0]*dt | |
| 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[dbydt_order-1].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 | |
| # Runge-Kutta 4th order method | |
| runge_function(state_vector[0], runge_vectors, 1) #k1=dt*(dx/dt) | |
| runge_function(state_vector[0], runge_vectors, 2) #k2=dt*d/dt(x+k1/2) | |
| runge_function(state_vector[0], runge_vectors, 3) #k3=dt*d/dt(x+k2/2) | |
| runge_function(state_vector[0], runge_vectors, 4) #k4=dt*d/dt(x+k3) | |
| for c1 in range(state_vector[1].rows): | |
| runge_vectors[4].data[c1][0]=(runge_vectors[0].data[c1][0]+2.0*runge_vectors[1].data[c1][0]+2.0*runge_vectors[2].data[c1][0]+runge_vectors[3].data[c1][0])*(1/6.0) | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+runge_vectors[4].data[c1][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, dt, 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() |
No comments:
Post a Comment