### Thread: Python problem using ode45...

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. 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.
3. 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....
4. 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'
5. 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.
6. 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.