원자 물리학에서 결합 된 PDE 시스템

Aug 19 2020

내 질문은 Mathematicas NDSolve 루틴에 결합 된 PDE 시스템의 구현에 관한 것입니다. 저는 원자 물리학에서 1 차원 장난감 모델을 고려하고 있습니다. 이 모델은 두 가지 필드를 설명합니다.$\psi =\psi(t,z)$ 과 $\sigma= \sigma(z;t)$ 서로 결합 ie $$ 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 $$ 추가로 주기적 경계 조건을 부과합니다. $\psi(-L/2,t) = \psi(L/2,t)$ 과 $\sigma(-L/2,t) = \sigma(L/2,t)$ 몇 가지 초기 조건을 설정합니다. $\psi(z,0)=f(z)$ 과 $\sigma(z,0)=g(z)$.

편집 :

다음은 현재 코드 버전입니다.

    (*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}];

불행히도 두 가지 문제가 있습니다. 첫 번째 문제는 Solver 자체에 관한 것입니다. 두 번째 필드에 대한 방정식에는 시간 미분이 없기 때문입니다. $\sigma$ 시스템을 DAE로 처리하고이 두 가지 경고를 제공합니다.

NDSolve :: pdord : 일부 함수의 미분 차수가 0이므로 방정식이 미분 대수 방정식 시스템으로 풀립니다. >>

NDSolve :: mconly : IDA 메소드의 경우 기계 실제 코드 만 사용할 수 있습니다. 복잡한 값 또는 부동 소수점 예외 이상으로 계속할 수 없습니다. >>

이것이 "진짜"문제인지 모르겠습니다 (Mathematica 9.x를 사용하고 있습니다). 두 번째 문제는 사용 된 그리드 포인트의 양과 관련이 있습니다. 이것은 주로 내가 추측하는 방정식 자체에서 비롯되며 허용 오차 범위 내에서 적절한 솔루션을 찾을 수 없다는 오류를 유발합니다.

NDSolve :: mxsst : 독립 변수 z에 대해 MaxPoints 또는 MinStepSize 옵션에서 허용하는 최대 그리드 포인트 수 500을 사용합니다. >>

NDSolve :: icfail : 지정된 공차 내에서 잔류 함수를 만족하는 초기 조건을 찾을 수 없습니다. 함수의 값과 미분 모두에 대한 초기 조건을 지정해보십시오. >>

나는 또한 오류 메시지에서 제안한 추가 초기 데이터를 제공하려고했지만 성공하지 못했습니다. 질문 내가 모르는 것은 내 코드를 개선 할 가능성이 있는지, Mathematica의 최신 버전으로 업그레이드하면 문제가 해결되는지, 최악의 경우 수치 처리를위한 "너무 못생긴"시스템인지 여부입니다.

답변

1 AlexTrounev Aug 26 2020 at 21:54

이런 종류의 문제를 해결하기 위해 파동 함수를 두 부분으로 나눌 수 있습니다. $\psi=\psi_1+i\psi_2$. 또한 NDSolve이 문제를 해결할 수 있도록 몇 가지 옵션을 사용 합니다. 가정하자$\sigma$ 진짜라면 우리는