Sistem PDE yang digabungkan dalam fisika atom

Aug 19 2020

Pertanyaan saya adalah tentang implementasi sistem PDE digabungkan ke rutinitas Mathematicas NDSolve. Saya sedang mempertimbangkan model mainan satu dimensi dalam fisika atom. Model menggambarkan dua bidang$\psi =\psi(t,z)$ dan $\sigma= \sigma(z;t)$ digabungkan satu sama lain yaitu $$ i \hbar \partial_t \psi = -\frac{\hbar^2 }{2 m} \psi_{zz} +V \psi +\frac{\hbar^2 \alpha_s }{m}\sigma^{-2} \left| \psi \right|^2 \psi+\frac{\hbar^2}{2m }\sigma^{-2}\psi+\frac{1}{2} m \omega_{\perp} \sigma^2 \psi +\frac{\hbar^2 }{2 m} \sigma^{-2}\sigma_z^2 \psi , \\ 0 =-\frac{\hbar^2}{4 m}\sigma \sigma_{zz}+\frac{\hbar^2 }{ m } \sigma^{-3} \sigma_z^2 -\frac{\hbar^2 }{4 m} \sigma \sigma_z \frac{1}{\left| \psi \right|^2} \left(\psi\psi_z^*+\psi^* \psi_z\right)+\frac{\hbar^2}{2 m }\sigma^{-3}-\frac{m \omega_{\perp}}{2} \sigma + 2 \frac{\hbar^2 \alpha_s}{m } \sigma^{-3} \left| \psi \right|^2 $$ Tambahan Saya memberlakukan persyaratan batas berkala untuk $\psi(-L/2,t) = \psi(L/2,t)$ dan $\sigma(-L/2,t) = \sigma(L/2,t)$ dan mengatur beberapa kondisi awal $\psi(z,0)=f(z)$ dan $\sigma(z,0)=g(z)$.

DIEDIT:

Ini adalah versi kode saya saat ini

    (*constants*)
h = 1; (* Planck constant *)
m = 1; (* particle mass *)
Subscript[\[Alpha], s] = 1; (* scattering length *)
\[Omega] = 1; (* frequency *)
V = 0; (* potential *)

(*ranges*)
L = 2; (*length of the box *)
tmin = 0;
tmax = 0.1;

(*equations*)
eqn1 = I  D[\[Psi][z, t], t] == -h^2/(2 m) D[\[Psi][z, t], z, z] + 
    V \[Psi][z, t] + 
    h^2 Subscript[\[Alpha], s]/
      m  \[Sigma][z, t]^(-2) Abs[\[Psi][z, t]]^2 \[Psi][z, t] + 
    h^2/(2 m) \[Sigma][z, t]^(-2) \[Psi][z, t] + 
    m \[Omega] /2 \[Sigma][z, t]^2 \[Psi][z, t] + 
    h^2/(2 m) \[Sigma][z, t]^(-2) D[\[Sigma][z, t], z]^2 \[Psi][z, t];

eqn2 = -h^2/(4 m) \[Sigma][z, t]  D[\[Sigma][z, t], z, z] ==  
   h^2/(2 m) \[Sigma][z, t]^(-3) D[\[Sigma][z, t], z]^2 -  
    h^2/(4 m) \[Sigma][z, t]   D[\[Sigma][z, t], z]  /
      Abs[\[Psi][z, t]]^2  ( \[Psi][z, t]  D[\[Psi][z, t], 
         z] + \[Psi][z, t] D[\[Psi][z, t], z]) + 
    h^2/(2 m) \[Sigma][z, t]^(-3)   - m \[Omega] /2 \[Sigma][z, t] + 
    2 h^2 Subscript[\[Alpha], s]/
      m \[Sigma][z, t]^(-3) Abs[\[Psi][z, t]]^2;

(*boundary conditions*)
bc = \[Psi][L/2, t] == \[Psi][-L/2, t];
bcwidth = \[Sigma][L/2, t] == \[Sigma][-L/2, t];

(*initial conditions*)
icwidth = \[Sigma][z, 0] == z^2 + 1;
icdwidth = D[\[Sigma][z, t], t] == 2 /. t -> 0;
icwave = \[Psi][z, 0] == Exp[-((z)^2)];

(*solve system*)
sol1 = NDSolve[{eqn1, eqn2, bc, bcwidth , icwave, icwidth, 
    icdwidth}, {\[Psi], \[Sigma]}, {z, -L/2, L/2}, {t, tmin, tmax}, 
   Compiled -> True, MaxSteps -> {500, Infinity}];

Sayangnya itu datang dengan dua masalah, yang pertama menyangkut Solver itu sendiri, karena tidak ada turunan waktu dalam persamaan saya untuk bidang kedua. $\sigma$ itu menangani sistem sebagai DAE dan memberikan dua Peringatan ini

NDSolve :: pdord: Beberapa fungsi memiliki orde diferensial nol, sehingga persamaan akan diselesaikan sebagai sistem persamaan diferensial-aljabar. >>

NDSolve :: mconly: Untuk metode IDA, hanya kode mesin yang tersedia. Tidak dapat melanjutkan dengan nilai kompleks atau di luar pengecualian floating-point. >>

Saya tidak tahu apakah ini adalah masalah "nyata" (saya menggunakan Mathematica 9.x). Yang kedua lebih bermasalah, ini menyangkut jumlah titik grid yang digunakan. Ini terutama berasal dari persamaan itu sendiri, saya kira dan menyebabkan kesalahan bahwa dia tidak dapat menemukan solusi yang tepat dalam batas toleransi.

NDSolve :: mxsst: Menggunakan jumlah titik kisi maksimum 500 yang diizinkan oleh opsi MaxPoints atau MinStepSize untuk variabel independen z. >>

NDSolve :: icfail: Tidak dapat menemukan kondisi awal yang memenuhi fungsi residual dalam toleransi yang ditentukan. Coba berikan kondisi awal untuk nilai dan turunan fungsi. >>

Saya juga mencoba memberinya data awal tambahan seperti yang disarankan oleh pesan kesalahan tetapi tidak berhasil. Pertanyaan Hal yang saya tidak tahu adalah apakah ada potensi untuk meningkatkan kode saya, atau jika peningkatan ke versi yang lebih baru dari Mathematica akan menyelesaikan masalah atau dalam kasus terburuk ini adalah sistem yang "terlalu jelek" untuk perawatan numerik.

Jawaban

1 AlexTrounev Aug 26 2020 at 21:54

Untuk mengatasi masalah seperti ini kita dapat membagi fungsi gelombang menjadi dua bagian $\psi=\psi_1+i\psi_2$. Juga kami menggunakan beberapa opsi untuk NDSolvemembuat masalah ini dapat dipecahkan. Anggap saja begitu$\sigma$ itu nyata, maka kita punya