Using Matlab to solve 1D Schrödinger Equation (Strange Eigenfunctions)

39 ビュー (過去 30 日間)
Niklas Kurz
Niklas Kurz 2023 年 5 月 31 日
コメント済み: Torsten 2023 年 6 月 2 日
Following this splendid tutorial (in Python) I tried recreating the code in Matlab
First of, here the Python Code:
% import numpy as np
%
% import matplotlib.pyplot as plt
%
% from scipy.linalg import eigh_tridiagonal
%
% N = 2000
%
% dy = 1/N
%
% y = np.linspace(0,1,N+1)
%
% def mL2V(y):
% return 1000*(y-1/2)**2
%
% d = 1/dy**2 + mL2V(y)[1:-1]
%
% e = -1/(2*dy**2) * np.ones(len(d)-1)
%
% w, v = eigh_tridiagonal(d,e)
%
% plt.plot(v.T[0])
Now my attempt converting this to Matlab:
%% Initialise
% Number of Steps
N = 100;
dy = 1/N;
y = linspace(0,1,N+1);
% Potential
V = @(y,m,L) m*L*(y-1/2).^2;
Vy = V(y,1,1);
% Differential Matrix
d = 1/dy.^2 + Vy(1:end-1);
n = length(d);
e = -1/(2*dy.^2).*ones(n,1);
A = spdiags([e d' e],[-1 0 1],n,n);
%% Solve
[v,w] = eig(full(A));
%% Plot
plot(y(1:end-1),v(:,2),'-','LineWidth',2)
axis([0 1 -1 1])
However the solutions from Matlab are totally off also comparing them with textbook-ones
How is this to explain?
___________
Background:
the goal is to solve the 1D SG:
whose solution is mainly determined by the potential V
setting we get:
writing the derivatives in discrete form brings:
This can be written in matrix-shape:
so basically a Matrix with a main diagonal consisting of:
and two side diagonals with a constants entry of:
The initial conditions are implemented already: respectively
Now the function Ψ is directly given by determing the eigenfunctions of the Matrix
_________________________________________________________________
- This is where the magic happens and where Matlab somehow finds solutions that don't suite the physics content.
  2 件のコメント
Torsten
Torsten 2023 年 5 月 31 日
Most probably because you didn't transfer the Python code correctly.
But since I don't have experience with Python, it would be easier if you include the mathematical description of your problem instead of code in a different computer language.
Niklas Kurz
Niklas Kurz 2023 年 6 月 1 日
I updated the physics. Me neither I'am not having many experiences in Python. I just thought it might help some of you.

サインインしてコメントする。

採用された回答

Torsten
Torsten 2023 年 6 月 1 日
編集済み: Torsten 2023 年 6 月 1 日
N = 2000;
dy = 1/N;
y = linspace(0,1,N+1).';
mL2V = @(y)1000*(y-0.5).^2;
e = -1/(2*dy^2)*ones(N-1,1);
d = 1/dy^2 + mL2V(y(2:end-1));
A = spdiags([e d e],-1:1,N-1,N-1);
A = full(A);
[v,w] = eig(A);
hold on
plot(v(:,1).^2)
plot(v(:,2).^2)
plot(v(:,3).^2)
hold off
grid on
  2 件のコメント
Niklas Kurz
Niklas Kurz 2023 年 6 月 2 日
Awesome, thank you! Apparently my matrix was slightly different but my potential clearly too weak, causing the shape to differ strongly. Now it's all fine and my faith in Matlab still strong.
Torsten
Torsten 2023 年 6 月 2 日
Especially the setting of the ylimits was not adequate.

サインインしてコメントする。

その他の回答 (1 件)

James Tursa
James Tursa 2023 年 6 月 1 日
編集済み: James Tursa 2023 年 6 月 1 日
Python:
return is misspelled retrun
The mL2V( ) function returns 1000*(y-1/2)**2
I don't see where y is defined prior to the line d = 1/dy**2 + mL2V(y)[1:-1]
MATLAB:
The V( ) function returns 1*1*(y-1/2).^2
I stopped looking after that. Maybe post a working Python code that we can run on the side before comparing that to MATLAB.
  1 件のコメント
Niklas Kurz
Niklas Kurz 2023 年 6 月 1 日
編集済み: Niklas Kurz 2023 年 6 月 1 日
I also updated the Python Code. Works like a charm now.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeCall Python from MATLAB についてさらに検索

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by