差异的离散化.方程给出了奇怪的结果 [英] Discretization of a diff. equation gives weird results

问题描述

我有一个 Python 代码，可以为不同的势绘制波函数和能量:

I have a code in Python that draws wave functions and energy for different potentials:

    # -*- coding: cp1250 -*-
from math import *
from scipy.special import *
from pylab import *
from scipy.linalg import *

firebrick=(178./255.,34./255.,34./255.)
indianred=(176./255.,23./255.,31./255.)
steelblue=(70./255.,130./255.,180./255.)
slategray1=(198./255.,226./255.,255./255.)
slategray4=(108./255.,123./255.,139./255.)
lavender=(230./255.,230./255.,230./255.)
cobalt=(61./255.,89./255.,171./255.)
midnightblue=(25./255.,25./255.,112./255.)
forestgreen=(34./255.,139./255.,34./255.)

#definiranje mreze
Nmesh=512
L=4.0
dx=L/Nmesh
Xmax=L
x=arange(-L,L+0.0001,dx)
Npts=len(x)
numwav=0 #redni broj valne funkcije koji se iscrtava

V=zeros([Npts],float)
for i in range(Npts):
    V[i]=x[i]**50

a=zeros([2,Npts-2],float)
wave=zeros([Npts],float)

wave1=zeros([Npts],float)
encor=3.0/4*(3.0/4)**(1.0/3)

#numericko rjesenje
for i in range(1,Npts-1,1):
    a[0,i-1]= 1.0/dx**2+V[i]     #dijagonalni elementi
    a[1,i-1]=-1.0/dx**2/2        #elementi ispod dijagonale
a[1,Npts-3]=-99.0                #element se ne koristi
eig,vec=eig_banded(a,lower=1)    #rutina koja dijagonalizira tridijagonalnu matricu

for i in range(1,Npts-1,1):
    wave[i]=vec[i-1,numwav]
wave[0]=0.0             #valna funkcija u prvoj tocki na mrezi ima vrijednost nula
wave[Npts-1]=0.0        #valna funkcija u zadnjoj tocki na mrezi ima vrijednost nula

for i in range(1,Npts-1,1):
    wave1[i]=(2.0/pi*(3.0/4)**(1.0/3))**0.25*exp(-(3.0/4)**(1.0/3)*x[i]**2)
wave1[0]=0.0             #valna funkcija u prvoj tocki na mrezi ima vrijednost nula
#wave1[Npts-1]=0.0        #valna funkcija u zadnjoj tocki na mrezi ima vrijednost nula


#wave1=omjer*150*wave1+encor
wave=150*wave+eig[numwav]

#graf potencijala
line=plt.plot(x,V)
plt.setp(line,color='firebrick',linewidth=2)

#crtanje odabranog nivoa i odgovarajuce valne funkcije
plt.axhline(y=eig[numwav],linewidth=2,color='steelblue')
#plt.axhline(y=encor,linewidth=2,color='midnightblue')


#crtanje tocaka valne funkcije
plt.plot(x,wave,"b-",linewidth=2,color='forestgreen')
#plt.plot(x,wave1,"-",linewidth=2,color='indianred')

plt.xlabel(r'$x$',size=14)
plt.ylabel(r'$V(x)$',size=14)
plt.title(r'Valna funkcija i energija 3. pobuđenog stanja za $V(x)=x^{50}$')
plt.axis([-3.0,3.0,-8.0,100.0]) #raspon x i y osi
plt.grid(True)
plt.legend((r'$V(x)$',r'$E_0$',r'$\psi_0$'))
plt.show()

忽略被忽略的行，它们在这种情况下并不重要，语言 :D

Ignore the ignored lines, they are not important for this case, and the language :D

无论如何，我有问题.如果我绘制电位(V 部分)，假设达到 x^20，它绘制得很好，就像这样 x^6:

Anyhow, I have a problem. If I draw the potentials (the V part), for let's say up to x^20 it draws nice, like this for x^6:

如果我把潜力，说x^50它变成这样:

If I put the potential, say x^50 it becomes this:

