[SciPy-user] ODE solver's in scipy

[Robert Kern]

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 at caltech.edu

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter
-------------- next part --------------
from scipy import *
from scipy import linalg, integrate
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):

    return 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):

    return 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):
    tmp = zeros((n,n),Float)
    tmp.flat[0:n:n+1] = (1.0-t)*(K.flat[0:n:n+1]-z[n]*M.flat[0:n:n+1])+t*P.flat[0:n:n+1]
    tmp[0:n,n-1]   = -(1-t)*dot(M,z[:n])+t*dot(Pl,z[:n])
    tmp[n-1,0:n]   = z[:n]
    return tmp

def f(P, M, K, t, z):
    tmp = dot(K,z[:n])-z[n]*dot(M,z[:n])-dot(P,z[:n])
    tmp[n-1] = 0.0
    return tmp

def F(z, t):
    p = P(t,z)
    pl = Pl(t,z)
    x = linalg.solve(A(p, pl, M, K, t, z), f(p, M, K, t, z))
    return x

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 z