Solving differential equation in Python with variable coefficients (I just know the coefficients numerically)

Question

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?

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.