I apologize if this shows up twice. On 11/9/06, David L Goldsmith <David.L.Goldsmith@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