DSolve - incorrect solution to PDE

Question

Bug introduced in 12.2, persisting through 12.3.1.

I'm trying to solve a simple Schrödinger equation in external field:

DSolve[{-Derivative[2, 0][psi][x,t]/2 + F Sin[π t]x  psi[x,t]==
    I Derivative[0, 1][psi][x,t],psi[x,0]==Exp[-x^2]}, psi, {x,t}, 
    Assumptions-> t ∈ Reals && F ∈ Reals && x ∈ Reals]

And the solution that Mathematica spits out appears to be incorrect:

E^(I x (x/(-I + 2 t) - F t Sin[π t]))/Sqrt[1 + 2 I t]

When I plug it to the original equation and FullSimplify I get:

(E^(I x (x/(-I+2 t)-F t Sin[π t])) F t (2 π (I-2 t) x Cos[π t]+Sin[π t] 
    (-4 x+F t (-I+2 t) Sin[π t])))/(Sqrt[1+2 I t] (-I+2 t))==0

which obviously doesn't appear to be zero for all F values (although it is for F=0).

Here's a screenshot from my notebook:

I'm using MMA 12.3.0 on MacOS.

what's psi^(2,0) (etc.)? I think something might have gone wrong while copying... — thorimur, May 25 '21 at 20:22
ah, Derivative[2,0][psi] and Derivative[0,1][psi], probably, i see. — thorimur, May 25 '21 at 20:24
This is partial derivative alternative notion. Basically the same as D[psi[x,t],{x,2}], also psi^(0,1) would be D[psi[x,t],{t,1}] — Ranza, May 25 '21 at 20:25
right, it just doesn't parse as mathematica input syntax, so I couldn't copy-paste your code! I'll edit it so others can (and I checked to make sure it gives the same answer) — thorimur, May 25 '21 at 20:29
@thorimur thanks! BTW, Is there a way to easily copy compatible code from MMA? — Ranza, May 25 '21 at 20:34
@thorimur: The edited version of the code doesn't return a result on MMA 12.1.1 on Mac OS. Does it return a result in more recent versions? — Michael Seifert, May 25 '21 at 20:37
@MichaelSeifert ah yes, sorry, I'm on 12.3 and it works. if you can edit it to something more backwards-compatible though, assuming it is the syntax that's the issue, that'd be great :) if it's a difference in DSolve itself, though, that'd be interesting. — thorimur, May 25 '21 at 20:38
@thorimur: No, the syntax looks correct to me. I suspect that the PDE methods have been changed/upgraded in the most recent versions. — Michael Seifert, May 25 '21 at 20:41
@Ranza there are a couple of ways! you can get technically copy-pastable code by right clicking the cell bracket, and clicking Copy As -> Input Text. I say "technically" because this can be very ugly, involving RowBoxes. A cleaner output could be gotten by right clicking the cell bracket and first selecting Convert To > InputForm, then copying! — thorimur, May 25 '21 at 20:43
@thorimur thanks! This time I've did it with -> Plain Text, but next time I'll try -> Input Text. I remember that years ago there was a package/script for uploading to stack exchange, is this still available? — Ranza, May 25 '21 at 20:45
@MichaelSeifert - MMA 12.2 (Windows 10) returns the OP's result. — MelaGo, May 25 '21 at 21:20
fyi, you can write the ode like this ode = -1/2*D[psi[x, t], {x, 2}] + F*Sin[Pi *t]*x* psi[x, t] == I* D[psi[x, t], t] which is simpler than using Derivative notation. Btw, Maple 2021 can't solve it. It says there is no solution. — Nasser, May 25 '21 at 21:26
The analytic solution can be obtained by method of characteristics and Fourier transform, see, however I wanted to solve it in MMA... — Ranza, May 25 '21 at 21:31

xzczd · Accepted Answer · 2021-05-27T01:36:30.190

Since DSolve cannot solve the problem up to now i.e. v12.3.0, I'd like to add a (imperfect) work-around. The idea is similar to that in the linked paper. We first transform the problem to a initial value problem of 1st order PDE with the new-in-12.3 BilateralLaplaceTransform. (Notice bilateral Laplace transform is essentially a Fourier transform with special coefficient. )

With[{psi = psi[x, t]}, 
  eq = -D[psi, x, x]/2 + F Sin[π t] x psi == I D[psi, t];
  ic = psi == Exp[-x^2] /. t -> 0];
tsys = BilateralLaplaceTransform[{eq, ic}, x, s] /. 
   HoldPattern@BilateralLaplaceTransform[h_[x, t_], __] :> h[s, t] /. psi -> Ψ
