Anyone have a "little" shooting-method function to share
![](https://secure.gravatar.com/avatar/dd7980fbcb8033634af99fca18f97062.jpg?s=120&d=mm&r=g)
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. After a little review, I believe I understand the problem: near t=0, the thing is like an advection-diffusion equation with diffusion as strong as advection (a case I'm sure is treated somewhere in the literature, but not as easily findable as:) near t=1, the thing is like the advection-diffusion equation I did easily find treated, namely one where advection dominates. Based on that treatment, I understand why an FD approach (with a uniform grid) would fail around t=0, and from there it is easy to accept that the fact that the thing changes its "nature" over the course of its "life," implies that (a uniform grid) FD is probably not a very good approach for this equation anywhere on its domain. I could try (and maybe will have to) a variable grid FD approach, but I'd first like to try a "shooting" method (for some reason my intuition is telling me that in this case this might be more efficient). This is where you all come in: I understand the algorithm and could (may have to) code it myself, but if anyone out there already has code for this (in Python using numpy, and preferably already under test), might you be willing to share? Thanks in advance, David Goldsmith PS: In the interim, I realized that if something like this "pre-exists" in a module, that module might be scipy. Sure enough, scipy has integrate.ode and integrate.odeint; Google-ing scipy: integrate.odeint help led me to Travis's 2004 "SciPy Tutorial," where I find "There are many optional inputs and outputs available when using odeint which can help tune the solver. These additional inputs and outputs are not needed much of the time, however..." Well, is one of those inputs an algorithm specifier? Alternatively, is there a "shooting" algorithm implemented elsewhere in scipy? ------------------------------------------------------------------------- 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 http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642
![](https://secure.gravatar.com/avatar/2fec026482d0ffd708974ff4c5584ebd.jpg?s=120&d=mm&r=g)
ke, 2006-11-08 kello 16:08 -0800, David L Goldsmith kirjoitti:
You could try to use some of the pre-existing www.netlib.org 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: http://www.iki.fi/pav/bvp 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]) pylab.savefig('solution.png') ---------------------------------------------------- BR, Pauli Virtanen ------------------------------------------------------------------------- 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 http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
![](https://secure.gravatar.com/avatar/dd7980fbcb8033634af99fca18f97062.jpg?s=120&d=mm&r=g)
Wow, thanks! DG Pauli Virtanen wrote:
-- HMRD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/emergencyresponse/> ------------------------------------------------------------------------- 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 http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
I apologize if this shows up twice. On 11/9/06, David L Goldsmith <David.L.Goldsmith@noaa.gov> wrote:
<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
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 11/9/06, Pauli Virtanen <pauli.virtanen@iki.fi> wrote:
<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
![](https://secure.gravatar.com/avatar/2fec026482d0ffd708974ff4c5584ebd.jpg?s=120&d=mm&r=g)
ke, 2006-11-08 kello 16:08 -0800, David L Goldsmith kirjoitti:
You could try to use some of the pre-existing www.netlib.org 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: http://www.iki.fi/pav/bvp 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]) pylab.savefig('solution.png') ---------------------------------------------------- BR, Pauli Virtanen ------------------------------------------------------------------------- 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 http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
![](https://secure.gravatar.com/avatar/dd7980fbcb8033634af99fca18f97062.jpg?s=120&d=mm&r=g)
Wow, thanks! DG Pauli Virtanen wrote:
-- HMRD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/emergencyresponse/> ------------------------------------------------------------------------- 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 http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
I apologize if this shows up twice. On 11/9/06, David L Goldsmith <David.L.Goldsmith@noaa.gov> wrote:
<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
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On 11/9/06, Pauli Virtanen <pauli.virtanen@iki.fi> wrote:
<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
participants (3)
-
Charles R Harris
-
David L Goldsmith
-
Pauli Virtanen