那么问题出在哪里?他为什么会犯这么大的错误?它应该是平滑的，从理论上讲，当我到达 V(x)=x^p 对于非常大的 p (p → ∞) 点时，潜力应该去著名的无限方井，长这样:

So what seems to be the problem? Why is he making such big mistake? It should be smooth, and from the theory as I reach the point V(x)=x^p for very large p (p → ∞) the potential should go to the famous infinite square well, which looks like this:

所以我怀疑对于更大的潜力，我需要更多的点来在给定的范围内绘制它.那么我应该增加 Nmesh(网格)的数量吗?因为他说 Npts=len(x) - 他的分数.我对吗?这似乎合乎逻辑，但我想确定一下.

So I'm suspecting that for bigger potentials I need more points to draw it in the given range. So should I just increase the number of the Nmesh (grid)? Since he says that the Npts=len(x) - the number of points he's taking. Am I right? That seems logical, but I want to be certain.

感谢您的建议和帮助

我尝试扩展 Nmesh 但在非常大的数字时，我要么认为网格太大，要么存在内存问题.

I tried expanding the Nmesh but at very large numbers I either get that grid is too big, or that there is memory problems.

如果我拍摄，比如说 2048，我会得到相同的图片，但只是稍微移动了一点，更窄了.

If I take, say 2048 I get the same picture but just shifted a bit and narrower.

推荐答案

Use select argument for eig_banded():

Use select argument for eig_banded():

#!/usr/bin/env python
from __future__ import division
import functools
import math
import sys

from timeit import default_timer as timer

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as linalg

def report_time(func):
    @functools.wraps(func)
    def wrapper(*args, **kwargs):
        start = timer()
        try: return func(*args, **kwargs)
        finally:
            print '%s takes %.2f seconds' % (func.__name__, timer()-start)
    return wrapper

@report_time
def calc(Nmesh,POWER,L ,numwav=0):
    #
    dx=L/Nmesh
    x = np.arange(-L,L+0.0001,dx)
    Npts=len(x)

    V = x**POWER

    #
    ai = np.empty((2,Npts))   # ai[:,i] = a[:,i-1]
    ai[0,:] = 1/dx**2 + V     #
    ai[1,:] = -1.0/dx**2/2    #
    ai[1,Npts-2] = -99.0      #
    a = ai[:,1:-1]
    f = report_time(linalg.eig_banded)
    eig, vec = f(a, lower=True,overwrite_a_band=True,
                 select='i',select_range=(0,numwav)
                 ) #

    wave = np.empty(Npts)
    wave[1:-1] = vec[:,numwav]
    wave[0] = 0             #
    wave[Npts-1] = 0        #
    wave = 150*wave + eig[numwav]
    return x, V, wave, eig[numwav]

def main():
    try: numwav = int(sys.argv[1])
    except (IndexError, ValueError):
        numwav = 0
    POWER=100
    L=4.0
    Nmesh = 512
    print 'Nmesh=%d' % Nmesh
    x, V, wave, y = calc(Nmesh, POWER, L,numwav)

    #
    line = plt.plot(x,V)
    plt.setp(line,color='firebrick',linewidth=2)

    #
    plt.plot(x,wave,"b-",linewidth=2,color='forestgreen')

    #
    plt.axhline(y=y,linewidth=2,color='steelblue')

    plt.xlabel(r'$x$',size=14)
    plt.ylabel(r'$V(x)$',size=14)
    plt.title(r'$V(x)=x^{%d}$, ' % POWER)
    plt.axis([-(abs(L)-1), abs(L)-1,min(min(wave), y, min(V))-1, max(max(wave), y)+1]) #
    plt.grid(True)
    plt.legend((r'$V(x)$',r'$E_%d$' % numwav,r'$\psi_%d$' % numwav))
    plt.savefig('V_%d_%d_%d.png' % (Nmesh, POWER, numwav))
    plt.show()

if __name__=="__main__":
    main()

输出

Nmesh=512
eig_banded takes 0.01 seconds
calc takes 0.01 seconds

这是一个带有 numwav=4 的变体:

Here's a variant with numwav=4:

这篇关于差异的离散化.方程给出了奇怪的结果的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！