(*
{(-(1/2)) s^2 Ψ[s, t] - F Sin[Pi t] Derivative[1, 0][Ψ][s, t] == 
  I Derivative[0, 1][Ψ][s, t], Ψ[s, 0] == E^(s^2/4) Sqrt[Pi]}
 *)

Then solve it with DSolve. DSolve spits out ifun warning, which isn't too surprising, because you've chosen the periodic F Sin[π t] as $E(t)$:

tsol = Ψ[s, t] /. DSolve[tsys, Ψ, {s, t}]
(*
{E^((π (I F + π s - I F Cos[π t])^2 - 
   I ArcSin[Cos[π t]] (F^2 - 2 (π s - I F Cos[π t])^2) + 
   1/2 I π (2 F^2 - 2 π^2 s^2 + 4 I F π s Cos[π t] + 
      F^2 Cos[2 π t]) + F (4 π s - 3 I F Cos[π t]) Sqrt[Sin[π t]^2])/(
  4 π^3)) Sqrt[π], 
 E^((π (I F + π s - I F Cos[π t])^2 + 
   I ArcSin[Cos[π t]] (F^2 - 2 (π s - I F Cos[π t])^2) - 
   1/2 I π (2 F^2 - 2 π^2 s^2 + 4 I F π s Cos[π t] + 
      F^2 Cos[2 π t]) + I F (4 I π s + 3 F Cos[π t]) Sqrt[Sin[π t]^2])/(
  4 π^3)) Sqrt[π]}
  *)

The tsol is probably valid only for certain interval of $t$ (maybe $0<t<1$, if I have to guess), but let's proceed anyway. The last step in principle is to transform back with

InverseBilateralLaplaceTransform[tsol, s, x]

But sadly, InverseBilateralLaplaceTransform cannot handle tsol. Given the paper doesn't include solution in the time domain either, I think this is acceptable.

"OK, but how do you know the method is correct? " A rigorous validation isn't easy, but experimental validation is. Since definition of inverse bilateral Laplace transform of $F(s)$ is $\frac{1}{2\pi\mathbb{i}} \int_{\gamma-\mathbb{i}\infty}^{\gamma+\mathbb{i}\infty}F(s)e^{st}ds$, we can validate for certain $(x,t)$ as follows:

integrand2[x_, t_] = 1/(2 Pi I) I tsol[[2]] Exp[I w x] /. s -> I w;
(validate the i.c.: )
Integrate[integrand2[x, 0], {w, -Infinity, Infinity}]
(* Exp[-x^2] *)
(* validate for x == 1, t == 1/2 )
eq /. psi -> integrand2 /. x -> 1 /. t -> 1/2 /. F -> 1 // Simplify;
Integrate[Subtract @@ %, {w, -Infinity, Infinity}]
( 0 *)

BTW, if a simpler $E(t)$ is chosen, the symbolic solution without integral can be found. For $E(t)=F t$:

(* Still, DSolve cannot correctly handle this: *)
With[{psi = psi[x, t]}, 
  eq = -D[psi, x, x]/2 + F t x psi == I D[psi, t];
  ic = psi == Exp[-x^2] /. t -> 0];
tsys = BilateralLaplaceTransform[{eq, ic}, x, s] /. 
   HoldPattern@BilateralLaplaceTransform[h_[x, t_], __] :> h[s, t] /. psi -> Ψ
tsol = Ψ[s, t] /. DSolve[tsys, Ψ, {s, t}]
{f1[x_, t_], f2[x_, t_]} = InverseBilateralLaplaceTransform[tsol, s, x]
(*
{1/(E^((8 I F^(5/2) t^6 + 9 (F t^2)^(5/2) + 180 F^(3/2) t^2 x + 240 I (F t^2)^(3/2) x - 
   360 I Sqrt[F] x^2)/(360 (-I Sqrt[F] + 2 Sqrt[F t^2]))) Sqrt[
  1 + (2 I Sqrt[F t^2])/Sqrt[F]]), E^((
 8 I F^(5/2) t^6 - 9 (F t^2)^(5/2) + 180 F^(3/2) t^2 x - 240 I (F t^2)^(3/2) x - 
  360 I Sqrt[F] x^2)/(360 (I Sqrt[F] + 2 Sqrt[F t^2])))/Sqrt[
 1 - (2 I Sqrt[F t^2])/Sqrt[F]]}
 *)

This isn't the end. Substituting the solution back to the PDE, only f1 turns out to be correct for $t>0$:

eq /. psi -> f1 // Simplify[#, t > 0]&
(* True *)
eq /. psi -> f2 // Simplify[#, t < 0]&
(* True *)

Nice solution, thank you! I'm accepting this as a workaround. — Ranza, May 26 '21 at 15:57

DSolve - incorrect solution to PDE

1 Answers1