0
$\begingroup$

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!

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ I think you need to transform your differential equation (substitutions) in order to get closer to a numerically solvable model. $\endgroup$
    – Roman
    Jun 3, 2019 at 11:30
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ulrich Neumann
    Jun 3, 2019 at 11:36
  • $\begingroup$ 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. $\endgroup$
    – Roman
    Jun 3, 2019 at 13:13
  • $\begingroup$ @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... $\endgroup$
    – Free_ion
    Jun 4, 2019 at 14:49
  • $\begingroup$ @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. $\endgroup$
    – Roman
    Jun 4, 2019 at 14:58

1 Answer 1

Reset to default
1
$\begingroup$

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

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Free_ion
    Jun 4, 2019 at 12:34
  • $\begingroup$ The change from D[..] to Derivative[..] changes a lot . Look at pot[y] in both versions. Only the second one can be evaluated/plotted! $\endgroup$
    – Ulrich Neumann
    Jun 4, 2019 at 12:43
  • $\begingroup$ 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. $\endgroup$
    – Free_ion
    Jun 4, 2019 at 12:55
  • $\begingroup$ 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...) $\endgroup$
    – Ulrich Neumann
    Jun 4, 2019 at 13:06
  • $\begingroup$ 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... $\endgroup$
    – Free_ion
    Jun 4, 2019 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.