patchGNU Octave - Patches: patch #10023, bvp4c and bvpinit


patch #10023: bvp4c and bvpinit

Submitter:  None
Submitted:  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:  -email is unavailable-
Open/Closed:  Open
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Mon 09 Aug 2021 02:14:56 PM UTC, comment #6: 

odepkg would be the logical candidate, but AFAIK it hasn't been maintained for several years

Nir Krakauer <nir_krakauer>
Sun 08 Aug 2021 05:45:23 PM UTC, comment #5: 

Well there are several other "core Matlab function" implementations currently living in Octave packages, so maybe the functions here can start their career in an OF package? (just curious)

Philip Nienhuis <philipnienhuis>
Group Member
Sun 08 Aug 2021 01:25:39 PM UTC, comment #4: 

bvp4c (along with bvpinit) is a core Matlab function currently not in Octave. So it belongs in core Octave.

Nir Krakauer <nir_krakauer>
Sun 08 Aug 2021 08:33:59 AM UTC, comment #3: 

Which OF package do you think should this go into?
Then we can notify the pertinent package maintainer. Otherwise this patch can sit here for ages.

Philip Nienhuis <philipnienhuis>
Group Member
Mon 22 Feb 2021 09:57:45 AM UTC, comment #2: 

Fixed a typo in the finite element method which would cause a crash if the number of equations does not equal 2.

(file #50905)

Fri 19 Feb 2021 08:44:57 AM UTC, comment #1: 

As stated in the original submission I have now implemented bvp4c (bpv4cfe.m) using the finite element method. The function is completely vectorized and I get similar performance to the finite difference method originally submitted.

The implication is the standard finite element method. I used a cubic Hermite spline as the interpolation function which makes the method a third-order method rather than a fourth-order method. I assume this choice makes the interpolated values returned with interp1 using "chip" exact to what the solver returns. The function and first derivative are continuous everywhere. The type and order of the interpolation function could be changed with quite a lot of effort as most of the complexity of the code is to avoid loops, specifically, the assembly procedure would need to be rewritten. (The looped version was unbearably slow like the old Lobatto implementation.)

I integrate each element using a five-point Gauss Legendre quadrature rule. This rule integrates any ninth or lower-order polynomial exactly. The number of points can be changed by changing the variable nip. I calculate the weights and positions each time the function is run. This could be hardcoded rather than calculated if we agree the integration accuracy is appropriate.

I use the same solving and error estimating method I used for the finite difference method.

Let me know which of the methods:

  • Finite Element,
  • Finite Difference, or
  • neither if we want to use Lobatto IIIA method

to continue developing. If we choose the Lobatto IIIA method I think will need to be implemented in c++ as I cannot see a way to vectorizes it.

(file #50882)

Wed 17 Feb 2021 02:03:33 PM UTC, original submission:  

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.



(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #50905:  bvp4cfe.m added by None (7KiB - text/plain)
file #50882:  bvp4cfe.m added by None (5KiB - text/plain)
file #50871:  bvp4cfd.m added by None (4KiB - text/plain)
file #50872:  bvpinit.m added by None (849B - text/plain)
file #50873:  bvp4cold.m added by None (6KiB - text/plain)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by nir_krakauer (Posted a comment)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.


    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-02-22 None Attached File- Added bvp4cfe.m, #50905
    2021-02-19 None Attached File- Added bvp4cfe.m, #50882
    2021-02-17 None Attached File- Added bvp4cfd.m, #50871
        Attached File- Added bvpinit.m, #50872
        Attached File- Added bvp4cold.m, #50873

    Back to the top

    Powered by Savane 3.13-c54d.
    Corresponding source code