CODEWITHSUNDEEP
pythonnumpy
Σ baryon
Σ baryon

Reputation: 246

Problem using integration scheme in numpy

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)

enter image description here

Upvotes: 1

Views: 97

Answers (1)

senderle
senderle

Reputation: 150947

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.

Upvotes: 1

Related Questions