[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