#1
  1. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Oct 2013
    Posts
    3
    Rep Power
    0

    Python problem using ode45...


    hi to everybody...i have this problem.....

    this is my code:
    Code:
    import numpy as np
    from numpy import *
    from pylab import *
    from scipy.integrate import ode
    import matplotlib.pyplot as plt
    import warnings
    u=([[],[],[],[],[],[],[],[],[],[],[]])
    def molle(t,y):
        N = 1
        x = u0.flatten()[0:3:6]
        y = u0.flatten()[1:3:6]
        z = u0.flatten()[2:3:6]
        vx = u0.flatten()[7 + 0:3:12]
        vy = u0.flatten()[7 + 1:3:12]
        vz = u0.flatten()[7 + 2:3:12]
        for i in range(0,N):
              Fx = 0
              Fy = 0
              Fz = 0
              for j in range(0,N):
                    if A[i, j] != 0:
                                d = ([(x[i] - x[j]) ** 2 + (y[i] - y[j]) ** 2 + (z[i] - z[j]) ** 2]) ** 0.5
                                Fx = Fx + k[i, j] * (abs(x[j] - x[i]) - l[i, j]) * (x[j] - x[i]) * A[i, j] / d
                                Fy = Fy + k[i, j] * (abs(y[j] - y[i]) - l[i, j]) * (y[j] - y[i]) * A[i, j] / d
                                Fz = Fz + k[i, j] * (abs(z[j] - z[i]) - l[i, j]) * (z[j] - z[i]) * A[i, j] / d
                   
        u[0 + 2 * (i-1)] = u0[2 * N + 0 + 2 * (i-1)]    # dx
        u[1 + 2 * (i-1)] = u0[2 * N + 1 + 2 * (i-1)]    # dy
        u[2 + 2 * (i-1)] = u0[2 * N + 2 + 2 * (i-1)]    # dz
        u[2 * N + 0 + 2 * (i-1)] = Fx / m[i]    # ax=F/m
        u[2 * N + 1 + 2 * (i-1)] = Fy / m[i]    # ay=F/m
        u[2 * N + 2 + 2 * (i-1)] = Fz / m[i]   # az=F/m
        return u
    
    
    x0 = array([-1, 1])
    
    y0 = array([0, 0])
    
    z0 = array([0, 0])
    
    Vx0 = array([1, 0])
    Vy0 = array([0, 0])
    Vz0 = array([0, 0])
    
    A =array([[0,1],[1,0]])
    print A
    l = A
    k = A
    m = array([1, 1])
    
    
    u0=array([[x0[0]],[x0[1]],[y0[0]],[y0[1]],[z0[1]],[z0[1]],[Vx0[0]],[Vx0[1]],[Vy0[0]],[Vy0[1]],[Vz0[0]],[Vz0[1]]])
    
    
    # Create the time samples for the output of the ODE solver.
    # I use a large number of points, only because I want to make
    # a plot of the solution that looks nice.
    
    
    t0 = 0
    
    t1 = 10
    
    #backend = 'vode'
    backend = 'dopri5'
    #backend = 'dop853'
    
    solver = ode(molle).set_integrator(backend, nsteps=1)
    solver.set_initial_value(u0, t0)
    # suppress Fortran-printed warning
    solver._integrator.iwork[2] = -1
    
    sol = []
    warnings.filterwarnings("ignore", category=UserWarning)
    while solver.t < t1:
        solver.integrate(t1, step=True)
        sol.append([solver.t, solver.y])
    warnings.resetwarnings()
    sol = np.array(sol)
    
    plt.plot(sol[:,0], sol[:,1], 'b.-')
    plt.show()

    and this is the problem:

    Traceback (most recent call last):
    File "C:\Python27\Francesco\test.py", line 77, in <module>
    solver.integrate(t1, step=True)
    File "C:\Python27\lib\site-packages\scipy\integrate\_ode.py", line 386, in integrate
    self.f_params, self.jac_params)
    File "C:\Python27\lib\site-packages\scipy\integrate\_ode.py", line 888, in run
    tuple(self.call_args) + (f_params,)))
    ValueError: setting an array element with a sequence.

    can anyone modify correctly my code?
  2. #2
  3. Contributing User
    Devshed God 1st Plane (5500 - 5999 posts)

    Join Date
    Aug 2011
    Posts
    5,859
    Rep Power
    509
    It might help if you present the equations you intend to solve in either mathematica or LaTeX forms. Meanwhile, you need, it appears, to study python more deeply.

    import numpy as np
    from numpy import *

    Why do you import numpy as np and also into the entire module namespace? Because you could just write np.array near the end of your program as array . I'd import numpy and prefix my numpy routines with numpy. but that's not the common way to use numpy.

    You should know that (your N is 1) range(0,1) provides just 1 object. I suspect you want 2 objects.

    >>> print(list(range(0,1)))
    [0]

    And you should know that

    >>> print('abcdefghijklmnopqrstuvwxyz'[0:3:6])
    a

    slice(0,3,6) represents indexes on the interval [0, 3) in steps of 6.

    Maybe if you fix these sorts of problems the other error will vanish and we won't have to actually figure out what is the problem.
    [code]Code tags[/code] are essential for python code and Makefiles!
  4. #3
  5. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Oct 2013
    Posts
    3
    Rep Power
    0
    thanks...
    but i need to re-write in python my matlab code....
    here there are my m.file....if you watch it you could understand and help me...
    Matlab FILE:

    test.m
    Code:
    % esempio di 3 molle
    clear all
    clc
    x0=[-1 0 1];
    y0=[0 1/2 0];
    z0=[0 0 0];
    
    Vx0=[0 0 0];
    Vy0=[0 1 0];
    Vz0=[0 0 0];
    
    A=[0 1 0; 1 0 1; 0 1 0]
    
    l = A ;
    
    k = A ;
    
    m = [1 1 1];
    
    u0 = [x0';y0';z0';Vx0';Vy0';Vz0'];
    tspan=[0 10];
    [t,u] = ode45(@molle,tspan,u0,[],A,l,k,m);
    
    
    x=u(:,1:3:3*size(A,1));
    y=u(:,2:3:3*size(A,1));
    z=u(:,3:3:3*size(A,1));
    vx=u(:,3*size(A,1)+1:3:6*size(A,1));
    vy=u(:,3*size(A,1)+2:3:6*size(A,1));
    vz=u(:,3*size(A,1)+3:3:6*size(A,1));
    
    maxX=max(max(x));
    maxY=max(max(y));
    maxZ=max(max(z));
    minX=min(min(x));
    minY=min(min(y));
    minZ=min(min(z));
    
    
    figure,
    for k = 1 : length(t)
        [aa bb]=view;
        plot3(x(k,:),y(k,:),z(k,:),'*')
        view(aa,bb)
        axis([minX maxX minY maxY minZ maxZ])
        grid on
        pause(0.5)
    end

    molle.m (spring)
    Code:
    function u = molle(tspan,u0,A,l,k,m)
    
    N=size(A,1);
    
    x=u0(1:3:3*N)
    y=u0(2:3:3*N);
    z=u0(3:3:3*N);
    vx=u0(3*N+1:3:6*N);
    vy=u0(3*N+2:3:6*N);
    vz=u0(3*N+3:3:6*N);
    
    for i = 1 : N
        Fx=0;
        Fy=0;
        Fz=0;
        for j = 1 : N
            if A(i,j) ~= 0
            d=((x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2)^0.5;
            Fx=Fx+k(i,j)*(abs(x(j)-x(i))-l(i,j))*(x(j)-x(i))*A(i,j)/d;
            Fy=Fy+k(i,j)*(abs(y(j)-y(i))-l(i,j))*(y(j)-y(i))*A(i,j)/d; 
            Fz=Fz+k(i,j)*(abs(z(j)-z(i))-l(i,j))*(z(j)-z(i))*A(i,j)/d;
            end
        end
        u(1+3*(i-1))=u0(3*N+1+3*(i-1)); % dx
        u(2+3*(i-1))=u0(3*N+2+3*(i-1)); % dy
        u(3+3*(i-1))=u0(3*N+3+3*(i-1)); % dz
        u(3*N+1+3*(i-1))=Fx/m(i) ;  % ax=F/m
        u(3*N+2+3*(i-1))=Fy/m(i) ;  % ay=F/m
        u(3*N+3+3*(i-1))=Fz/m(i) ;  % az=F/m
    end
    
    u=u';

    where is my error??help me please....
  6. #4
  7. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Oct 2013
    Posts
    3
    Rep Power
    0
    here is my new code with some error correct:

    Code:
    import numpy as np
    from numpy import *
    from pylab import *
    from scipy.integrate import ode
    import matplotlib.pyplot as plt
    from scipy.special import gamma, airy
    import warnings
    
    
    def molle(u0,time):
        N = 1
        x = u0.flatten()[0:4:2]
        y = u0.flatten()[1:4:2]
        vx = u0.flatten()[4 + 0:9:2]
        vy = u0.flatten()[4 + 1:9:2]
        for i in range(0,N):
              Fx = 0
              Fy = 0
              for j in range(0,N):
                    if A[i, j] != 0:
                                d = ([(x[i] - x[j]) ** 2 + (y[i] - y[j]) ** 2 + (z[i] - z[j]) ** 2]) ** 0.5
                                Fx = Fx + k[i, j] * (abs(x[j] - x[i]) - l[i, j]) * (x[j] - x[i]) * A[i, j] / d
                                Fy = Fy + k[i, j] * (abs(y[j] - y[i]) - l[i, j]) * (y[j] - y[i]) * A[i, j] / d           
        u[0 + 2 * (i-1)] = u0[2 * N + 0 + 2 * (i-1)]    # dx
        u[1 + 2 * (i-1)] = u0[2 * N + 1 + 2 * (i-1)]    # dy
        u[2 * N + 0 + 2 * (i-1)] = Fx / m[i]    # ax=F/m
        u[2 * N + 1 + 2 * (i-1)] = Fy / m[i]    # ay=F/m
        return u
    
    
    x0 = array([-1, 1])
    print x0
    y0 = array([0, 0])
    print y0
    
    Vx0 = array([1, 0])
    Vy0 = array([0, 0])
    
    
    A =array([[0,1],[1,0]])
    print A
    l = A
    k = A
    m = array([1, 1])
    
    
    u0=array([[x0[0]],[x0[1]],[y0[0]],[y0[1]],[Vx0[0]],[Vx0[1]],[Vy0[0]],[Vy0[1]]])
    
    print u0
    # Create the time samples for the output of the ODE solver.
    # I use a large number of points, only because I want to make
    # a plot of the solution that looks nice.
    
    
    
    t0=0
    t1=10
    time = arange(0,10, 1)
    ychk = airy(time)[0]
    print ychk[:36:6]
    solver = ode(molle).set_integrator('vode', nsteps=1)
    print solver
    solver.set_initial_value(u0, t0)
    print solver
    sol=([[],[]])
    while solver.successful() and solver.t<t1:
        solver.integrate(t1)
    
    ###backend = 'vode'
    ##backend = 'dopri5'
    ###backend = 'dop853'
    ##
    ##solver = ode(molle).set_integrator(backend, nsteps=1)
    ##solver.set_initial_value(u0, t0)
    ###suppress Fortran-printed warning
    ##solver._integrator.iwork[2] = -1
    ##
    ##while solver.t<t1:
    ##    solver.integrate(t1, step=True)
    ##   
    ##sol = []
    ###warnings.filterwarnings("ignore", category=UserWarning)
    ##while solver.t < t1:
    ##    solver.integrate(t1, step=True)
    ##    sol.append([solver.t, solver.y])
    ##warnings.resetwarnings()
    ##sol = np.array(sol)
    ##
    ##plt.plot(sol[:,0], sol[:,1], 'b.-')
    ##plt.show()
    but i have this error:

    Traceback (most recent call last):
    File "C:\Python27\Francesco\test.py", line 67, in <module>
    solver.integrate(t1)
    File "C:\Python27\lib\site-packages\scipy\integrate\_ode.py", line 386, in integrate
    self.f_params, self.jac_params)
    File "C:\Python27\lib\site-packages\scipy\integrate\_ode.py", line 692, in run
    args[5:]))
    File "C:\Python27\Francesco\test.py", line 12, in molle
    x = u0.flatten()[0:4:2]
    AttributeError: 'float' object has no attribute 'flatten'
  8. #5
  9. Contributing User
    Devshed God 1st Plane (5500 - 5999 posts)

    Join Date
    Aug 2011
    Posts
    5,859
    Rep Power
    509
    What's z doing in there? This looks like a planar problem. Never mind, it was 3D in a prior post. And why do you copy Airy functions from the example into your code? You need to understand what problem you're solving.

    It's 3 springs with initial conditions.
    Last edited by b49P23TIvg; October 23rd, 2013 at 12:31 AM.
    [code]Code tags[/code] are essential for python code and Makefiles!
  10. #6
  11. Contributing User
    Devshed God 1st Plane (5500 - 5999 posts)

    Join Date
    Aug 2011
    Posts
    5,859
    Rep Power
    509
    Code:
                M0  S1  M1  S2   M2   S3   M3
    
                #-VVVVVV-#-VVVVVV-#-VVVVVVV-#
    The masses Mi have mass Mi
    The springs Sj have a natural length Lj and spring constant Kj and provide force Fj

    None of the springs are anchored.
    The positions Xi and velocity Vi of the masses are variable.
    Thus the number of degrees of freedom is twice the product of the spacial dimensionality by the number of masses.

    Mr. Newton said the sum of all forces at a mass equals the time rate of change of its momentum, which as it happens is relativistically correct, but we'll dumb it down a bit and just use mass times acceleration.

    The force generated by a spring is, to linear approximation which is pretty good for low stress on chunks of steel, the negative of amount it's stretched multiplied by the spring constant.

    The amount the spring is stretched is the natural length less the square root of the dot product of some vector with itself. That "some vector" is, and now I'll switch notation, for S[j] with j==i, X[i]-X[i-1].

    We're working with an integrator. The integrator provides the values of the variables at time t. Our function must compute the derivatives based on the variables. In this case we do not have a forcing function. None of the variables depend explicitly on time. After accounting for time in the argument list, we ignore time. The variables are a 1D vector. Their order doesn't matter to the ODE solver. However, we need to be consistent between them and the derivatives; they must match.

    The force on M2 is F2-F3. Force, position, velocity, acceleration are all vectors. The spring constants and natural lengths and masses are scalar. As is time, which we agreed to ignore.
    [code]Code tags[/code] are essential for python code and Makefiles!

IMN logo majestic logo threadwatch logo seochat tools logo