Anyone have a "little" shooting-method function to share

Pauli Virtanen pauli.virtanen at
Thu Nov 9 11:47:17 CST 2006

ke, 2006-11-08 kello 16:08 -0800, David L Goldsmith kirjoitti:
> Hi!  I tried to send this earlier: it made it into my sent mail folder, 
> but does not appear to have made it to the list.
> I need to numerically solve:
>     (1-t)x" + x' - x = f(t), x(0) = x0, x(1) = x1
> I've been trying to use (because it's the approach I inherited) an 
> elementary finite-difference discretization, but unit tests have shown 
> that that approach isn't working.  

You could try to use some of the pre-existing codes for
solving BVPs. This only requires writing a wrapper.

However, I have written a wrapper for the COLNEW code, which is a
finite-difference method written by Ascher & Bader in '87, having an
adaptive mesh selection. You can find it here:

As I understand, COLNEW does not have any special handling for singular
coefficients, but nevertheless it seems to be able to solve your
problem. The code goes as follows:

import scipy as N
import bvp

u0 = 1
u1 = 2

def f(t):
    return N.sin(2*t)

def fsub(t, z):
    u, du = z
    return N.array([(f(t) + u - du)/(1-t)])

def gsub(z):
    u, du = z
    return N.array([u[0] - u0, u[1] - u1])

tol = [1e-5, 1e-5]
boundary_points = [0, 1]

solution = bvp.colnew.solve(boundary_points, [2], fsub, gsub,
    is_linear=True, tolerances=tol,
    vectorized=True, maximum_mesh_size=300)

import pylab

x = solution.mesh
pylab.plot(x, solution(x)[:,0])


	Pauli Virtanen

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Digitaalisesti allekirjoitettu viestin osa
Url : 
-------------- next part --------------
Using Tomcat but need to do more? Need to support web services, security?
Get stuff done quickly with pre-integrated technology to make your job easier
Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo
-------------- next part --------------
Numpy-discussion mailing list
Numpy-discussion at

More information about the Numpy-discussion mailing list