|
|
12
|
|
Hi All,
As someone with a keen interest in mathematics and software development I am interested in getting involved with Octave's development. As a starting place I have been looking into the odepkg. It seems pretty comprehensive apart from the bvp4c function ( http://octave.sourceforge.net/odepkg/function/bvp4c.html). The examples on this webpage fail to run and, from previous discussions on this mailing list, I am under the impression this function needs some work. This seems like a good starting place for me.
The current specification for bvp4c is:
bvp4c (odefun, bcfun, solinit).
In MATLAB solinit can be provided by the function:
solinit = bvpinit(x, yinit, params)
However, in Octave, I am not aware of an analogous bvpinit. Is this something that should be implemented before further work is done on bvp4c?
Best Regards,
Luke
|
|
|
Den 08-06-2016 kl. 16:22 skrev lakerluke:
> As someone with a keen interest in mathematics and software development I am
> interested in getting involved with Octave's development. As a starting
> place I have been looking into the odepkg. It seems pretty comprehensive
> apart from the bvp4c function
> ( http://octave.sourceforge.net/odepkg/function/bvp4c.html). The examples on
> this webpage fail to run and, from previous discussions on this mailing
> list, I am under the impression this function needs some work. This seems
> like a good starting place for me.
>
> The current specification for bvp4c is:
>
> bvp4c (odefun, bcfun, solinit).
>
> In MATLAB solinit can be provided by the function:
>
> solinit = bvpinit(x, yinit, params)
>
> However, in Octave, I am not aware of an analogous bvpinit. Is this
> something that should be implemented before further work is done on bvp4c?
I think having a working implementation of bvpinit would be very helpful
when moving code from Matlab to Octave.
Søren
|
|
|
Il 08/Lug/2016 12:53, "Carlo de Falco" <carlo.defalco@polimi.it> ha scritto:
>
>>
>> I would like to work on improving this and so I have started by trying to understand the code that is already written. There are certain hard coded values in this file that give me the impression that this is just an implementation of Simpson's rule. Am I right in thinking this? For example:
>>
>> line 69:
>>
>> s = 3;
>>
>> line 162:
>>
>> persistent B= [1/6 2/3 1/6 ];
>>
>> I am aware that matlab's implementation of bvp4c uses a three stage Lobatto IIIA formula
The two are actually equivalent in this context according to the authors of bvp4c :
http://www.orcca.on.ca/TechReports/TechReports/2001/TR-01-02.pdf>>and I was wondering what are the main things to look at in Octave's current bvp4c file in order to progress it.
the current implementation is very crude , there's lots to be done ...
>>
>> I am currently using this as a resource for how to improve this code:
>>
>> http://homepage.math.uiowa.edu/~ljay/publications.dir/Lobatto.pdf>>
>> If you know of any other useful resources for this project please let me know.
>>
>> Best regards,
>>
>> Luke
>>
>>
Thank you for the response and the resource, Carlo. I am now looking through that pdf and analysing the code in bvp4c.m to look at ways of improving the script. If anyone else has any resources they think would be useful for this project, please let me know.
Luke
|
|
|
On 12 Jul 2016, at 12:59, lakerluke < [hidden email]> wrote:
> Thank you for the response and the resource, Carlo. I am now looking through
> that pdf and analysing the code in bvp4c.m to look at ways of improving the
> script. If anyone else has any resources they think would be useful for this
> project, please let me know.
Hi Luke,
Thanks again for your interest in contributing to
bvp4c for Octave.
Others on this list have expressed their
interest in improving this function,
so if anyone has done anything in that
direction I hope he/she jumps into
this discussion to discuss the current status
and avoid duplication of effort.
As for my initial implementation, as I already
said it is extremely crude and I don't think
it reflects the correct algorithm described in
the Shampine paper [1].
All it does at the moment is set up the residual
for a nonlinear problem using Simpson's rule and
let fsolve do the rest.
The error estimate for mesh refinement is even more
naive as it just compares y' computed with first order
finite differences to the evaluation of f(t,y) at interval
midpoints instead of applying the integral-norm-based
approach of [1].
HTH,
c.
[1] Kierzenka, Jacek, and Lawrence F. Shampine.
"A BVP solver based on residual control and the Maltab PSE."
ACM Transactions on Mathematical Software (TOMS) 27.3 (2001): 299-316.
|
|
|
Hi Carlo,
I have been doing some investigation into bvp4c and would like to ask you some questions relating to the Shampine paper you linked. I believe that the first step in implementing bvp4c, in the manner of the algorithm outlined in the Shampine paper, is to form the residual approximation that appears at the middle of page 5 (I will refer to this equation as [1]). Then, the integral norm must be taken of this residual in order to base error estimation and mesh selection.
The equation for the residual [1] depends on the 4th derivative of the solution y(x). My question is how best to form this term? Since dy/dx = f(t, y) we can use the chain rule and take continuous derivatives of this governing equation to form the 4th derivative of y in terms of f and derivatives of f. However, this will need recalculating at each iteration and I don't believe this is the best method to use.
The Shampine paper goes on to explain on page 6 that we can apply the Simpson method to system (4.1) and further explains some simplification that can be done for the Jacobian terms. However, applying the Simpson method results in a cubic spline and hence if this is used as the approximation to y(x) then upon taking the 4th derivative of the spline we will get a constant function of zero at all points?
Finally, at the top of page 9, the paper expresses the residual in terms of the solution to the system (4.1), (I will refer to this expression of the residual at the midpoint as [2]). Since the residual at the midpoint is an indication of how well our equations were satisfied, then this expression [2] is the residual we wish to minimize. It then goes on to say that [2] is used by the quadrature formula for the norm of the residual (i.e. [1] mentioned earlier). I fail to see how [2] can be applied to [1] since they are evaluated at different points.
Apologies for the length email but I'm hoping someone could help me out with my questions:
1) How to best form the term for the 4th derivative of the solution evaluated at the midpoint in order to form [1]?
2) How expression [2] is used in the evaluation of [1]?
Best Regards,
Luke
|
|
|
On 3 Aug 2016, at 08:53, lakerluke < [hidden email]> wrote:
> Hi Carlo,
>
> I have been doing some investigation into bvp4c and would like to ask you
> some questions relating to the Shampine paper you linked. I believe that the
> first step in implementing bvp4c, in the manner of the algorithm outlined in
> the Shampine paper, is to form the residual approximation that appears at
> the middle of page 5 (I will refer to this equation as [1]). Then, the
> integral norm must be taken of this residual in order to base error
> estimation and mesh selection.
>
> The equation for the residual [1] depends on the 4th derivative of the
> solution y(x). My question is how best to form this term? Since dy/dx = f(t,
> y) we can use the chain rule and take continuous derivatives of this
> governing equation to form the 4th derivative of y in terms of f and
> derivatives of f. However, this will need recalculating at each iteration
> and I don't believe this is the best method to use.
>
> The Shampine paper goes on to explain on page 6 that we can apply the
> Simpson method to system (4.1) and further explains some simplification that
> can be done for the Jacobian terms. However, applying the Simpson method
> results in a cubic spline and hence if this is used as the approximation to
> y(x) then upon taking the 4th derivative of the spline we will get a
> constant function of zero at all points?
>
> Finally, at the top of page 9, the paper expresses the residual in terms of
> the solution to the system (4.1), (I will refer to this expression of the
> residual at the midpoint as [2]). Since the residual at the midpoint is an
> indication of how well our equations were satisfied, then this expression
> [2] is the residual we wish to minimize. It then goes on to say that [2] is
> used by the quadrature formula for the norm of the residual (i.e. [1]
> mentioned earlier). I fail to see how [2] can be applied to [1] since they
> are evaluated at different points.
>
> Apologies for the length email but I'm hoping someone could help me out with
> my questions:
>
> 1) How to best form the term for the 4th derivative of the solution
> evaluated at the midpoint in order to form [1]?
>
> 2) How expression [2] is used in the evaluation of [1]?
>
> Best Regards,
>
> Luke
Hi Luke,
Thanks for working on this!
I am traveling and I am unfortunately not able
to give a detailed answer right now, will try to
do so when I am back at the end of next week.
Anyway I wanted to give a quick hint as I am afraid
you might spend too much time workin in the wrong direction...
One thing you have to keep in mind is that
you will never need to compute any high order
derivatives, this is an important feature of the
method.
I can only guess as I don't have the paper at hand
right now, but I am quite sure that, if the formula
you refer to includes a fourth-order derivative, than
that formula is meant for the (theoretical) analysis
of the error, it is not an expression to be used in
the implementation.
I'll come back to this as soon as I'm back,
c.
|
|
|
Thanks for the pointer! That really helps. So don't worry about answering my questions from the previous post. In light of your comment I think I understand the approach intended by the Shampine paper. We use Newton's method to look for roots for the system (4.1) which gives us our solution at mesh points. As pointed out in the paper in order to apply this method we must compute the jacobian of the terms in the system (4.1). This in turn requires us to compute the partial derivatives of the f and g functions. I think the paper then goes on to explain that in order to form these partial derivative terms we consider the Jacobians of f (I will refer to these Jacobians as [1]):
J_{i} = (partial f_{i}) / (partial y)
J_{i - 1/2} = (partial f_{i - 1/2}) / (partial y)
We then compute these using finite difference to then use in the global Jacobian (printed in the middle of page 7):
(partial phi_{i}) / (partial y_{i})
Now in order to form the Jacobian terms [1] the paper says they use the function numjac to calculate the partial derivatives. I believe Octave does not have such a function...
Does anyone know of any equivalent method in octave to numjac? Do you think this is something I am going to have to implement first?
Luke
|
|
|
Here are a few thoughts on numjac.
numjac implements an algorithm for finite difference (FD) step size selection based on this paper by Salane:
D. E. Salane, Adaptive routines for forming Jacobians numerically, Tech. Report SAND86 1319, Sandia National Laboratories, Albuquerque, NM, 1986.
I have done some experimentation with this algorithm and numjac and am much less positive on this approach than either Salane or the numjac authors; I have seen cases where the algorithm prescribes very tiny step sizes that lead to large roundoff errors.
A second part of numjac is that, for sparse matrices (which I assume applies to the bvp4c algorithm) it is usually necessary to make far fewer than n (number of equations) calls to the residual function when evaluating the jacobian by FD. This substantially improves performance. A discussion of this can be found in this paper:
T. F. Coleman, B. S. Garbow, and J. J. More, Software for estimating sparse Jacobian matrices, ACM Transactions on Mathematical Software, 11 (1984), pp. 329-345.
Unfortunately, implementing this is not trivial.
A possible alternative to both of these issues is the approach used by the algorithm for calculation of jacobians by FD when the matrix is banded (I'm assuming this applies to bvp4c). The routine is idaDlsBandDQJac and it is discussed somewhat in the Sundials documentation. A downside of this (other than the implementation being in C) is that, as far as I am aware, Octave (and MATLAB) has no specific support for banded matrices.
|
|
|
I have been using Sundials and (indirectly) idaDlsBandDQJac for the past several months and it seems to beworking fine. If a band-matrix representation is acceptable in your bvp4c implementation, leveraging Sundials might be a good idea.
For the next phase of my work, I need a FD jacobian routine for general sparse matrices and, unfortunately, Sundials doesn't have this. So I've been looking further at the sparse jacobian paper and implementation by Coleman, et. al. Since my last post, I discovered that this approach
If I find out more, I'll pass it along.
Bill
|
|
|
On 6 Aug 2016, at 15:38, lakerluke < [hidden email]> wrote:
> Now in order to form the Jacobian terms [1] the paper says they use the
> function numjac to calculate the partial derivatives. I believe Octave does
> not have such a function...
>
> Does anyone know of any equivalent method in octave to numjac? Do you think
> this is something I am going to have to implement first?
There are functions to compute numerical jacobians in the optim package,
I believe the code of one such functions is replicated in ode23s in the
package odepkg to avoid placing a dependency of odepkg on optim.
I highly recommend having an m-file only implementation of bvp4c, at
least as a first step.
c.
|
|
|
As a proof-of-concept, I did a quick and dirty C++ implementation of bvp4c. One of the main purposes was to test the Sundials/Kinsol nonlinear algebraic solver for this application. In that regard, I would say my experience has been mixed. My impression is that, even with line search enabled, it is not as effective at coping with poor initial guesses as MATLAB bvp4c.
This implementation computes a dense jacobian by finite difference (FD). I had hoped to try FD computation of a sparse jacobian but have not done that yet. But even this dense jacobian code has better performance than I expected.
The main limitation (that I know of) of my little test code is that it doesn't compute discretization errors and refine the mesh; it simply computes a solution with the input mesh.
If you have any interest in looking at the code, it is here:
Bill Greene
|
|
|
This post was updated on .
Thanks Bill, that is very useful.
I have another silly question to ask. I'm failing to understand a fundamental element related to implementing bvp4c from the Shampine paper.
We form the system [4.1] (pg. 6) by applying the Simpson method to the original system. We then proceed to solve this by use of Newton's method to find the roots of the system [4.1]. This leads to an iterative method where we use the value of the solution at the boundary to form an approximation to the solution at the next value of the mesh.
bvp4c_issue1.jpgI have written this up on some paper to hopefully better explain my question. If someone could correct my understanding, that would be much appreciated.
Best regards,
Luke
|
|
|
Glad the code was useful, Luke.
Not sure I understand your question, though. Yes, equation 4.1 is the key to the algorithm and, yes, it needs to be solved for Y by some Newton-like method.
The user is required to provide the initial values for *all* the components in Y as solinit.y. That initial guess and the function defining Phi is all the KINSOL solver needs to find the Y that satisfies 4.1.
Bill
|
|
|
Bill Greene-3 wrote
The user is required to provide the initial values for *all* the
components in Y as solinit.y. That initial guess and the function defining
Phi is all the KINSOL solver needs to find the Y that satisfies 4.1.
Bill
Right, I see. So, for example, in the simple SHM ODE I pictured in my previous post, only the derivative of y is known at the initial value, I must therefore, as the user, provide a guess for the entire system at all mesh values. In matlab, this is typically performed by the helper function bvpinit.
However, in forming the global jacobian in 4.1, we need to evaluate f at mid point mesh values, x_{i - 1/2}. Does this mean we also need to provide a guess for the solution at these mid points? How else can we find the value of f_{i - 1/2}?
|
|
|
The last equation on page 6 in the paper shows how to evaluate f at the mid-point ( f_{i - 1/2}). It includes an expression for y_{i - 1/2}. They call this a "condensed" formulation.
|
|
|
The last equation on page 6 gives f_{i - 1/2} in terms of f evaluated at x_{i - 1} + h/2 .
Does this not mean that we need an expression for f evaluated at this mid way point between mesh-values?
|
|
|
Not sure what you mean by "expression". Implementation-wise, f(...) means call the user's odefun with the arguments in (). So that is what the code does with the last two equations of page 6. Similarly, g(...) means call the user's bcfun.
Bill
|
|
|
Maybe I've misunderstood this but from page 6, f_{i - 1/2} is given by odefun evaluated at the point (*):
[ (x_{i - 1} + h/2) , 0.5*(y_{i - 1} + y_{i}) - (h / 8) * (f_{ i } - f_{ i - 1}) ]
However,in order to evalue f at this point requires us to know how to evaluate the components of f at x = (x_{i - 1} + h/2), which we do not know.
Take for example demo 1 given here:
http://octave.sourceforge.net/odepkg/function/bvp4c.htmlf = @(t,u) [ u(2); -abs(u(1)) ];
Seeing as u is a function of t, in order to evaluate f at (*) do we not need to know u as a function of t (e.g. the solution) in order to evaluate it? Or, seeing as in this example, f does not explicitly depend on t, do we just evaluate the u components of f at the second component of the point (*) - effectively treating f = f(u)?
I hope this makes sense. My point essentially comes down to my interpretation of evaluating f at the point (*). I would think that in order to do this we must know the solution, u(t), in order to evaluate u at point (x_{i - 1} + h/2).
Bill Greene-3 wrote
Not sure what you mean by "expression". Implementation-wise,
f(...) means call the user's odefun with the arguments in ().
So that is what the code does with the last two equations of page 6.
Similarly, g(...) means call the user's bcfun.
Bill
On Mon, Aug 22, 2016 at 12:33 PM, lakerluke < [hidden email]>
wrote:
> The last equation on page 6 gives f_{i - 1/2} in terms of f evaluated at
> x_{i
> - 1} + h/2 .
>
> Does this not mean that we need an expression for f evaluated at this mid
> way point between mesh-values?
>
>
>
> --
> View this message in context: http://octave.1599824.n4.
> nabble.com/Working-on-bvp4c-tp4677540p4679405.html
> Sent from the Octave - Maintainers mailing list archive at Nabble.com.
>
>
|
|
|
I have recently added the calculation of the residual error and the capability to adaptively refine the mesh based on this error measure. I believe that now my code basically corresponds to the algorithm described in the Kierzenka/Shampine paper.
If there is someone out there who has some experience with the bvp4c function, has access to MATLAB, and could do some comparisons between bvp4c and this code, I would be much-obliged.
If you look at "Release 1" on this page, there are prebuilt binaries for Octave 4.0 on windows, macos, and linux. Alternatively, the Makefile I used to build the code for these three platforms is also included with the source code. (I happen to call the function bvp1d to distinguish it from bvp4c but the intent is that is closely corresponds to bvp4c.)
Thanks.
Bill Greene
|
12
|
Locked
