ระบบ PDE ควบคู่ในฟิสิกส์อะตอม
คำถามของฉันเกี่ยวกับการนำระบบ PDE มาใช้กับรูทีน Mathematicas NDSolve ฉันกำลังพิจารณาแบบจำลองของเล่นมิติเดียวในฟิสิกส์อะตอม แบบจำลองอธิบายสองฟิลด์$\psi =\psi(t,z)$ และ $\sigma= \sigma(z;t)$ คู่กันคือ $$ 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}];
น่าเสียดายที่มันมาพร้อมกับสองปัญหาปัญหาแรกเกี่ยวข้องกับตัวแก้เองเนื่องจากไม่มีอนุพันธ์ของเวลาในสมการของฉันสำหรับฟิลด์ที่สอง $\sigma$ มันจัดการระบบเป็น DAE และให้คำเตือนสองข้อนี้
NDSolve :: pdord: ฟังก์ชันบางอย่างมีลำดับความแตกต่างเป็นศูนย์ดังนั้นสมการจะได้รับการแก้ไขเป็นระบบสมการเชิงอนุพันธ์ - พีชคณิต >>
NDSolve :: mconly: สำหรับวิธี IDA จะมีเฉพาะรหัสจริงของเครื่องเท่านั้น ไม่สามารถดำเนินการต่อด้วยค่าที่ซับซ้อนหรือเกินกว่าข้อยกเว้นทศนิยม >>
ฉันไม่รู้ว่านี่เป็นปัญหา "จริง" (ฉันใช้ Mathematica 9.x) อย่างที่สองเป็นปัญหามากกว่าซึ่งเกี่ยวข้องกับจำนวนจุดกริดที่ใช้ สิ่งนี้ส่วนใหญ่มาจากสมการเองที่ฉันเดาและทำให้เกิดข้อผิดพลาดที่เขาไม่สามารถหาทางออกที่เหมาะสมภายในขอบเขตความอดทนได้
NDSolve :: mxsst: การใช้จำนวนจุดกริดสูงสุด 500 ที่อนุญาตโดยตัวเลือก MaxPoints หรือ MinStepSize สำหรับตัวแปรอิสระ z >>
NDSolve :: icfail: ไม่พบเงื่อนไขเริ่มต้นที่ตรงตามฟังก์ชันที่เหลือภายในค่าความคลาดเคลื่อนที่ระบุ ลองกำหนดเงื่อนไขเริ่มต้นสำหรับทั้งค่าและอนุพันธ์ของฟังก์ชัน >>
ฉันยังพยายามให้ข้อมูลเบื้องต้นเพิ่มเติมแก่เขาตามที่แนะนำโดยข้อความแสดงข้อผิดพลาด แต่ไม่ประสบความสำเร็จ คำถามสิ่งที่ฉันไม่รู้คือมีศักยภาพในการปรับปรุงโค้ดของฉันหรือไม่หรือการอัปเกรด Mathematica เป็นเวอร์ชันที่ใหม่กว่าจะช่วยแก้ปัญหาได้หรือในกรณีที่เลวร้ายที่สุดระบบ "น่าเกลียดเกินไป" สำหรับการคำนวณเชิงตัวเลข
คำตอบ
ในการแก้ปัญหาประเภทนี้เราสามารถแบ่งฟังก์ชันคลื่นออกเป็นสองส่วน $\psi=\psi_1+i\psi_2$. นอกจากนี้เรายังใช้ตัวเลือกบางอย่างNDSolve
เพื่อทำให้ปัญหานี้สามารถแก้ไขได้ สมมติว่า$\sigma$ เป็นเรื่องจริงแล้วเราก็มี
(*constants*)h = 1;(*Planck constant*)m = 1;(*particle mass*)
Subscript[\[Alpha],
s] = 1;(*scattering length*)\[Omega] = 1;(*radial frequency*)V = \
0;(*longitudinal potential*)(*ranges*)L = 2;(*length of the box*)tmin \
= 0;
tmax = 0.1;
(*equations*)
eqn1 = { D[\[Psi]1[z, t], t] == -h^2/(2 m) D[\[Psi]2[z, t], z, z] +
V \[Psi]2[z, t] +
h^2 Subscript[\[Alpha], s]/
m \[Sigma][z,
t]^(-2) (\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2) \[Psi]2[z, t] +
h^2/(2 m) \[Sigma][z, t]^(-2) \[Psi]2[z, t] +
m \[Omega]/2 \[Sigma][z, t]^2 \[Psi]2[z, t] +
h^2/(2 m) \[Sigma][z, t]^(-2) D[\[Sigma][z, t], z]^2 \[Psi]2[z,
t], - D[\[Psi]2[z, t],
t] == -h^2/(2 m) D[\[Psi]1[z, t], z, z] + V \[Psi]1[z, t] +
h^2 Subscript[\[Alpha], s]/
m \[Sigma][z,
t]^(-2) (\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2) \[Psi]1[z, t] +
h^2/(2 m) \[Sigma][z, t]^(-2) \[Psi]1[z, t] +
m \[Omega]/2 \[Sigma][z, t]^2 \[Psi]1[z, t] +
h^2/(2 m) \[Sigma][z, t]^(-2) D[\[Sigma][z, t], z]^2 \[Psi]1[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]/(\[Psi]1[z, t]^2 + \[Psi]2[z,
t]^2) (D[(\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2), 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) (\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2);
(*boundary conditions*)
bc = {\[Psi]1[L/2, t] == \[Psi]1[-L/2, t], \[Psi]2[L/2,
t] == \[Psi]2[-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]1[z, 0] == Exp[-((z)^2)], \[Psi]2[z, 0] == 0};
(*solve system*)
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{Psi1, Psi2, S} =
NDSolveValue[{eqn1, eqn2, bc, bcwidth, icwave,
icwidth}, {\[Psi]1, \[Psi]2, \[Sigma]}, {z, -L/2, L/2}, {t,
tmin, tmax},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 81, "MaxPoints" -> 81,
"DifferenceOrder" -> "Pseudospectral"}}},
EvaluationMonitor :> (currentTime = t;)];];
การแสดงผลการแก้ปัญหาเชิงตัวเลข
{Plot3D[Psi1[z, t], {z, -L/2, L/2}, {t, tmin, tmax}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic,
PlotLabel -> "Re\[Psi]"],
Plot3D[Psi2[z, t], {z, -L/2, L/2}, {t, tmin, tmax}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic,
PlotLabel -> "Im\[Psi]"],
Plot3D[S[z, t], {z, -L/2, L/2}, {t, tmin, tmax}, Mesh -> None,
ColorFunction -> Hue, AxesLabel -> Automatic,
PlotLabel -> "\[Sigma]", PlotRange -> All]}
