Working on bvp4c

classic Classic list List threaded Threaded
23 messages Options
12
Reply | Threaded
Open this post in threaded view
|

Working on bvp4c

lakerluke
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
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Søren Hauberg
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


Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
In reply to this post by lakerluke
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
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Carlo de Falco-3

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.




Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
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
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Carlo de Falco-3

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.




Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
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
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Bill Greene-3
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.

On Sat, Aug 6, 2016 at 9:38 AM, lakerluke <[hidden email]> wrote:
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



--
View this message in context: http://octave.1599824.n4.nabble.com/Working-on-bvp4c-tp4677540p4679048.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.


Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
Thank you for the summary and information Bill. I have been looking into Sundials via the link you mentioned. I couldn't actually find a reference to idaDlsBandDQJac (presumably either in IDA or IDAS) in the following user guides linked to from here:

http://computation.llnl.gov/projects/sundials-suite-nonlinear-differential-algebraic-equation-solvers/sundials-software

However I did find the code here:

http://sundials.sourcearchive.com/documentation/2.4.0/idas__direct_8c-source.html

Given your background research do you think the best way to write numjac is via the Sundials approach? Did you get anywhere with this?

Luke
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Bill Greene-3
I have been using Sundials and (indirectly) idaDlsBandDQJac for the past several months and it seems to be
working 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 Sun, Aug 7, 2016 at 2:28 PM, lakerluke <[hidden email]> wrote:
Thank you for the summary and information Bill. I have been looking into
Sundials via the link you mentioned. I couldn't actually find a reference to
idaDlsBandDQJac (presumably either in IDA or IDAS) in the following user
guides linked to from here:

http://computation.llnl.gov/projects/sundials-suite-nonlinear-differential-algebraic-equation-solvers/sundials-software

However I did find the code here:

http://sundials.sourcearchive.com/documentation/2.4.0/idas__direct_8c-source.html

Given your background research do you think the best way to write numjac is
via the Sundials approach? Did you get anywhere with this?

Luke



--
View this message in context: http://octave.1599824.n4.nabble.com/Working-on-bvp4c-tp4677540p4679058.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.


Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Carlo de Falco-3
In reply to this post by lakerluke

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.





Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Bill Greene-3
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 

On Fri, Aug 12, 2016 at 6:56 PM, c. <[hidden email]> wrote:

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.






Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
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.jpg

I 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
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Bill Greene-3
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

On Sat, Aug 20, 2016 at 12:43 PM, lakerluke <[hidden email]> wrote:
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.

I have written this up on some paper to hopefully better explain my
question. If someone could correct my understanding, that would be much
appreciated.

<http://octave.1599824.n4.nabble.com/file/n4679349/bvp4c_issue1.jpg>

Best regards,

Luke



--
View this message in context: http://octave.1599824.n4.nabble.com/Working-on-bvp4c-tp4677540p4679349.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.


Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
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}?
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Bill Greene-3
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.

On Mon, Aug 22, 2016 at 3:44 AM, lakerluke <[hidden email]> wrote:
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}?



--
View this message in context: http://octave.1599824.n4.nabble.com/Working-on-bvp4c-tp4677540p4679376.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.


Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
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?
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Bill Greene-3
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.


Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

lakerluke
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.html

f = @(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.
>
>
Reply | Threaded
Open this post in threaded view
|

Re: Working on bvp4c

Bill Greene-3
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.

The source code is here: https://github.com/wgreene310/bvp1d.
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

On Tue, Aug 23, 2016 at 5:29 AM, lakerluke <[hidden email]> wrote:
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.html

f = @(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 <

> luke_purnell@

> >
> 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.
>>
>>





--
View this message in context: http://octave.1599824.n4.nabble.com/Working-on-bvp4c-tp4677540p4679432.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.


12