-2
$\begingroup$

Currently, I am trying to apply a Crank-Nicolson method on a function that I want to evolve. However, I am only facing one drawback and it is the step when I normalize the initial function using the scipy integrate function.

It tells me that it runs into a

runtime error: invalid value encountered in dive

This only happens when I change the sigma parameter in the function (it works when sigma=0.05 but I require to be sigma=5e-7). If anyone could give me some insight about this I will appreciate it. Here is the relevant part of the code:

import scipy as sp
import numpy as np
from scipy import integrate, sparse, linalg,special
import scipy.sparse.linalg
import pylab as pl
from multiprocessing import Pool
import multiprocessing

def crank_nicol(tf):
    ax=-10. #punto inicial
    bx=10. #punto final
    nx= 500#numero de puntos en el x-grid
    Lx=bx-ax #longitud total
    delx= Lx/nx #x-step size=Lx/nx
    dt= 1.e-3#T/mt # time step delta t
    mt=int((2*np.pi*tf/dt)/400)
    nux= 1j*dt/(4.0*(delx**2)) #nu q son iguales a mis c_i=b_i 

    gridx = sp.zeros(nx) #space grid
    igridx = sp.array(range(nx)) # time grid
    psi = sp.zeros(nx)
    pot = sp.zeros(nx)
    gridx = delx*(igridx - nx/2)
    x=gridx
##############
    sigma=0.05
    alpha=2.66
    R_ave=2.75
    r_e=.969
    xmin=(1/2.66)*np.arccosh(0.5*np.exp(2.66*1.012/2.0))
    ####################################################
    psi = sp.pi**(-0.25)*sp.exp(-(0.5/sigma**2)*(gridx+xmin)**2) #-1j*3*(gridx-10))

    psi /= sp.sqrt(sp.integrate.simps(psi*np.conjugate(psi), dx=delx))# THIS IS THE PART WITH THE ERROR

    print(nux)


    return x, psi*np.conjugate(psi),psi

times=np.linspace(10,200,20)
if __name__ == '__main__':
    n_times = 10
    pool = multiprocessing.Pool(processes=20)
    results = pool.map(crank_nicol, times)
$\endgroup$

closed as off-topic by Kirill, GoHokies, Christian Clason, nicoguaro Jan 15 at 1:00

This question appears to be off-topic. The users who voted to close gave this specific reason:

If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ You expect to get a sensible value for exp(-1e11)? Sum up the vector of zeros and then divide by this zero result? $\endgroup$ – LutzL Jan 12 at 11:45
  • $\begingroup$ is that I don't understand why it gives me a vector of zeros when I only change the sigma value? I believe that python should be able to perform the calculationeven if the number is very small (if this is not the case, please explain me the reason why) $\endgroup$ – Oliver Jan 14 at 20:16
  • $\begingroup$ Because exp(-1e11) has about 1e10 zeros after the decimal point. The floating point representation of that result can only be 0. Do not forget, if sigma is very small, 0.5/sigma**2 is very large. $\endgroup$ – LutzL Jan 14 at 20:24