Re: Questions about BVP solver 1.1

[John Salvatier]

Hi Michael, 

I'm happy to hear you've found bvp_solver useful!

The version change does not indicate any major code changes, except for the fix of a long standing compilation bug. I don't actively work on the package much anymore, so I wanted to say it's officially 'mature'. 

I am unfortunately not the boundary value problem expert. I just took Boisvert, Muir and Spiteri's BVP_SOLVER Fortran package and wrapped it up in an easy to use Python package so it would be more widely accessible. Muir helped me out by switching to LAMPAK (for licensing issues; COLROW uses an overly restrictive license).

If you forward me the global error analysis thing, I'd be interested in taking a peek. I would also happily accept patches adding the functionality (the repository is here). 

Cheers,

John

P.S. I've forwarded this to the help mailing list (if that's a problem, I can delete the post no problem).

On Tue, Oct 23, 2012 at 11:51 PM, Michael Steinke <m.st...@lzh.de> wrote:

Hey John,

 

My name is Michael and I am working in Hannover, Germany in a laser development department, where I am (beside the work in the lab) trying to simulate the behavior of so called fibre amplifiers. One of my colleagues was using your bvp solver (version 0.3.0) to solve problems related to fibre lasers and it turned out, that the solver was doing some really nice work. Due to the fact that I am also dealing with boundary value problems, I also tried to use your solver for my simulations and up to some point I got also nice results. Nevertheless I am running into some problems, when the whole simulation becomes too complex (up to 500 coupled differential equations), resulting in failures during the solving. I am pretty sure, that I know why the solver is failing, so the reason, why I am writing this mail is not to get any support or hints, but for some questions about the newest version (1.1) of the solver:

 

I was reading the paper from Boisvert et al about the implementation of the global error analysis in the original fortran program, even though I am just a physicist and not a real mathematician and therefor I am pretty sure, that I just understand half of what was written in the paper. Nevertheless I thought that maybe the newer version of the solver could solve my problems and therefor I installed the new version of your solver, to see if it´s maybe working better (with respect to my problem) then the older version. Unfortunately I was quite surprised to see, that apparently there seems to be no implementation of this whole global error analysis in your solver?  Can I ask your, what real changes you made with respect to the older version? Is it just, that you switched from COLROW to LAMPAK or are there any other significant (with respect to stability of the solver) changes? Unfortunately I was not able to find any deeper information about the changes in the online-documentation of the solver.

 

I am looking forward to your answers and best regards

 

Michael

 

*********************************************************
Laser Zentrum Hannover e.V.
Dipl. Phys. Michael Steinke
Laser Development Department
Single-Frequency Lasers Group
Hollerithallee 8
D-30419 Hannover

Germany
Tel.: +49 511 2788-521
Fax: +49 511 2788-100
E-mail: m.st...@lzh.de
Internet: www.lzh.de

 

Geschäftsführender Vorstand/Executive Board Members:

Dr. rer. nat. Dietmar Kracht, Dipl.-Verw. (FH) Klaus Ulbrich

Vereinsregister-Nr.: 82 VR 5404, Amtsgericht Hannover

 

 

 

[Hwan]

John,

I just did the global error analysis. It looked sensible to me. scikits.bvp_solver does not save estimated yprime . But scipy.misc.derivative can be used to generate estimated yprime data. For instance,

from scipy.misc import derivative
def residuals(solution, t):
      yprime = derivative(solution, t, 0.001, 1)
      return dydt(t, solution(t)) - derivative(solution, t, 0.001, 1)

t = np.linspace(a, b, 100)
pylab.plot(t, residuals(solution, t)[0])   #for the first unknown
pylab.plot(t, residuals(solution, t)[1])   #for the second unknown
......
Note that "derivative()" uses central difference..
So, at each end of [a, b], one-sided difference should be applied. and this should be straightforward.

The default error tolerance is 1e-6. I wonder if scikits.bvp.solver allows to change the default tolerance. If not, the next version should make the error tolerance adjustable.

In general, scikits.bvp_solver is a very useful solver.

[John Salvatier]

Nice! Glad you got that worked out! 

You can change the tolerance with by using the tolerance argument: solve(...., tolerance = 1e-10)    

See the docstring for more: https://github.com/jsalvatier/scikits.bvp_solver/blob/master/scikits/bvp_solver/solver.py 

[Hwan]

John,

Great! Thanks for the information.

Hwan