Solving dynamical ODE System with Python
up vote
0
down vote
favorite
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!
python ode solver
add a comment |
up vote
0
down vote
favorite
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!
python ode solver
add a comment |
up vote
0
down vote
favorite
up vote
0
down vote
favorite
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!
python ode solver
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!
python ode solver
python ode solver
asked Nov 18 at 14:19
Mark
366
366
add a comment |
add a comment |
1 Answer
1
active
oldest
votes
up vote
0
down vote
accepted
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).
add a comment |
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
up vote
0
down vote
accepted
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).
add a comment |
up vote
0
down vote
accepted
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).
add a comment |
up vote
0
down vote
accepted
up vote
0
down vote
accepted
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).
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).
answered Nov 18 at 14:46
9769953
1,241311
1,241311
add a comment |
add a comment |
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53361870%2fsolving-dynamical-ode-system-with-python%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown