[Numpy-discussion] Anyone have a "little" shooting-method function to share

[Charles R Harris]

On 11/9/06, Pauli Virtanen <pauli.virtanen at iki.fi> wrote:
>
> 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.

<snip>

I thought I would try a Chebyshev method on this problem. My solution with
order ten (degree 9) Chebyshev polynomials goes as follows.

In [111]: import chebyshev as c

In [112]: t = c.modified_points(10,0,1)   # use 10 sample points

In [113]: D = c.modified_derivative(10,0,1) # derivative operator

In [114]: op = (1.0 - t)[:,newaxis]*dot(D,D) + D - eye(10) # differential
equation

In [115]: op[0] = 0 # set up boundary condition y(0) = y0

In [116]: op[0,0] = 1

In [117]: op[9] = 0 # set up boundary condition y(1) = y1

In [118]: op[9,9] = 1

In [119]: opinv = alg.inv(op) # invert the operator

In [120]: f = exp(t) # try  f(t) = exp(t)

In [121]: f[0] = 2 # y0 = 2

In [122]: f[9] = 1 # y1 = 1

In [123]: soln = dot(opinv,f) # solve equation

In [124]: plot(t,soln)
Out[124]: [<matplotlib.lines.Line2D instance at 0x976eaec>]

The plot is rather rough with only 10 points. Replot with more.

In [125]: tsmp = linspace(0,1)

In [126]: interp = c.modified_values(tsmp, 10, 0, 0, 1)

In [127]: plot(tsmp, dot(interp, soln))
Out[127]: [<matplotlib.lines.Line2D instance at 0x977ca4c>]


Looks OK here. You can save opinv as it doesn't change with f. Likewise, if
you always want to interpolate the result, then save dot(interp, opinv).
I've attached a plot of the solution I got along with the chebyshev module I
use.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20061217/f3e26293/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: solution.jpg.zip
Type: application/zip
Size: 16196 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20061217/f3e26293/attachment.zip>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: chebyshev.py.zip
Type: application/zip
Size: 3997 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20061217/f3e26293/attachment-0001.zip>