BVP5C Nmax setting and the number of mesh points of the solution

I am using the BVP5C function to solve a nonlinear ODE BVP (consisting of 4 first-order ODEs and 1 unknown parameter). To aid the solver, I supply analytical partial derivatives of the vector field as well as the BCs (with respect to both the four state variables and the unknown parameter). I am interested in high accuracy solution (RelTol=100*eps()). I seek a continuation solution with respect to one of the input state variables, say, theta, i.e., I am solving the BVP at each value of theta as it varies in an interval [0, Pi]. While my code works for most smaller values of theta, the accuracy breaks down as theta increases beyond 2. I have increased the number of mesh points by setting 'Nmax=100000'. However, the actual number of mesh points in the output is much less than that. See the warning that I get:

Warning: Unable to meet the tolerance without using more than 100000 mesh points. The last mesh of 16892 points and the solution are available in the output argument. The maximum error is 2.59365e-08, while requested accuracy is 2.22045e-14.

I am wondering why did it not give me the solution with 100000 mesh points. Do I need to increase Nmax even further?

Another thing I am wondering is MATLAB's bvpset doc says the default Nmax value is floor(1000/n), which in my case would be 250. This does not make any sense since I have obtained a valid solution (i.e., meeting RelTol requirement) which had more than 250 mesh points when using the default Nmax setting for my problem (at smaller values of theta).

[Saurabh Das]

Hi Navendu,

Do the usage of analytical partial derivatives help to increase the ACCURACY of the solution or does it just help with reducing computational time?

Thanks,

Saurabh

Answer 1

[Benjamin Larin]

I have encountered a very similar issue and obtained the following insights.

First, the documentation for bvpset is unspecific or wrong. It always refers to bvp4c and it is not clear whether all options and descriptions of default settings apply equally to bvp5c. In case of NMax the default is actually floor(10000/n), as can be seen in bvparguments.m, line 148. In your case that would imply a value of 2500 instead of 250. The documentation should be adjusted accordingly.

Second, I guess that the algorithm tries a solution with an (endogenously chosen) larger mesh, but stops, if this is greater than NMax and therefore returns the last solution where the upper bound for the mesh size is not violated. One could take a deeper look at bvp4c.m in order to fully grasp that.