This was an error I am surprised I was even able to locate to the solver because I thought it was something related to the circuit solver. So far, the ODE is solved as an entire state vector at a time. This means, all states are calculated together and when one state needs the information of another state it used the previous value of the state. So even if a more updated value of a state is available, with this method the previous value is used.
This actually defeats the purpose of solving ODEs in a backward manner. That is we solve from the last state to the first state since the matrices have been made upper triangular.
So, the ODE was changed such that a particular state is calculated completely, its value is updated and this updated value is available for the next state to be calculated.
Here is the code (click on "view raw" below the code box to view the code in another window):
This actually defeats the purpose of solving ODEs in a backward manner. That is we solve from the last state to the first state since the matrices have been made upper triangular.
So, the ODE was changed such that a particular state is calculated completely, its value is updated and this updated value is available for the next state to be calculated.
Here is the code (click on "view raw" below the code box to view the code in another window):
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, ode_vectors, sys_loops, list_of_nodes): | |
| """ Solves the ODE using Runge Kutta 4th order method. """ | |
| def runge_function4(x_plus_dxdt, ddt_mat, dbydt_order, ode_row): | |
| """ Defines the function dx/dt=f(x) for 4th order | |
| Runge Kutta solver. """ | |
| #Calculate d/dt vector in reverse. | |
| #for c1 in range(matrix_e.rows-1, -1, -1): | |
| if matrix_e.data[ode_row][ode_row]: | |
| 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[ode_row][0]=matrix_b.data[ode_row][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[ode_row][0]=0.0 # Added on Dec. 29, 2012. | |
| for c2 in range(matrix_b.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0] | |
| if (dbydt_order==2): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][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[ode_row][0]-=matrix_a.data[ode_row][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[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0] | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*x_plus_dxdt.data[c3][0] | |
| ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]*dt | |
| for c3 in range(ode_row+1, matrix_e.columns): | |
| # for c3 in range(matrix_e.columns): | |
| if not ode_row==c3: | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_e.data[ode_row][c3]*ode_vectors[4].data[c3][0] | |
| ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0] | |
| # ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]/matrix_e.data[ode_row][ode_row] | |
| else: | |
| if not (matrix_a.data[ode_row][ode_row]): | |
| #The variable has no dynamics and can't even be calculated statically! | |
| # May be due to a redundant loop. | |
| ddt_mat[dbydt_order-1].data[ode_row][0]=0.0 | |
| state_vector[0].data[ode_row][0]=0.0 | |
| state_vector[1].data[ode_row][0]=0.0 | |
| else: | |
| ddt_mat[dbydt_order-1].data[ode_row][0]=0.0 | |
| state_vector[0].data[ode_row][0]=0.0 | |
| state_vector[1].data[ode_row][0]=0.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[ode_row][0]+=matrix_b.data[ode_row][0]*input_vector | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][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[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0] | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0] | |
| for c2 in range(matrix_a.columns): | |
| if not (ode_row==c2): | |
| state_vector[0].data[ode_row][0]-=matrix_a.data[ode_row][c2]*state_vector[1].data[c2][0] | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0] | |
| state_vector[0].data[ode_row][0]=state_vector[0].data[ode_row][0]/matrix_a.data[ode_row][ode_row] | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0] | |
| return | |
| def runge_function5(x_plus_dxdt, ddt_mat, dbydt_order, ode_row): | |
| """ Defines the function dx/dt=f(x) for higher order | |
| ODE solver """ | |
| #Calculate d/dt vector in reverse. | |
| #for c1 in range(matrix_e.rows-1, -1, -1): | |
| if matrix_e.data[ode_row][ode_row]: | |
| 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[ode_row][0]=matrix_b.data[ode_row][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[ode_row][0]=0.0 # Added on Dec. 29, 2012. | |
| for c2 in range(matrix_b.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0] | |
| if (dbydt_order==2): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*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[ode_row][0]-=matrix_a.data[ode_row][c3]*0.5*ddt_mat[dbydt_order-3].data[c3][0] | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][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[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-4].data[c3][0]*(14.0/64.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-3].data[c3][0]*(5.0/64.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]*(-3.0/64.0) | |
| if (dbydt_order==5): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-5].data[c3][0]*(-12.0/96.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-4].data[c3][0]*(-12.0/96.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-3].data[c3][0]*(8.0/96.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]*(64.0/96.0) | |
| if (dbydt_order==6): | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-5].data[c3][0]*(-9.0/64.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-4].data[c3][0]*(5.0/64.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-3].data[c3][0]*(16.0/64.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]*(36.0/64.0) | |
| for c3 in range(matrix_a.columns): | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*x_plus_dxdt.data[c3][0] | |
| ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]*dt | |
| for c3 in range(matrix_e.columns): | |
| if not ode_row==c3: | |
| ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_e.data[ode_row][c3]*ode_vectors[4].data[c3][0] | |
| ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0] | |
| else: | |
| if not (matrix_a.data[ode_row][ode_row]): | |
| #The variable has no dynamics and can't even be calculated statically! | |
| # May be due to a redundant loop. | |
| state_vector[0].data[ode_row][0]=0.0 | |
| state_vector[1].data[ode_row][0]=0.0 | |
| else: | |
| ddt_mat[dbydt_order-1].data[ode_row][0]=0.0 | |
| state_vector[0].data[ode_row][0]=0.0 | |
| state_vector[1].data[ode_row][0]=0.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[ode_row][0]+=matrix_b.data[ode_row][0]*input_vector | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][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[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0] | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0] | |
| for c2 in range(matrix_a.columns): | |
| if not (ode_row==c2): | |
| state_vector[0].data[ode_row][0]-=matrix_a.data[ode_row][c2]*state_vector[1].data[c2][0] | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0] | |
| state_vector[0].data[ode_row][0]=state_vector[0].data[ode_row][0]/matrix_a.data[ode_row][ode_row] | |
| state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0] | |
| return | |
| def trapezoidal_mthd(x_plus_dxdt, ddt_mat, dbydt_order): | |
| """ Defines the function dx/dt=f(x) for Trapezoidal | |
| method. """ | |
| #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] | |
| 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] | |
| ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0]*dt | |
| for c3 in range(matrix_e.columns): | |
| if not c1==c3: | |
| ddt_mat[dbydt_order-1].data[c1][0]-=matrix_e.data[c1][c3]*ode_vectors[2].data[c3][0] | |
| ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0] | |
| else: | |
| if not (matrix_a.data[c1][c1]): | |
| #The variable has no dynamics and can't even be calculated statically! | |
| # May be due to a redundant loop. | |
| 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 | |
| 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 | |
| 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] | |
| 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] | |
| 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 | |
| # Trapezoidal method | |
| ## trapezoidal_mthd(state_vector[0], ode_vectors, 2) #k1=dt*(dx/dt) | |
| ## ode_vectors[2].data[c1][0]=(ode_vectors[0].data[c1][0]+ode_vectors[1].data[c1][0])/2.0 | |
| ## state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[2].data[c1][0] | |
| ## ode_vectors[0].data[c1][0]=ode_vectors[1].data[c1][0] | |
| # Runge-Kutta 4th order method | |
| for c1 in range(matrix_e.rows-1, -1, -1): | |
| runge_function4(state_vector[0], ode_vectors, 1, c1) #k1=dt*(dx/dt) | |
| runge_function4(state_vector[0], ode_vectors, 2, c1) #k2=dt*d/dt(x+k1/2) | |
| runge_function4(state_vector[0], ode_vectors, 3, c1) #k3=dt*d/dt(x+k2/2) | |
| runge_function4(state_vector[0], ode_vectors, 4, c1) #k4=dt*d/dt(x+k3) | |
| ode_vectors[4].data[c1][0]=(ode_vectors[0].data[c1][0]+2.0*ode_vectors[1].data[c1][0]+2.0*ode_vectors[2].data[c1][0]+ode_vectors[3].data[c1][0])/6.0 | |
| state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[4].data[c1][0] | |
| # Runge-Kutta 5th order method | |
| # for c1 in range(state_vector[1].rows): | |
| # runge_function5(state_vector[0], ode_vectors, 1, c1) #k1=dt*(dx/dt) | |
| # runge_function5(state_vector[0], ode_vectors, 2, c1) #k2=dt*d/dt(x+k1) | |
| # runge_function5(state_vector[0], ode_vectors, 3, c1) #k3=dt*d/dt(x+(k1+k2)/2) | |
| # runge_function5(state_vector[0], ode_vectors, 4, c1) #k4=dt*d/dt(x+(14k1+5k2-3k3)/64) | |
| # runge_function5(state_vector[0], ode_vectors, 5, c1) #k4=dt*d/dt(x+(-12k1-12k2+8k3+64k4)/96) | |
| # runge_function5(state_vector[0], ode_vectors, 6, c1) #k4=dt*d/dt(x+(-9k2+5k3+16k4+36k5)/64) | |
| # ode_vectors[4].data[c1][0]=(ode_vectors[0].data[c1][0]+2.0*ode_vectors[1].data[c1][0]+2.0*ode_vectors[2].data[c1][0]+ode_vectors[3].data[c1][0])/6.0 | |
| # state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[4].data[c1][0] | |
| ## ode_vectors[6].data[c1][0]=7.0*ode_vectors[0].data[c1][0] | |
| ## ode_vectors[6].data[c1][0]+=12.0*ode_vectors[2].data[c1][0] | |
| ## ode_vectors[6].data[c1][0]+=7.0*ode_vectors[3].data[c1][0] | |
| ## ode_vectors[6].data[c1][0]+=32.0*ode_vectors[4].data[c1][0] | |
| ## ode_vectors[6].data[c1][0]+=32.0*ode_vectors[5].data[c1][0] | |
| ## ode_vectors[6].data[c1][0]=ode_vectors[6].data[c1][0]/90.0 | |
| ## state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[6].data[c1][0] | |
| return |
No comments:
Post a Comment