2
$\begingroup$

I've been working on trying to get Mathematica to solve this system for a while now, and I have a couple of questions relating to various parts of the system:

  1. Interpolating functions: As a number of functions involved are complicated or piecewise or otherwise annoying, I make interpolating functions out of them so that Mathematica can just reference the values for later, like a look-up table, as opposed to calculating them again every time. Is this sound thinking/at all how these work?

  2. Non-numerical Errors: I get two errors at the end of this that don't make sense to me. One says that I'm asking for data outside the domain of the interpolating function (seems to be UfitInt here), which doesn't seem to be true. The other says there's a non-numerical derivative, which doesn't make sense to me given the NumericQs and defining everything. Is it related to InterpolatingFunctions?

  3. Consistency: Much less important, but after getting a "boundary conditions and initial conditions are not consistent" error earlier, I tweaked those a bit (fixedM and fixedT), and the error went away. I noticed however that at the point of intersection, they do not even satisfy this differential equation (by a significant amount, a couple powers of 10), even though Mathematica doesn't complain! (I noticed this to be the case before encountering the current trouble--I got some plots of a less precise version even with this being the case). Why would this be?

I've been staring at this problem for a while now, and I've noticed I've been starting to make some silly mistakes. It could be that that is the root of the issue here :)

Thanks!

Code: Define a foundational function

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];    

Define the Interpolating Function

Acub[t_] := (p1 + q1 z[t] + r1 z[t]^2 + s1 z[t]^3) //. {p1 -> .45, q1 -> -0.417, r1 -> 0.0403, s1 -> 0.0112}
Bcub[t_] := (p2 + q2 z[t] + r2 z[t]^2 + s2 z[t]^3) //. {p2 -> 0.861, q2 -> -0.674, r2 -> 0.16, s2 -> -0.0102}
Ccub[t_] := (p3 + q3 z[t] + r3 z[t]^2 + s3 z[t]^3) //. {p3 -> 0.637, q3 -> -0.603, r3 -> .17, s3 -> -0.0132}
U[M_, t_] := Min[Acub[t]/((M/M0)^Bcub[t] + (M/M0)^Ccub[t]), 1] //. {M0 -> 10^7.5}
uList = {}
For[i = 5, i < 10.235, i += 0.035, 
 For[j = 2.5, j < 9.26, j += 0.035, 
  AppendTo[uList, {10^i, 10^j, U[10^i, 10^j]}]]];
UfitInt[m_?NumericQ, t_?NumericQ] := Evaluate[Interpolation[uList]];

Other Functions:

lamavg[t_] := Min[1, 0.01 + 0.07 zedInt[t]];
e = 1/100;
l = (126/100)*10^(31);
sol = 3*^8;
secInYear = 31536000;

DifEq = D[P[M, t], t] == -secInYear* M l/sol^2 D[(1 - e)/e lamavg[t]*3.93242*UfitInt[M, t] 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));

Finally:

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}]

Errors:

  1. input value{183231, 1000} lies outside the range of data in the interpolating function. extrapolation will be used

  2. encountered non-numerical value for a derivative at M==183231

$\endgroup$
3
  • 1
    $\begingroup$ Approximating a discontinuous function with an InterpolatingFunction is not necessarily a good idea, as can be seen from Plot[{UnitStep[x], Interpolation[Table[{x, UnitStep[x]}, {x, -1, 1, .1}]][x]}, {x, -1, 1}]. Of course, using InterpolationOrder -> 0 eliminates this error at the discontinuity but is not very accurate for functions that are not piecewise constant. $\endgroup$
    – bbgodfrey
    Commented Jul 26, 2016 at 18:03
  • $\begingroup$ @bbgodfrey If you're referring to zed, if I've done everything correctly it should be continuous (if not differentiable) everywhere. Thanks though! $\endgroup$ Commented Jul 26, 2016 at 18:31
  • $\begingroup$ Seems to be a combination of simple mistakes. The mistake I can figure out is: 1. The definition of z[t] is missing. (Maybe it's zed?) 2. UfitInt is wrong, you've mixed up function and function relationship. Anyway, without the definition of z[t] this question is unanswerable. $\endgroup$
    – xzczd
    Commented Jan 4, 2017 at 11:35

0