Robert Kern schrieb:
On Mon, Nov 04, 2002 at 12:34:36PM +0100, Nils Wagner wrote:
Hi,
How can I solve the following ODE in scipy
A(t,z) \dot{z} = f(z), z(0) = z_0
where A is a real matrix depending on t and z; \dot denotes derivative with respect to time t.
from scipy import linalg, integrate
def A(t,z): ...
def f(z): ...
def F(z, t): return linalg.solve(A(t,z), f(z))
t = arange(0, 100, 0.1) # or whatever z0 = array([...])
z = integrate.odeint(F, z0, t)
Thank you for your prompt reply. I have tried to transfer your suggestion to my problem. However, I get some errors. Traceback (most recent call last): File "ivp.py", line 46, in F return linalg.solve(A(P, Pl, M, K, t,z), f(P, M, K, z)) File "ivp.py", line 37, in A A[0:n,0:n] = (1.0-t)*(K[0:n,0:n]-z[n]*M[0:n,0:n])+t*P[0:n,0:n] TypeError: unsubscriptable object odepack.error: Error occured while calling the Python function named F So, what has to be amended in my program ivp.py ? Nils
-- Robert Kern Ruddock House President kern@caltech.edu
"In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter _______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
from scipy import * from scipy import linalg, integrate # # A homotopy method for a nonlinear eigenvalue problem # # h = (1-t) [K x-lambda M x] + t P(lambda) x = 0 # 0.5*(x^T x -1) = 0 # n = 3 # # Linear eigenvalue problem : K x = lambda M x # K = array(([1.0,-1.0,0.0], [-1.0,2.0,-1.0], [0.0,-1.0,2.0])) # M = array(([2.0,1.0,0.0], [1.0,4.0,1.0], [0.0,1.0,4.0]))/6.0 # w, vr = linalg.eig(K,M) def P(t,z): P = array(([z[n]/tan(z[n]),-z[n]/sin(z[n]),0.0], [-z[n]/sin(z[n]),2.0*z[n]/tan(z[n]),-z[n]/sin(z[n])], [0.0,-z[n]/sin(z[n]),2.0*z[n]/tan(z[n])])) def Pl(t,z): Pl = array(([cos(z[n])-z[n],z[n]*cos(z[n])-sin(z[n]),0.0], [z[n]*cos(z[n])-sin(z[n]),2.0*(cos(z[n])-z[n]),z[n]*cos(z[n])-sin(z[n])], [0.0,z[n]*cos(z[n])-sin(z[n]),2.0*(cos(z[n])-z[n])]))/sin(z[n])**2 def A(P, Pl, M, K, t,z): A = zeros((n+1,n+1),Float) A[0:n,0:n] = (1.0-t)*(K[0:n,0:n]-z[n]*M[0:n,0:n])+t*P[0:n,0:n] A[n:0,n] = -(1-t)*dot(M,z[0:n])+t*dot(Pl,z[0:n]) A[n:0,n] = z[0:n] def f(P, M, K, z): f = zeros(n+1,Float) f[0:n] = dot(K,z[0:n])-z[n]*dot(M,z[0:n])-dot(P,z[0:n]) def F(z, t): return linalg.solve(A(P, Pl, M, K, t,z), f(P, M, K, z)) t = arange(0, 0.1, 0.05) # or whatever z0=zeros(n+1,Float) z0[0:n] = vr[:,0] z0[n] = abs(w[0]) print z0 z = integrate.odeint(F, z0, t) print 'Eigenvector',vr[:,0] print 'Eigenvalue',w[0] print z