[Octave-patch-tracker] [patch #10023] bvp4c and bvpinit

[anonymous]

URL:
  <https://savannah.gnu.org/patch/?10023>

                 Summary: bvp4c and bvpinit
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Wed 17 Feb 2021 02:03:33 PM UTC
                Category: Core : new function
                Priority: 5 - Normal
                  Status: None
                 Privacy: Public
             Assigned to: None
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any

    _______________________________________________________

Details:

I have implemented bvp4c using a finite difference method which is different
from the Lobatto IIIA methods used in Matlab.

= The Theoretical Basis of the Method =

In my method to obtain an estimate of the derivative, I use the nearest five
points to calculate the coefficients of a fourth-order polynomial passing
through the potentially unevenly spaced points. This part of the method in the
subfunction finitedifferencecoef1d follows the normal way of deriving
finite-difference rules with more than two points. It matches coefficients of
the fit with the Taylor series expansion about the point we want to take the
derivative.

I use the inbuilt solver fsolve to obtain the values of the function with a
residue calculated from my finite difference derivative minus the
user-supplied function. Boundary conditions are appended.

I estimate the error by comparing my functions derivative with the
user-supplied derivative at the midpoint between two adjacent points. If the
difference is greater than the tolerance I add that point and resolve until
all midpoints are within the required tolerance.

= Why Not to use Matlab's Method =

The Lobatto IIIA method used by Matlab is an implicit Runge–Kutta method.
Runge–Kutta methods estimate the value of the function using recursive
substitution of weighted derivatives (different locations) into the derivative
function. This recursive substitution requires loops within the residue
function which may be called many times during the solution process. This
inability to vectorize the residue will make any pure Octave implementation
very slow. I found an implementation that I think was once part of Octave's
ode package that uses the Lobatto IIIA method. This implementation is
unbearably slow.

= Comparison of the Methods =

The main advantage of my method is the implementation is almost completely
vectorized. The residual function is completely vectorized. Profiling the
runtime shows almost all of the time is spent in the solver with \ the most
time-consuming operation.

My method has no subpoints between the points returned by the method. This
results in a larger system of equations needing to be solved.

As the method uses the nearest five points to calculate the function and
derivative, the function and the derivative are not continuous.

= Current State of the Code =

The code can solve odes efficiently. Two demos included.

It is missing the following:

* Ability to apply constraints in the middle of the domain
* Solve for unknown parameters
* Use analytical derivatives

My method is bvp4cfd.m it uses bvpinit.m to prepare the initial guess. The
method I found which I think was part of the ode package is bvp4cold.m.
(bvp4cold.m is not my code and I am not suggesting to include it. It
highlights the difficulties of using the Runge–Kutta method.)

I have not prepared the documentation or followed the style guide I will fix
these issues after I finalize the code.

= Information I Would Like Before Finishing the Code =

Is it acceptable to use a different method to what Matlab uses to solve the
same problem?

I'm planning to write a finite element method for bvp4c which will allow a
continuous function and first derivative. It contains subpoints so the system
of equations the solver needs to solve should be smaller. I think vectorized
code should be possible. If I can get the finite element method to a
comparable speed to this implementation would there be a preference for finite
difference or finite element method?

Let me know if any of my input or output is different from Matlab.



    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Wed 17 Feb 2021 02:03:33 PM UTC  Name: bvp4cfd.m  Size: 4KiB   By: None

<http://savannah.gnu.org/patch/download.php?file_id=50871>
-------------------------------------------------------
Date: Wed 17 Feb 2021 02:03:33 PM UTC  Name: bvpinit.m  Size: 849B   By: None

<http://savannah.gnu.org/patch/download.php?file_id=50872>
-------------------------------------------------------
Date: Wed 17 Feb 2021 02:03:33 PM UTC  Name: bvp4cold.m  Size: 6KiB   By: None

<http://savannah.gnu.org/patch/download.php?file_id=50873>

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/patch/?10023>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/