Working on bvp4c

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

Re: Working on bvp4c

On Sat, Aug 27, 2016 at 2:11 PM, Bill Greene <[hidden email]> wrote:
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.

It may be more feasible to post sample scripts/ commands for someone to run on Matlab and report back the results. (either a cut/paste of output, or a mat file, or something that you could compare). 

Reply | Threaded
Open this post in threaded view

Re: Working on bvp4c

In reply to this post by lakerluke

On 23 Aug 2016, at 11:29, 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:
> 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)?

If I understand correctly what formula you refer to, it is the one on line 5 at page 306 of the article.

In that formula, f is a function of two inputs, both of which are computable at the current step in the algorithm so I really don't understand what your question is about.

Maybe it helps if you consider the formula programmatically, i.e. divide the evaluation of f_{i - 1/2} into 5 steps:

1) f_{i}       = f (x_{i}, y_{i}, p)

2) f_{i - 1}   = f (x_{i - 1}, y_{i - 1}, p)

3) arg1        = (x_{i - 1} + h_{i - 1} / 2)

4) arg2        = 0.5 * (y_{i - 1} + y_{i}) - ...
                 (h_{i - 1} / 8) * (f_{i} - f_{i - 1});

5) f_{i - 1/2} = f (arg1, arg2, p)

at each of the 5 steps all of the quantities appearing on the RHS are known so the quantitities on the LHS are computable.
Which of the above steps is unclear to you?


Reply | Threaded
Open this post in threaded view

Re: Working on bvp4c

Yes, that is the correct formula. I think it is clear now.

However, I still want to check that I understand how the algorithm works programmatically before coding anything up. To this end I want to run through how the algorithm will operate on a simple example. This is what I have attempted to do with the two attached files in order to form the first approximation term. I would really appreciate it if you could see if this looks like the correct idea.

[1] - Set up the first order system, calculate domain, provide guess for solution.


[2] - Plug in to the iterative formula to find the first non trivial term of Phi. Then use Newton's method to calculate next term in approximation to solution.