[SciPy-user] ODE solver's in scipy

[Nils Wagner]

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 at 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 at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user
-------------- next part --------------
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