Problem using integration scheme in numpy

Question

I am trying to perform integration of a function in order to find the mean position. However when I perform the integrating with quad I get problems of dimensions not matching. When I run the function on its own it works without a problem. However, when the function is used by the quad integration scheme it gives the error of dimenions mismatch. I will post the completely functional code below and the error message and I hope someone can tell me whats going wrong and how I can fix it.

Please let me know if anything is unclear so I can add more information.

import numpy as np
import time
from scipy.sparse.linalg import eigsh
from scipy.sparse.linalg import spsolve
from scipy.sparse import diags
from scipy.sparse import identity        
import scipy.sparse as sp        
from scipy.integrate import quad
x0     = 15
v      = 16
sigma2 = 5
tmax   = 4       
        
def my_gaussian(x, mu=x0, var=5):
    
    return np.exp(-((x - mu)**2 / (2*var))+(1j*v*x/2))

L = 100
N = 3200
dx = 1/32
x = np.linspace(0, L, N)
func = lambda x: my_gaussian(x)*my_gaussian(x).conjugate()
C,e = quad(func, 0, L)

def task3(x):
    
    psi_0 = (C**-(1/2))*my_gaussian(x)

    H = (dx**-2)*diags([-1, 2, -1], [-1, 0, 1], shape=(N, N))
    H = sp.lil_matrix(H)
    H[0,0]=0.5*H[0,0]
    H[N-1,N-1]=0.5*H[N-1,N-1]
    lam, phi = eigsh(H, 400, which="SM")

    a = phi.T.dot(psi_0)
    psi = phi.dot(a*np.exp(-1j*lam*tmax))
    return psi*psi.conjugate()

func = lambda x: task3(x)*x
N1,e = quad(func, 0, L)

senderle · Accepted Answer · 2020-11-06 22:04:40Z

1

Initially this seems like a pretty straightforward shape mismatch. If you multiply a 2-d array by a 1-d array, numpy will only broadcast them automatically if the last dimension of the 2-d array is the same size as the 1-d array. Otherwise you have to explicitly reshape the 1-d array to be broadcastable with the 2-d array.

To fix the problem we need to know which arrays have the mismatched shapes. I added this to your code to see:

try:                                                                        
    psi = phi.dot(a * np.exp(-1j*lam*tmax).reshape(-1, 1))
except ValueError:
    print('phi.shape:', phi.shape)
    print('a.shape:', a.shape)
    print('lam.shape:', lam.shape)
    raise

The result was

phi.shape: (3200, 400)
a.shape: (400, 3200)
lam.shape: (400,)

So we need to reshape the lam term to be a column vector:

np.exp(-1j*lam*tmax).reshape(-1, 1))

This fixes the shape problem. But it doesn't fix the whole problem, because... the output of task3 is then a (3200,3200) array! This is, of course, not useful for an integration routine that expects a function that returns a single scalar.

This last problem is something you'll have to work out on your own, since I have no way of knowing what your goal is.

answered Nov 6, 2020 at 22:04

senderle

152k36 gold badges218 silver badges244 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Σ baryon Over a year ago

thanks for the answer, I have since found a better way to approach this problem :), I was supposed to get a single value after doing the integration. I think my logic was off

Collectives™ on Stack Overflow

Problem using integration scheme in numpy

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related