I am trying to implement a routine to solve a differential equation in Python. Basically the kind of equation that I am interested in solving is of the form:

$$\displaystyle \frac{d}{dx^2} \left(x y(x) \right) = 2 x ((U(x)-a) y(x)+ 2 b y(x)^3)$$

where $a$ is an unknown constant and $b$ is a known constant and $U(x)$ is a function that depends on $x$ but that I only know numerically (I mean, I don't have an explicit form of $U(x)$ in terms of $x$).

I need to find the value of $a$ that fulfills my initial and boundary conditions ($y(0)=y_0$ and $y(x\rightarrow \infty)=0$, $a<1$).

For that purpose, I was considering using a shooting method (a secant method to be precise) by solving several times the above equation with RK4 (using scipy.integrate.ode).

The problem that I have is that I don't know how to introduce the numerical value of $U(x)$ in my equations given the fact that ode only asks for values at the initial conditions.

Is there a way to solve my equation with scipy.integrate.ode or with another solver?

  • $\begingroup$ Could you elaborate more on your statement that ode only asks for values at the initial conditions? Looking at the documentation here, you would need to specify f(x)as an argument, not just at the initial conditions but as a function that can be evaluated at any x. More detail would definitely be helpful in answering this. $\endgroup$
    – emprice
    Commented Apr 6, 2019 at 1:22
  • $\begingroup$ I suppose that you mean about $U(x)$. Basically what I understand about how to use it is that you provide some initial conditions (at $x=x_0$) and then ode will compute a posterior step (at $x=x_0+dx$) by considering the initial condition that you provided (it will calculate the value at $x_0+dx$ by doing several internal steps until a given accuracy is obtained). Then, it will take that new value as the initial condition and it will compute the next step (at $x=x_0+2*dx$). Then, basically the problem that I have is that for $U(x)$ I don't have a functional form for that expression. $\endgroup$ Commented Apr 7, 2019 at 3:11

1 Answer 1


First of all, you need to rewrite your non-linear differential equation in a set of first order differential equations since the solve_ivp routine solves problems of the form $u'(x) = f(x,u), \ u(0) = u_{0}$ where $u(x)$ can be a vector. So define $u_{1}(x) = y(x)$ and $u_{2}(x) = y'(x)$ and you can write your original equation as \begin{array}{lll} u_{1}'(x) &=& u_{2}(x) \\ u_{2}'(x) &=& 2(U(x)-a)u_{1}(x)+4b(u_{1}(x))^{3} - \frac{2}{x}u_{2}(x)\end{array} which is in the form required.

You will need to set your initial conditions $u_1(0) = y_{0}$ and $u_{2}(0) = y'(0)$; the latter which you do not specify in your original problem statement. But looking at the second equation, it would need to be zero since otherwise the last term gets a division by zero.

Next, you would need a routine that returns $U(x)$ as a function of $x$. You say that you have no analytical form but only numerical data. So you need to write a routine that interpolates the available numerical data (or an another type of approximation you see to be fit).

Then, you could launch a series of calculations with different values of your parameter a to see when the solution tends to zero for $x\to\infty$. And wrap these in a secant solver (or any other solver) to find values for a.

  • $\begingroup$ Thank you very much for your response. Basically the step that I had no idea of how to proceed is about the routine for the interpolation for $U(x)$. I will try to continue with this new information and I will let you know if I have any problem. Again, thank you very much for your reply! $\endgroup$ Commented Apr 13, 2019 at 10:50

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.