# NDSolve Step Size Effectively Zero

I'm trying to solve a PDE for a function of 2 variables. The most accurate parameterizations of this equation are very unwieldy and involve numerous piecewise elements, and so right now I'm trying to solve an approximate form. I'm getting step size effectively 0 errors, unless I increase the starting step size, in which case the solution I get is clearly unphysical. Since this is an astrophysics problem, I'm working over many orders of magnitude, which I suspect might be part of the issue.

Code:

z[x_]:= 2458.31 - 100.087 x + 1.23213 x^2 - 0.0046743 x^3

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

Ufit2[M_,t_]:=0.50519 + 3.127*10^10/M^2 - 274337./M + 2.12127*10^-10 M
-1.92858*10^-20 M^2 - 6.20762*10^-11 t;

fixedM[x_] :=
0.00006192808740866853/
591.5967816994163 (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)

DifEq = D[P[M, t],t] == -M l /sol^2 D[(1 - e)/e
lamavg[t] *3.93242 * Ufit2[M, t] P[M, t], M] //. {e ->
0.1,l -> 1.26 * 10^(31), sol -> 3*^8 };

soln = NDSolve[{DifEq, P[M, 8.15240949872944*^8] ==
1.629652976233948/M, P[1*^5 .263, t] == fixedM[t]}, P[M,
t], {t, 1000, 8.15240949872944*^8}, {M, 1*10^5.263, 10000000000},
Method -> { "StiffnessSwitching", "NonstiffTest" -> False}]


In terms of what I expect, this is a black hole mass function over the evolution of the universe, so P shouldn't fall off as time decreases, and there should be a difference of several orders of magnitude (especially at large t) between the low and high mass end.

Also I shortened the interval because the part here is the most important part, but ideally I would be able to continue the integration up to t of about 10 billion (not important at the moment though).

Thanks!

• "I'm working over many orders of magnitude" - is there really no choice of units that will give quantities of comparable size? – J. M. will be back soon Jul 23 '16 at 1:57
• What is the definition for fixedM – Young Jul 23 '16 at 2:00
• Do you have two accounts? You should not need the review of others to edit your own question. – bbgodfrey Jul 23 '16 at 2:17
• Can the time boundary be moved to one of the extremes of t – Young Jul 23 '16 at 2:23
• Please look here to have your accounts merged. – J. M. will be back soon Jul 23 '16 at 2:23

This evaluates without errors:

z[x_] := 2458.31 - 100.087 x + 1.23213 x^2 - 0.0046743 x^3

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

Ufit2[M_, t_] := 0.50519 + 3.127*10^10/M^2 - 274337/M + 2.12127*10^-10 M -
1.92858*10^-20 M^2 - 6.20762*10^-11 t;

e = 1/100;
l = (126/100)*10^(31);
sol = 3*^8;

DifEq = D[P[M, t],t] == -M l/ sol^2 D[(1 - e)/e lamavg[t]*3.93242*Ufit2[M, t] P[M, t], M];

fixedM[x_] := 0.00006192808740866853/ 591.5967816994163
(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)

soln = NDSolve[{
DifEq,
P[M, 8.15240949872944*^8] == 1.629652976233948/M,
P[0.263*^5, t] == fixedM[t]},
P[M, t], {M, 0.263*^5, 10000000000}, {t, 1000, 8.15240949872944*^8},
Method -> "StiffnessSwitching"]

Plot3D[P[M, t] /. soln, {M, 0.263*^5, 10000000000}, {t, 1000, 8.15240949872944*^8}]


• Maybe you meant to try Plot3D[P[M, t] /. soln,...] – Michael E2 Jul 23 '16 at 4:25
• @MichaelE2 ... oops ... your right – Young Jul 23 '16 at 4:27
• What about the rest of the interval? (I can get it to work with some trickiness, but your approach would be simpler.) – Michael E2 Jul 23 '16 at 4:47
• @MichaelE2 The OP said that I could set max t to equal the boundary condition of 8.15240949872944*^8 – Young Jul 23 '16 at 4:51
• This is about the right scale but not really the behavior I'm expecting...blackholes shouldn't just disappear! – basementDweller Jul 23 '16 at 5:03