1
$\begingroup$

I am trying to solve a boundary value problem on $[0, \infty]$, using scipy's scipy.integrate.solve_bvp and I am seeing that the solutions are not converging even with high number of iterations, and that they depend strongly on my initial guess.

The problem I have is complex but persists even in the following minimal example. Consider

$$(rh')' = 4r h^3$$

with boundary conditions $h(0)= h(r_\mathrm{max})=0$ for some finite $r_\mathrm{max}$. Clearly $h=0$ is a solution. However,scipy.integrate.solve_bvp fails to find this trivial solution for generic initial guesses. For example:

from scipy.interpolate import interp1d

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy.integrate
matplotlib.rc("mathtext",fontset="cm") 

####################################################################

def solve_h(tol=1e-7,nb_rvals=500000,max_rval=100,max_nodes=1000000):

    x = np.logspace(-7,0,nb_rvals)*max_rval

    def fun(r,y):
        return y[1]/r, r*(4*y[0]**3)

    def bc(ya,yb):
        return ya[0], yb[0]
    
    # guess for vector y at points x  
    y=[x, np.zeros(len(x))]

    #solve
    res=scipy.integrate.solve_bvp(fun, bc, x, y, tol=tol, max_nodes=max_nodes)
    print(res.message,"  ", len(res.x))

    ##################################################
    # compute derivatives of radial mode profile
    h = res.sol(x)[0]
    rms_residuals = res.rms_residuals
    
    return x, h, rms_residuals

x, h, rms_residuals = solve_h(tol=1e-3,nb_rvals=500000,max_rval=10,max_nodes=1000000)

plt.plot(x,h)

I am interested to know why the solver is not converging, and if I should instead be using some more numerically stable method. On the second point, the problem I'm trying to solve actually does not have Dirichlet boundary conditions but instead mixed conditions like $h'(0)=0, h(\infty)=0$,

$\endgroup$

1 Answer 1

3
$\begingroup$

You should use the S parameter for the singular part of the equation. With $g=h'$ the equation can be written as system $$ h''+\frac{h'}r=4h^3\\~\\ \pmatrix{h'\\g'}=\underbrace{\pmatrix{0&0\\0&-1}}_{=S}\pmatrix{h/r\\g/r}+\pmatrix{g\\4h^3} $$ The initial guess should not be too wild, it should come close to both boundary conditions. Remember that the internal step of the solver is the solution of a high-dimensional non-linear system resulting from the collocation method that is used as multiple-shooting method. Try $$h(r)=r·e^{-r}, ~~h'(r)=(1-r)e^{-r}.$$

The solution that crosses the axis at $r_\max$ is likely to diverge after that explosively, meaning it is far from the asymptotic solution. This will also lead to strange solution behavior before that point. See " Boundary value problem with singularity and boundary condition at infinity " for an example of spurious wild oscillations.

At $r=0$ you have a regular singularity with a double root at zero for the characteristic polynomial. This give a power series for the finite solution. With the limit of the equation at $r\to 0$ giving $h'(0)=4h(0)^3$ the solution starts as $h(r)=h_0+h_0^3r^2+...$ at $r\approx 0$. Indeed one would expect it to increase fast enough that the negative term $-h'/r$ never becomes relevant for the shape of the curve.

For large $r$ the equation can be approximated as $h''=4h^3$. Multiplying with $2h'$ and integrating gives $$ h'^2=2h^4+C $$ with $C=0$ for an asymptotic convergence to zero. For the solutions converging to zero the suitable square root is $$ h'=-\sqrt2|h|h, $$ which also can be solved to $$ h=\frac{a}{1+br},~~~ h'=-\frac{ab}{(1+br)^2}=-\frac{b}{a}h^2, ~~\frac{b}{a}=\sqrt2·{\rm sign}(h). $$

Assembled into code this gives

S = [[0,0],[0,-1]]
def f(r,u): h,g = u; return [g, 4*h**3]
def bc(u0, uR): return [ u0[0], uR[1]+2**0.5*abs(uR[0])*uR[0] ]

R=100; r = np.logspace(-3,0,15)*R; r[0]=0; u = [ r*np.exp(-r), (1-r)*np.exp(-r) ]
res = solve_bvp(f,bc,r,u,S=S,tol=1e-4); res.message
#>>> 'The algorithm converged to the desired accuracy.'
max(abs(res.y[0]))
#>>> 1.1812227123833282e-08
##  with tol=1e-6 the max is 1.0067373977626852e-12

This is the zero solution, within the bounds of the error tolerance.

$\endgroup$
1
  • $\begingroup$ thanks. unfortunately your suggested changes don't solve the more complex problem I have, so this minimal example does not accurately capture the issues. I will open a new question with the more complex problem and tag you if you could take a look! $\endgroup$
    – phyphy
    Oct 9, 2022 at 2:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.