1
$\begingroup$

I just solved a differential equation with NDSolve, and the behavior of the solution was not what I was expecting. In trying to track down where things went wrong, I plotted the solution function evaluated at the time that the initial condition is given, and it doesn't match what I had specified as the initial condition. I tried the solution proposed here with no luck.

Code:

a[t_] := Piecewise[{{e 0.000204499/Sqrt[47000] Sqrt[t], t < 46995}, {f t^3 + g t^2 + h t + i, 46995 <= t < 47005}, {c ((t)^0.667)/(6.39143*10^6) + d, 47005 <= t < 9800000000}, {q E^(Sqrt[1.989*10^(-20)/3] (t)) + b, t >= 9800000000}}] //. {q -> 0.292367892161671, b -> 0.10064726726594109, c -> 1.0386606919493655, d -> 0.00006529554347973175, e -> 1.3579559739239264, f -> -4.3633807877894195*10^(-13), g -> 6.152604549062771*10^(-8), h -> -0.002891832817142504, i -> 45.30731397090777, j -> 1.1404135497820182*10^-23, k -> 1.1236120846288418*10^-13, u -> -0.0054880391504018525, v -> 3.225813286569839*10^7}; zed[t_] := 1/a[t] - 1; zTable = {}; zTable = Table[{10^T, zed[10^T]}, {T, 3, 10, 0.0001}]; zedInt = Interpolation[zTable]; lamavg[t_] := Min[1, 0.01 + 0.07 zedInt[t]] eps = .1; l = (126/100)*10^(31); sol = 3*^8; moDot = 1.989*10^30 DifEq = D[P[M, t], t] == -31536000* M l/(moDot*sol^2) (1 - eps)/eps lamavg[t] D[P[M, t], M]; fixedM[x_] := 2.396946971556801*^-7*(32.59434080693661 - 2.0297457454952188*^-7 x + 1.2887918124478197*^-15 x^2 - 2.5861699533219344*^-25 x^3 + 1.941139394441828*^-35 x^4 - 5.134757798851362*^-46 x^5); fixedT[M_] := 26 (M^(-1)) E^(-4.7 M*10^(-10)); soln = NDSolve[{DifEq, P[M, 8.15240949872944*^8] == fixedT[M], P[10^5.263, t] == fixedM[t]}, P[M, t], {M, 10^5.263, 10000000000}, {t, 1000, 8.15240949872944*^8}, Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> "Projection"}}]

Plot of initial condition and problematic solution:

ϕ[m_, T_] := P[M, t] /. soln //. {M -> m, t -> T} Plot[{ϕ[M, 8.15240949872944*^8] /. soln, fixedT[M]}, {M, 10^6, 10^10}, PlotRange -> {0, 10^-7}]

Plot of initial condition and solution at that time

Thanks!

P.S. Also note that the interpolation for a and zedInt should be pretty good, a is continuous but not differentiable.

Edits: e-> eps, added a factor to the previously enormous coefficient

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ When I copy and past your NDSolve code in version 10.4, it returns with a message and unevaluated. $\endgroup$ Commented Jul 28, 2016 at 12:30
  • $\begingroup$ I just edited it to include the definition of lamavg[t], sorry about that $\endgroup$ Commented Jul 28, 2016 at 21:00

1 Answer 1

Reset to default
1
$\begingroup$

This is both an extended comment and an a qualitative answer.

First, because e is used for two different quantities, it is prudent to execute Clear[e] at the beginning of the code.

Second, the definition of lamavg is missing. From a related question, it probably is

lamavg[t_] := Min[1, 0.01 + 0.07 zedInt[t]]

Third, the NDSolve Method used seems to have no influence on the solution and can safely be omitted.

Fourth, the discrepancy between the initial condition and the value of the solution there does not seem so severe when plotted as follows.

LogLinearPlot[{ϕ[M, 8.15240949872944*^8] /. soln, fixedT[M]}, {M, 1.83 10^5, 10^10}, PlotRange -> All]

enter image description here

Because the scale length of the discrepancies is of order 0.0001 of the range of M, it seems likely that they arise from inadequate resolution. However, the appropriate resolution would be at least 10000 points, significantly increasing the running time of the code.

Fifth, because lamavg[t] is equal to 1 for t < 2.56*10^8, 

Quiet@LogLinearPlot[lamavg[t], {t, 10^3, 8.15240949872944*^8}]

enter image description here

the advective differential equation there is given by 

D[P[M, t], t] == -1.71882*10^24 M D[P[M, t], M]

which can be solved analytically to yield

f[1.71882*10^24 t - Log[M]]

where f is an arbitrary function chosen to fit the initial and boundary conditions. Further, we would expect that the solution should be essentially independent of M for t < 2.56*10^8, and that indeed is the case.

Plot3D[P[M, t] /. soln, {M, 1.83 10^5, 10^10}, {t, 1000, 8.15240949872944*^8}, AxesLabel -> {M, t, P}]

enter image description here

Finally, I would recommend that the variable M be replaced by a variable proportional to Log[M], which would accommodate better resolution of M by NDSolve, perhaps improving solution accuracy for very large t.

Addendum: More general analytical solution

The actual differential equation in the question

D[P[M, t], t] == -1.71882*10^24 lamavg[t] M D[P[M, t], M]

also can be solved analytically,

f[1.71882*10^24 s - Log[M]]

although the relationship between s and t, obtained by solving 

D[s[t], t] == lamavg[t]

is rather cumbersome. After the boundary conditions have been recast as function of s and Log[M], f can be determined everywhere.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks so much for the super-detailed answer! I have a couple questions that might be a bit newb-ish: 1. Why is it obvious that the solution isn't dependent on M at low t? 2. I checked the documentation for NDSolve options, but it's not clear to me how I would add more points/increase resolution. I can afford to let it run a bit longer at the moment. $\endgroup$ Commented Jul 28, 2016 at 22:58
  • $\begingroup$ Also I dunno if this is fair to ask, but the reason I wasn't expecting this behavior was that I expected the solution to only increase with time, which makes sense since my initial condition has a negative derivative everywhere, leading to a positive time derivative everywhere. This is clearly not what happens. I attributed this to the oscillations that appeared in the solution near the initial condition; is it possible that this behavior will change with higher resolution? $\endgroup$ Commented Jul 28, 2016 at 23:01
  • $\begingroup$ @basementDweller I did not expect the solution to depend strongly on M, because the coefficient of D[P[M, t], M] is so large, forcing the M-derivative to be exceedingly small to satisfy the PDE. The exception is very near t == 8.15240949872944*^8, where the initial condition forces a rapid variation in M, which must be balanced by a rapid variation in t. So, tiny steps in both M and t are needed to resolve the variations. Use MaxStepSize to set the maximum size of the steps in these two variables. $\endgroup$ Commented Jul 29, 2016 at 0:00
  • $\begingroup$ @basementDweller The oscillations at t == 8.15240949872944*^8 are a symptom, not a cause of the problem at large t. The problem is lack of resolution. Away from that point, the shape of the solution is determined by the shape of the boundary condition at M == 10^5.263. (Remember that information propagates in from the boundary for an advective equation, and the propagation velocity is very high here.) By the way, are you sure that the constants determining the coefficient in the PDE are correct? $\endgroup$ Commented Jul 29, 2016 at 0:10
  • $\begingroup$ the size of the coefficient bothered me for a really long time, and I finally figured out what was wrong. It is indeed missing a factor as written. I'll update it shortly. Changes the behavior significantly too. $\endgroup$ Commented Jul 29, 2016 at 0:24

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.