import numpy as np
from scipy import integrate
#import ompc
from matplotlib.pylab import *
from numpy import *
from numpy import matrix
from numpy import linalg
import scipy
from matplotlib import pyplot
def molle(tspan, u0, A, l, k, m):
N = 3
x = u0[0:N*3:3]
y = u0[1:N*3:3]
z = u0[2:N*3:3]
vx = u0[3*N+0:6*N:3]
vy = u0[3*N+1:6*N:3]
vz = u0[3*N+2:6*N:3]
# Reshape u0 for ease
u0 = u0.reshape( [u0.shape[0]*u0.shape[1],1],order='F' )
# Create output array
u = np.zeros(u0.shape)
for i in xrange(0,N):
Fx = 0
Fy = 0
Fz = 0
for j in xrange(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
#-------------------------------------------------------------------------------------------
# No more meaningful changes to this function...
u(1 + 3 * (i - 1)).lvalue = u0[3 * N + 1 + 3 * (i - 1)-1] # dx
u(2 + 3 * (i - 1)).lvalue = u0[3 * N + 2 + 3 * (i - 1)-1] # dy
u(3 + 3 * (i - 1)).lvalue = u0[3 * N + 3 + 3 * (i - 1)-1] # dz
u(3 * N + 1 + 3 * (i - 1)).lvalue = Fx / m[i] # ax=F/m
u(3 * N + 2 + 3 * (i - 1)).lvalue = Fy / m[i] # ay=F/m
u(3 * N + 3 + 3 * (i - 1)).lvalue = Fz / m[i] # az=F/m
u = u.cT
x0=matrix([-1,0,1])
y0=matrix([0,0.5,0])
z0=matrix([0,0,0])
Vx0=matrix([0,0,0])
Vy0=matrix([0,1,0])
Vz0=matrix([0,0,0])
A=matrix([[0,1,0],[1,0,1],[0,1,0]])
print A
l=A
k=A
m=([1,1,1])
u0=array([transpose(x0),transpose(y0),transpose(z0),transpose(Vx0),transpose(Vy0),transpose(Vz0)])
tspan=matrix([0,10])
[t,u]=scipy.integrate.ode(molle(tspan, u0, A, l, k, m)).set_integrator('dopri5')
import numpy as np
from scipy.integrate import ode
import scipy
from pylab import *
def molle(tspan,u0,A,k,l,m):
#
u=([[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]])
N = 2
x = u0.flatten()[0:3:8]
y = u0.flatten()[1:3:8]
z = u0.flatten()[2:3:8]
vx = u0.flatten()[8 + 0:3:17]
vy = u0.flatten()[8 + 1:3:17]
vz = u0.flatten()[8 + 2:3:17]
for i in range(0,N-1):
Fx = 0
Fy = 0
Fz = 0
print "ciccio"
for j in range(0,N-1):
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
print "ciao"
return transpose(u)
x0 = array([-1, 0, 1])
y0 = array([0, 0.5, 0])
z0 = array([0, 0, 0])
Vx0 = array([0, 0, 0])
Vy0 = array([0, 1, 0])
Vz0 = array([0, 0, 0])
A =array([[0,1,0],[1,0,1],[0,1,0]])
l = A
k = A
m = array([1, 1, 1])
u0=array([[transpose(x0)],[transpose(y0)],[transpose(z0)],[transpose(Vx0)],[transpose(Vy0)],[transpose(Vz0)]])
tspan=linspace(0,10,100)
u = scipy.integrate.ode(molle(tspan,u0,A,l,k,m)).set_integrator('dopri5')
anastasi.fr wrote:it'all python code but it not works
....now can you help me?
anastasi.fr wrote:i am writing in matlab a code to simulate a mass spring system
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, time)
print solver
sol=([[],[]])
while solver.successful():
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()
Users browsing this forum: magic_coding and 3 guests