0

I am trying a system of two differential equations that looks pretty simple, but with a potential varying a lot on a short period centered on 0. The definition of the potential followed by the attempt at solving the differential equation is

a0[y_] := 10^(-22) y^2/k^2 + Sqrt[1 + y^2/k^2]; 
a2[y_] := D[D[a0[y], y], y]; 
f[y_] := 10^(-18)* a2[y]/a0[y]^3; 
g[y_] := 1 + f[y]; 
pot[y_] := D[D [g[y], y], y]/g[y]; (* Potential *)

yi = -300;
yf = 100;
k = 10^(-20);
sol = NDSolve[{x1'[y] == (pot[y] - 1) - x1[y]^2 + x2[y]^2, 
x2'[y] == -2*x1[y]*x2[y], x1[yi] == 0, x2[yi] == -1} , {x1,x2}, {y, yi, yf}, WorkingPrecision -> 30, MaxSteps -> Infinity][[1]]
Plot[Evaluate[x1[y] /. sol], {y, yi, yf}, PlotRange -> All, PlotPoints -> 500, AxesLabel -> {y, "x1(y)"}]
Plot[Evaluate[x2[y] /. sol], {y, yi, yf}, PlotRange -> All, PlotPoints -> 500, AxesLabel -> {y, "x2(y)"}]

And it appears that the solutions $x1$ and $x2$ are just equal to the initial conditions. I was expecting something growing after 0, since the potential there is about $10^{40}$, but I can't even see a little deviation from the initial values. Any comment, any clue on what I am doing wrong is greatly appreciated!

Free_ion
  • 58
  • 7
  • 1
    I think you need to transform your differential equation (substitutions) in order to get closer to a numerically solvable model. – Roman Jun 03 '19 at 11:30
  • 1
    If you change a2[y]:=... to a2[y_] := Derivative[2][a0][y]; and pot[y]:=... to pot[y_] := Derivative[2][g][y]/g[y] the potential can be evaluated and evaluates to Ò[10^-68]. That might explain the behavior. – Ulrich Neumann Jun 03 '19 at 11:36
  • What I mean is that your pot[y] is so singular at y=0 that you could simplify your differential equation dramatically by approximating pot[y_] = q*DiracDelta[y] (for some amplitude q) and trying to find an analytic solution of the differential equation. This analytic solution could then be used as a base to find the exact solution by appropriate transformations. – Roman Jun 03 '19 at 13:13
  • @Roman I tried what you suggested and, unsurprisingly, NDSolve encountered a step size effectively zero (for y around pi). I then tried to solve the system using DSolve with and without initial conditions, but the output MMA gave me was the input I gave. I am puzzled here... – Free_ion Jun 04 '19 at 14:49
  • @Free_ion for the Dirac $\delta$-suggestion you could solve the equations analytically, it may be easier than NDSolve. Are you familiar with the techniques for doing this in the presence of $\delta$-potentials? Essentially you solve the differential equation for $pot=0$ and then connect two different regions ($y<0$ and $y>0$) by a kink given by the $\delta$-function amplitude. – Roman Jun 04 '19 at 14:58
  • @Roman With pot=0, I got the following results for my functions: x1[y] -> (300 + y)/(90001 + 600 y + y^2) and 1/2 (-1 - Sqrt[(89999 + 600 y + y^2)^2/(90001 + 600 y + y^2)^2]) . The first problem here is that the plot is peaked at y=-300, when it should have been 0. Then, I tried to connect the two regions as you said using sol1[y_] := Piecewise[{{(300 + y)/(90001 + 600 y + y^2), yi <= y < 0 || 0 < y <= yf}, {10^(41)*DiracDelta[y], y == 0}}] but the potential doesn't seem to have any effect at all... – Free_ion Jun 13 '19 at 15:06
  • That's not what I meant by connecting them. You need to integrate the differential equations across the singularity to get rid of the Dirac-$\delta$. Maybe ask for more details at the math stackexchange, as I don't have the time to solve this right now. – Roman Jun 13 '19 at 18:36

1 Answers1

1

Just to illustrate my comment:

Clear[a2, pot]
a2[y_] := Derivative[2][a0][y];
pot[y_] := Derivative[2][g][y]/g[y];(*Potential*)


Plot[pot[y], {y, yi, yf}, PlotRange -> All]

enter image description here

and

Plot[pot[y], {y, yi, yf}, PlotRange -> All,PlotPoints -> {Automatic, {0}}]

enter image description here

Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55
  • Well, the potential itself is well-defined, and is actually similar to this one link, with different ranges in x and y. But changing the definition of the derivative does not give a different behaviour for x1 and x2, unfortunately. – Free_ion Jun 04 '19 at 12:34
  • The change from D[..] to Derivative[..] changes a lot . Look at pot[y] in both versions. Only the second one can be evaluated/plotted! – Ulrich Neumann Jun 04 '19 at 12:43
  • But is this really an effect of using Derivative[..] instead of D[..]? It seems to me the only difference is in the number of points used to make the plot. – Free_ion Jun 04 '19 at 12:55
  • Try to evaluate pot[0] and you'll get, expecting a number, a strange result \!\( \*SubscriptBox[\(\[PartialD]\), \(0\)]\( \*SubscriptBox[\(\[PartialD]\), \(0\)]\((1 + \*FractionBox[\( \*SubscriptBox[\(\[PartialD]\), \(0\)]\( \*SubscriptBox[\(\[PartialD]\), \(0\)]1\)\), \ \(1000000000000000000\)])\)\)\)/(1 + \!\( \*SubscriptBox[\(\[PartialD]\), \(0\)]\( \*SubscriptBox[\(\[PartialD]\), \(0\)]1\)\)/1000000000000000000) (sorry I don't know how to show the formated output...) – Ulrich Neumann Jun 04 '19 at 13:06
  • I guess I see your point. I always used to do pot[y] /. y -> 0 // N, with a number as output. But the output I get from NDSolve does not change whether I use one or the other definition for the derivative, i.e., I still get the same constant plots... – Free_ion Jun 04 '19 at 14:07