payload: {"action":"created","issue":{"url":"https://api.github.com/repos/scipy/scipy/issues/9832","repository_url":"https://api.github.com/repos/scipy/scipy","labels_url":"https://api.github.com/repos/scipy/scipy/issues/9832/labels{/name}","comments_url":"https://api.github.com/repos/scipy/scipy/issues/9832/comments","events_url":"https://api.github.com/repos/scipy/scipy/issues/9832/events","html_url":"https://github.com/scipy/scipy/issues/9832","id":411234771,"node_id":"MDU6SXNzdWU0MTEyMzQ3NzE=","number":9832,"title":"integrate.solve_bvp does not check for convergence of boundary conditions","user":{"login":"scbarton","id":1451344,"node_id":"MDQ6VXNlcjE0NTEzNDQ=","avatar_url":"https://avatars3.githubusercontent.com/u/1451344?v=4","gravatar_id":"","url":"https://api.github.com/users/scbarton","html_url":"https://github.com/scbarton","followers_url":"https://api.github.com/users/scbarton/followers","following_url":"https://api.github.com/users/scbarton/following{/other_user}","gists_url":"https://api.github.com/users/scbarton/gists{/gist_id}","starred_url":"https://api.github.com/users/scbarton/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/scbarton/subscriptions","organizations_url":"https://api.github.com/users/scbarton/orgs","repos_url":"https://api.github.com/users/scbarton/repos","events_url":"https://api.github.com/users/scbarton/events{/privacy}","received_events_url":"https://api.github.com/users/scbarton/received_events","type":"User","site_admin":false},"labels":[{"id":37334597,"node_id":"MDU6TGFiZWwzNzMzNDU5Nw==","url":"https://api.github.com/repos/scipy/scipy/labels/defect","name":"defect","color":"e11d21","default":false},{"id":37342805,"node_id":"MDU6TGFiZWwzNzM0MjgwNQ==","url":"https://api.github.com/repos/scipy/scipy/labels/scipy.integrate","name":"scipy.integrate","color":"207de5","default":false}],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":13,"created_at":"2019-02-17T21:03:37Z","updated_at":"2019-03-04T14:43:24Z","closed_at":null,"author_association":"NONE","body":"I’d like to report an issue with `scipy.integrate.solve_bvp`. This is an excellent piece of code but it unfortunately does not check for convergence of boundary conditions. I can demonstrate with the following example code with nonlinear boundary conditions. The code solves Laplace's equation, which is linear and converges in only one iteration.\r\n\r\n### Reproducing code example:\r\n```\r\nimport numpy as np\r\nfrom scipy.integrate import solve_bvp\r\n\r\n# Parameters\r\nkappa= 1.64\r\nioA, ioC= 0.01, 1.0e-4\r\nV, UA, UC= 0.5, 0.0, 1.0\r\nf= 38.9\r\nL=0.1\r\n\r\n## Governing Equation- Laplace Eq. for potential\r\ndef deq(x,y):\r\n    return np.stack([y[1], np.zeros_like(x)])\r\n\r\n## Boundary Conditions\r\ndef bc(uA,uC):\r\n\r\n    phiA, phipA= uA\r\n    phiC, phipC= uC\r\n\r\n    # Butler-Volmer Kinetics at Anode\r\n    hA=0-phiA-0\r\n    iA= ioA * ( np.exp(f*hA) - np.exp(-f*hA))\r\n    res0= iA + kappa*phipA\r\n\r\n    # Butler-Volmer Kinetics at Cathode\r\n    hC= V - phiC - UC\r\n    iC= ioC * (np.exp(f*hC) - np.exp(-f*hC))\r\n    res1= iC - kappa*phipC\r\n\r\n    return np.array([  res0, res1 ] )\r\n        \r\n# Initial Guess\r\nxinit=np.linspace(0,L)\r\nuinit=np.array([0,0]) [:,None]   # a column vector\r\nyinit=uinit * np.ones([len(uinit),len(xinit)])\r\n\r\n# Find Solution\r\nsol = solve_bvp(deq, bc, xinit, yinit, verbose=2)\r\n\r\n# Calculate Boundary residuals\r\nbc_res = bc(sol.y[:, 0], sol.y[:, -1])\r\nprint('Boundary Condition Residuals: ',bc_res)```\r\n```\r\n\r\nUsing scipy v1.2.0, the above code generates the following output:\r\n\r\n```\r\n   Iteration    Max residual    Total nodes    Nodes added  \r\n       1          3.39e-14          50              0       \r\nSolved in 1 iterations, number of nodes 50, maximum relative residual 3.39e-14.\r\nBoundary Condition Residuals:  [ 5.25230866e-02 -2.30274472e+02]\r\n```\r\n\r\nAfter only one iteration, the BC residuals are much greater than the tolerance (the default tolerance, 1e-3).\r\n\r\n### Proposed solution\r\nI have added code to check for convergence, modifying lines 1081-1111 of `_bvp.py` as follows:\r\n\r\n```\r\n        bc_res = bc_wrapped(y[:, 0], y[:, -1], p)  \r\n        \r\n        # This relation is not trivial, but can be verified.\r\n        r_middle = 1.5 * col_res / h\r\n        sol = create_spline(y, f, x, h)\r\n        rms_res = estimate_rms_residuals(fun_wrapped, sol, x, h, p,\r\n                                         r_middle, f_middle)\r\n        max_rms_res = np.max(rms_res)\r\n        max_bc_res = np.max(np.abs(bc_res))\r\n        max_res = np.max([max_rms_res, max_bc_res])\r\n\t\t\r\n        if singular:\r\n            status = 2\r\n            break\r\n\r\n        insert_1, = np.nonzero((rms_res > tol) & (rms_res < 100 * tol))\r\n        insert_2, = np.nonzero(rms_res >= 100 * tol)\r\n        nodes_added = insert_1.shape[0] + 2 * insert_2.shape[0]\r\n\r\n        if m + nodes_added > max_nodes:\r\n            status = 1\r\n            if verbose == 2:\r\n                nodes_added = \"({})\".format(nodes_added)\r\n                print_iteration_progress(iteration, max_res, m,\r\n                                         nodes_added)\r\n            break\r\n\r\n        if verbose == 2:\r\n            print_iteration_progress(iteration, max_res, m, nodes_added)\r\n\r\n        if nodes_added > 0:\r\n            x = modify_mesh(x, insert_1, insert_2)\r\n            h = np.diff(x)\r\n            y = sol(x)\r\n        elif max_res <= tol:\r\n            status = 0\r\n            break\r\n```\r\n\r\nThis code calculates `max_res` for both the domain and the boundary, and only breaks the parent `while` loop when that `max_res` falls below `tol`.\r\n\r\nUsing this version of `solve_bvp`, the same Laplace Equation code results in multiple iterations to reduce the boundary residuals to less than the tolerance.\r\n\r\n```\r\n   Iteration    Max residual    Total nodes    Nodes added  \r\n       1          2.30e+02          50              0       \r\n       2          2.86e+00          50              0       \r\n       3          6.18e-03          50              0       \r\n       4          7.25e-06          50              0       \r\nSolved in 4 iterations, number of nodes 50, maximum relative residual 7.25e-06.\r\nBoundary Condition Residuals:  [ 2.04326239e-07 -7.25180010e-06]\r\n\r\n```\r\n\r\nI've attached the modified verion of `_bvp.py`, and I'll be glad to make a pull request if needed.\r\n\r\n### Scipy/Numpy/Python version information:\r\n\r\n```\r\nimport sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)\r\n\r\n1.2.0 1.16.1 sys.version_info(major=3, minor=7, micro=2, releaselevel='final', serial=0)\r\n```\r\n[_bvp-scb.py.txt](https://github.com/scipy/scipy/files/2873413/_bvp-scb.py.txt)\r\n"},"comment":{"url":"https://api.github.com/repos/scipy/scipy/issues/comments/469276911","html_url":"https://github.com/scipy/scipy/issues/9832#issuecomment-469276911","issue_url":"https://api.github.com/repos/scipy/scipy/issues/9832","id":469276911,"node_id":"MDEyOklzc3VlQ29tbWVudDQ2OTI3NjkxMQ==","user":{"login":"scbarton","id":1451344,"node_id":"MDQ6VXNlcjE0NTEzNDQ=","avatar_url":"https://avatars3.githubusercontent.com/u/1451344?v=4","gravatar_id":"","url":"https://api.github.com/users/scbarton","html_url":"https://github.com/scbarton","followers_url":"https://api.github.com/users/scbarton/followers","following_url":"https://api.github.com/users/scbarton/following{/other_user}","gists_url":"https://api.github.com/users/scbarton/gists{/gist_id}","starred_url":"https://api.github.com/users/scbarton/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/scbarton/subscriptions","organizations_url":"https://api.github.com/users/scbarton/orgs","repos_url":"https://api.github.com/users/scbarton/repos","events_url":"https://api.github.com/users/scbarton/events{/privacy}","received_events_url":"https://api.github.com/users/scbarton/received_events","type":"User","site_admin":false},"created_at":"2019-03-04T14:43:24Z","updated_at":"2019-03-04T14:43:24Z","author_association":"NONE","body":"The second strategy passes all of the tests, which as you know includes several problems.  I did need to lower the tolerance on my new test, `test_nonlin_bc`, to be equal to the assertion tolerance of 1e-5 on `sol(x)`.\r\n\r\nDo you think we need to test against more problems? I can think of a few."}}
other: {"actor":{"display_login":"scbarton"}}