CODEWITHSUNDEEP
pythonnumpymathscipydifferential-equations
rhombidodecahedron
rhombidodecahedron

Reputation: 7922

Solving a BVP on a fixed non-uniform grid in python without interpolating

I am aware of scipy.solve_bvp but it requires that you interpolate your variables which I do not want to do.

I have a boundary value problem of the following form:

y1'(x) = -c1*f1(x)*f2(x)*y2(x) - f3(x)

y2'(x) = f4(x)*y1 + f1(x)*y2(x)

y1(x=0)=0, y2(x=1)=0

I have values for x=[0, 0.0001, 0.025, 0.3, ... 0.9999999, 1] on a non-uniform grid and values for all of the variables/functions at only those values of x.

How can I solve this BVP?

Upvotes: 1

Views: 633

Answers (1)

hpaulj
hpaulj

Reputation: 231530

This is a new function, and I don't have it on my scipy version (0.17), but I found the source in scipy/scipy/integrate/_bvp.py (github).

The relevant pull request is https://github.com/scipy/scipy/pull/6025, last April.

It is based on a paper and MATLAB implementation, 

       J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
       Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
       Number 3, pp. 299-316, 2001.

The x mesh handling appears to be:

while True:
    ....
    solve_newton
    ....
    insert_1, = np.nonzero((rms_res > tol) & (rms_res < 100 * tol))
    insert_2, = np.nonzero(rms_res >= 100 * tol)
    nodes_added = insert_1.shape[0] + 2 * insert_2.shape[0]

    if m + nodes_added > max_nodes:
        status = 1
        if verbose == 2:
            nodes_added = "({})".format(nodes_added)
            print_iteration_progress(iteration, max_rms_res, m,
                                     nodes_added)
    ...
    if nodes_added > 0:
        x = modify_mesh(x, insert_1, insert_2)
        h = np.diff(x)
        y = sol(x)

where modify_mesh add nodes to x based on:

insert_1 : ndarray
    Intervals to each insert 1 new node in the middle.
insert_2 : ndarray
    Intervals to each insert 2 new nodes, such that divide an interval
    into 3 equal parts.

From this I deduce that

  • you can track the addition of nodes with the verbose parameter

  • nodes are added, but not removed. So the out mesh should include all of your input points.

  • I assume nodes are added to improve resolution in certain segments of the problem

This is based on reading the code, and not verified with test code. You may be the only person to be asking about this function on SO, and one of the few to have actually used it.

Upvotes: 1

Related Questions