Solving an ODE using odeint in Python and continuing the integration

Question

The following relates to the linked question:

Scattering of waves in a symmetrical potential (using python)

I have attempted to solve the problem for $U(r)$ using odeint. From this, I need to find $\mathrm{arctan}(U(r))$ and continue the integration to a large $r$

I am aware that for $\mathrm{arctan}$ I should use np.arctan

My issue is that I am getting dimensional errors with my code and hence cannot solve for $U(r)$

I am also unsure as to how to continue the integration to a large value.

I have attached my code (I am new to Python so apologies for the obvious errors.)

I really could do with help here thanks.

"""
Code to integrate ODE to numerically solve phase-shifts for given
values of k and r
"""
import numpy as np
import scipy.special as sp
import matplotlib.pyplot as plt
from scipy.integrate import odeint


def pend(U, r, l, k, A):
    theta,  omega  =  U
    dUdr = [r,  (-k*(A**2)/(r**2))([sp.spherical_jn(l, z, derivative=False) \
           -(sp.spherical_yn(l, z,  derivative=False)*U(r))])**2]
    return dUdr


# Set limits and step 
rmin = 0
rmax = 100
dr = 1 

# Set up array of r-values
r = np.arange(rmin, rmax+dr, dr);
N = len(r);

# Set Constants
A = 35.3
l = [0.0, 1.0]
k = 0.5
r = 7
z = k*r

# set initial conditions
U0  =  [0.0]

# Solve
sol = odeint(pend, U0, r, args=(l, k, A))

# Plot
plt.legend(loc='best')
plt.xlabel('r')
plt.ylabel('U(r)')
plt.grid()
plt.show()
```

Answer 1

since you are new to python, I suggest you look at the basic examples on how to use the odeint package and build your solution from there. It will take some time but your code will become more robust. That being said, there are many mistakes in the code. The error you are complaining about is due to the fact that you have assigned the vector of evaluation points to the variable r and then typed

r=7

which re-assignes r to 7 which is a scalar. The interpreter is telling you that you need a vector to compute the time step dt.

Since you have not specified the mathematical model, I will link only a pseudo-code

import numpy as np
import scipy.special as sp
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Set limits, step and array
r = np.linspace(1,100,1000)

# Set Constants
A, l, k, rc = 35.3, [0.0, 1.0], 0.5, 7
z = k * rc

def pend(U, r):
    global A, l, rc, k
    theta, omega = U
    d_theta = ...
    d_omega = ...
    return [d_theta, d_omega]

sol = odeint(pend, [0.0,0.0], r)

I think you should use this structure and add the expression of the derivative. When you write the function to integrate you should remember that it takes two parameters, U and r. r is the point at which you are evaluating the function U, both U and r are scalars. Typing something like U(r) makes no sense since the variable U is already U(r).

Also note that you are working with vectors and lists (i.e. simple [1,2,3]) do not implement the operations you where trying to perform. I suggest you either use numpy or calculate the derivatives of the components like I did in the pseudo-code. I suggest that you only perform scalar operations, so express everything like it. They are easier to debug and modify.

As far as the integration is concerned, it really depends on the model you are trying to discretize. If you know that for larger values of r the solution becomes more smooth then you can take larger and larger steps (i.e. build a vector with dr large). Some models require a change of variable, in your case it could be 1/r.