ODE solver's in scipy
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. Nils
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) -- 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
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
On Mon, Nov 04, 2002 at 01:51:19PM +0100, Nils Wagner wrote: [snip]
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 ?
First, you must return values from functions, not assign values to the name of the function. Incorrect: def A(): A = zeros(...) Correct: def A(): return zeros(...) or def A(): tmp = zeros(...) tmp[n,:] = ... return tmp Secondly, I don't think you are indexing correctly in A(). A[0:n,0:n] references the whole array, not just the diagonal. To reference a diagonal, A.flat[0:n:n+1] works if A is still contiguous. (Side-question to everybody else: Is there a convenience function somewhere in Numeric or scipy_base that allows setting the k-diagonal of a matrix? A.flat[0:n:n+1]?) Third, (this is what's causing the exception) you aren't actually calling P(t,z) and Pl(t,z) so Python is trying to index a function. I've attached my modifications. By the way, is z[n] supposed to change? Because the return value of F() has n elements, not n+1. I'm surprised the code runs, but it does. If it doesn't change, then you shouldn't have to recalculate P or Pl all the time. I have to hit the sack now, so I'll try to get back to this tomorrow.
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
Robert Kern schrieb:
On Mon, Nov 04, 2002 at 01:51:19PM +0100, Nils Wagner wrote:
[snip]
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 ?
First, you must return values from functions, not assign values to the name of the function.
Incorrect:
def A(): A = zeros(...)
Correct:
def A(): return zeros(...)
or
def A(): tmp = zeros(...) tmp[n,:] = ... return tmp
Secondly, I don't think you are indexing correctly in A(). A[0:n,0:n] references the whole array, not just the diagonal.
To reference a diagonal, A.flat[0:n:n+1] works if A is still contiguous.
(Side-question to everybody else: Is there a convenience function somewhere in Numeric or scipy_base that allows setting the k-diagonal of a matrix? A.flat[0:n:n+1]?)
Third, (this is what's causing the exception) you aren't actually calling P(t,z) and Pl(t,z) so Python is trying to index a function.
I've attached my modifications. By the way, is z[n] supposed to change? Because the return value of F() has n elements, not n+1. I'm surprised the code runs, but it does. If it doesn't change, then you shouldn't have to recalculate P or Pl all the time.
I have to hit the sack now, so I'll try to get back to this tomorrow.
Nils
--
Actually, A is a so-called bordered matrix. A = [A_11 a_12; a_21 0] where A_11 = (1-t)*(K-lambda*M)+ t*P a_12 = -(1-t)*M*x + t*Pl*x s_21 = x^T s_22 = 0 The first elements of z should contain the time-varying eigenvector x(t) The last element of z should be the eigenvalue, which varies with time, that is lambda=lambda(t). 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
------------------------------------------------------------------------
ivp_kern.pyName: ivp_kern.py Type: Plain Text (text/plain)
participants (2)
-
Nils Wagner -
Robert Kern