Nsryjdtyk

I am trying to solve a dynamical system with three state variables V1,V2,I3 and then plot these in a 3d Plot. My code so far looks as follows:

from scipy.integrate import ode

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import math

def ID(V,a,b):

    return a*(math.exp(b*V)-math.exp(-b*V))





def dynamical_system(t,z,C1,C2,L,R1,R2,R3,RN,a,b):



    V1,V2,I3 = z

    f = [(1/C1)*(V1*(1/RN-1/R1)-ID(V1-V2,a,b)-(V1-V2)/R2),(1/C2)*(ID(V1-V2,a,b)+(V1-V2)/R2-I3),(1/L)*(-I3*R3+V2)]

    return f 



# Create an `ode` instance to solve the system of differential

# equations defined by `dynamical_system`, and set the solver method to 'dopri5'.

solver = ode(dynamical_system)

solver.set_integrator('dopri5')



# Set the initial value z(0) = z0.

C1=10

C2=100

L=0.32

R1=22

R2=14.5

R3=100

RN=6.9

a=2.295*10**(-5)

b=3.0038

solver.set_f_params(C1,C2,L,R1,R2,R3,RN,a,b)

t0 = 0.0

z0 = [-3, 0.5, 0.25] #here you can set the inital values V1,V2,I3

solver.set_initial_value(z0, t0)





# Create the array `t` of time values at which to compute

# the solution, and create an array to hold the solution.

# Put the initial value in the solution array.

t1 = 25

N = 200 #number of iterations

t = np.linspace(t0, t1, N)

sol = np.empty((N, 3))

sol[0] = z0



# Repeatedly call the `integrate` method to advance the

# solution to time t[k], and save the solution in sol[k].

k = 1

while solver.successful() and solver.t < t1:

    solver.integrate(t[k])

    sol[k] = solver.y

    k += 1





xlim = (-4,1)

ylim= (-1,1)

zlim=(-1,1)

fig=plt.figure()

ax=fig.gca(projection='3d')

#ax.view_init(35,-28)

ax.set_xlim(xlim)

ax.set_ylim(ylim)

ax.set_zlim(zlim)



print sol[:,0]

print sol[:,1]

print sol[:,2]

ax.plot3D(sol[:,0], sol[:,1], sol[:,2], 'gray')



plt.show()

Printing the arrays that should hold the solutions sol[:,0] etc. shows that apparently it constantly fills it with the initial value. Can anyone help? Thanks!

Use from __future__ import division.

I can't reproduce your problem: I see a gradual change from -3 to -2.46838127, from 0.5 to 0.38022886 and from 0.25 to 0.00380239 (with a sharp change from 0.25 to 0.00498674 in the first step). This is with Python 3.7.0, NumPy version 1.15.3 and SciPy version 1.1.0.

Given that you are using Python 2.7, integer division may be the culprit here. Quite a number of your constants are integer, and you have a bunch of 1/<constant> integer divisions in your equation.

Indeed, if I replace / with // in my version (for Python 3), I can reproduce your problem.

Simply add from __future__ import division at the top of your script to solve your problem.

Also add from __future__ import print_function at the top, replace print <something> with print(<something>) and your script is fully Python 3 and 2 compatible).