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

Charles R Harris charlesr.harris at gmail.com
Sun Dec 17 12:31:33 CST 2006


I apologize if this shows up twice.

On 11/9/06, David L Goldsmith <David.L.Goldsmith at noaa.gov> wrote:
>
> Wow, thanks!
>
> DG
>
> Pauli Virtanen 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://projects.scipy.org/pipermail/numpy-discussion/attachments/20061217/6142c5bb/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: solution.zip
Type: application/zip
Size: 30179 bytes
Desc: not available
Url : http://projects.scipy.org/pipermail/numpy-discussion/attachments/20061217/6142c5bb/attachment-0001.zip 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: chebyshev.py
Type: text/x-python
Size: 18572 bytes
Desc: not available
Url : http://projects.scipy.org/pipermail/numpy-discussion/attachments/20061217/6142c5bb/attachment-0001.py 


More information about the Numpy-discussion mailing list