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$

closed as off-topic by xzczd, corey979, m_goldberg, Feyre, MarcoB Jan 4 '17 at 23:05

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – xzczd, corey979, m_goldberg, Feyre, MarcoB
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 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 Jul 26 '16 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$ – basementDweller Jul 26 '16 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 Jan 4 '17 at 11:35