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!










share|improve this question


























    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!










    share|improve this question
























      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!










      share|improve this question













      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






      share|improve this question













      share|improve this question











      share|improve this question




      share|improve this question










      asked Nov 18 at 14:19









      Mark

      366




      366
























          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).






          share|improve this answer





















            Your Answer






            StackExchange.ifUsing("editor", function () {
            StackExchange.using("externalEditor", function () {
            StackExchange.using("snippets", function () {
            StackExchange.snippets.init();
            });
            });
            }, "code-snippets");

            StackExchange.ready(function() {
            var channelOptions = {
            tags: "".split(" "),
            id: "1"
            };
            initTagRenderer("".split(" "), "".split(" "), channelOptions);

            StackExchange.using("externalEditor", function() {
            // Have to fire editor after snippets, if snippets enabled
            if (StackExchange.settings.snippets.snippetsEnabled) {
            StackExchange.using("snippets", function() {
            createEditor();
            });
            }
            else {
            createEditor();
            }
            });

            function createEditor() {
            StackExchange.prepareEditor({
            heartbeatType: 'answer',
            convertImagesToLinks: true,
            noModals: true,
            showLowRepImageUploadWarning: true,
            reputationToPostImages: 10,
            bindNavPrevention: true,
            postfix: "",
            imageUploader: {
            brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
            contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
            allowUrls: true
            },
            onDemand: true,
            discardSelector: ".discard-answer"
            ,immediatelyShowMarkdownHelp:true
            });


            }
            });














             

            draft saved


            draft discarded


















            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

























            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).






            share|improve this answer

























              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).






              share|improve this answer























                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).






                share|improve this answer












                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).







                share|improve this answer












                share|improve this answer



                share|improve this answer










                answered Nov 18 at 14:46









                9769953

                1,241311




                1,241311






























                     

                    draft saved


                    draft discarded



















































                     


                    draft saved


                    draft discarded














                    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





















































                    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







                    Popular posts from this blog

                    Costa Masnaga

                    Fotorealismo

                    Sidney Franklin