0
$\begingroup$

I'm new to Mathematica and I'm trying to convert the output of a NDSolve problem into an exponential one. In particular I obtain the W(x) function ,but in the end I want to plot exp[W(x)]=Y(x). I naively tried the following with no results though

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

Mx = 100;
gx = 2;
gs = 106.75;
g = 100;
Mp = 10^19;
Y0 = 10^-20;
W0 = Log[Y0];
\[Sigma]v = 10^-6;

Clear[neq, f];
neq[x_Real, m_] := 
  NIntegrate[(gx/(2 Pi Pi)) Exp[-x Sqrt[1 + p p/(m m)]] p p, {p, 0, 
    Infinity}, MaxRecursion -> 50];
stot[x_] := (2 Pi Pi/45) gs Mx^3 /(x^3);
Yeq[x_, m_] := neq[x, m]/stot[x];
Weq[x_, m_] := Log[Yeq[x, m]];
RHS[x_, m_, \[Sigma]v_] := (0.264 Mp Mx gs \[Sigma]v/(Sqrt[g])) (Exp[
       2 Weq[x, m] - W[x]] - Exp[W[x]])/(x^2);
xmin = 10^-2; xmax = 10^2;
xdom = {x, xmin, xmax};
ODE = {W'[x] == RHS[x, Mx, \[Sigma]v]};
BC = {W[xmin] == W0};
Wsol = NDSolve[{W'[x] == RHS[x, Mx, \[Sigma]v], BC}, W, xdom, 
  Method -> "StiffnessSwitching"]

LogLogPlot[Evaluate[Y[x] /. Log[Wsol]], {x, 10^-2, 10^2}, 
 PlotRange -> Automatic]
$\endgroup$
3
  • $\begingroup$ Could you provide some background? Your equation seems a bit weird. Also you have constants that span 30 orders of magnitude, you should try to reparameterize the equations in order to have smaller constants. $\endgroup$
    – mattiav27
    May 2, 2022 at 11:54
  • $\begingroup$ Hi, I just added another question because your comment made me notice I wasn't looking at the right solution. So thanks! And also, if you want to take a look to the "general" problem of mine I'd be great ahah. Thanks in any case! $\endgroup$
    – Fredrigo6
    May 2, 2022 at 13:19
  • 1
    $\begingroup$ There must be something wrong with the magnitude of the parameter. You get Derivative[1][W][1/100] == 5.27765*10^34 as intial condition, which is utopically high. $\endgroup$
    – Akku14
    Jun 1, 2022 at 17:45

1 Answer 1

1
$\begingroup$

Leaving your specific problem aside for a moment, I just wanted to show you that you can obtain a function of your solution directly using NDSolveValue, with syntax practically identical to NDSolve.

For instance, the example below solves a simple differential equation modeling exponential decay, and then reports the logarithm of that solution as output.

sol = 
  NDSolveValue[
    {a'[t] == -1/2 a[t], a[0] == 1},
    Log[a[#]] &,
    {t, 0, 20}
  ]

Solution (Log of interpolating function) expressed as a pure function

You can use that solution as a pure function in further processing, e.g. for plotting:

Plot[sol[t], {t, 0, 10}]

plot of the solution from above, showing a linear decay with time

$\endgroup$

Your Answer

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

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