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:
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?
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?
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:
input value{183231, 1000} lies outside the range of data in the interpolating function. extrapolation will be used
encountered non-numerical value for a derivative at M==183231
InterpolatingFunction
is not necessarily a good idea, as can be seen fromPlot[{UnitStep[x], Interpolation[Table[{x, UnitStep[x]}, {x, -1, 1, .1}]][x]}, {x, -1, 1}]
. Of course, usingInterpolationOrder -> 0
eliminates this error at the discontinuity but is not very accurate for functions that are not piecewise constant. $\endgroup$z[t]
is missing. (Maybe it'szed
?) 2.UfitInt
is wrong, you've mixed up function and function relationship. Anyway, without the definition ofz[t]
this question is unanswerable. $\endgroup$