12

The problem is a combination of classic stiffness and weak singularities due to linear interpolation of the data. The stiffness needs a stiff solver plus a higher working precision, as @Ulrich has observed previously and we will explain below. For the weak singularites, NIntegrate has a method "InterpolationPointsSubdivision" to deal with ...


9

Stiffness and singularities Around 1960, things became completely different and everyone became aware that the world was full of stiff problems. (G. Dahlquist in Aiken 1985) Shampine's "flame" model $y'=y^2-y^3$, $y(0)=10^{-3}$, $0 \le x \le 2/y(0)$, integrated with a relative precision goal of $10^{-4}$. The red dots show the steps and ...


9

This is an old and rather often discussed issue. The output of CompiledFunctionTools`CompilePrint[myFun] contains MainEvaluate[ Hold[G1][ R0, R1]] and this tells us that G1 has not been inlined correctly. Best way to circumvent this issue is to use With for inlining: G1[x_, y_] := Block[{ks = x y, id = Sin[x y]}, {(id - ks), (id + ks)}] myFun = With[{f = ...


6

Clear["Global`*"] Use UnitConvert on all of the constants G = Quantity[1, "GravitationalConstant"] // UnitConvert; Ms = Quantity[1, "SolarMass"] // UnitConvert; Rs = Quantity[1, "SolarRadius"] // UnitConvert; P = Quantity[6, "Days"] // UnitConvert; Tdur = Quantity[4, "Hours"] // UnitConvert; You ...


6

I think that ODE systems with $\| X'\| \sim O(\| X \|^2)$ tend to be unstable, that is, a small rounding error has a chance to cause a solution to blow up. First, boundary-value problems (BVPs) are not guaranteed to have solutions, and without a proof or evidence that a solution exists, difficulty in solving one should raise the question whether there is a ...


6

What you call g[T] in your introduction seems to be equal to g[t]==-rate[t]/(1 + 2.2*10^-7 t^6/mass^2)^2 jacobian[t]. This values are of order 10^26! You need to increase accuracy inside NDSolve dramatically. Try F = NDSolveValue[{D[f[T], T] jacobian[T]^-1 == fDynamics[f[T], T], f[Tini] == 0}, f, {T, Tfinal, Tini}, WorkingPrecision -> 50] //Quiet ...


4

I will use the definition of g given in the answer to question you link to. g[t_, x_] := t^3 - t + x^2 f[t_?NumericQ] := NIntegrate[g[t, x], {x, -1, 1}] minPt = With[{min = FindMinimum[f[t], t ∈ Interval[{0, 1}]]}, {min[[2, 1, 2]], min[[1]]}] {0.57735, -0.103134} Plot[f[t], {t, 0, 1}, Epilog -> {Red, AbsolutePointSize[8], Point @ minPt}] This ...


4

The message General::munfl only occurs for machine numbers. So, you can fix the issue by using arbitrary precision numbers internally, and then converting to a machine number: function[t_, V_, l_] := N @ ReleaseHold[ SetPrecision[Hold[1/Sin[t]*LegendreP[-1/2+V*I,l+1,Cos[t]]], 30] ] Then: function[.00001, 10, 50] 1.34228*10^-189 + 0. I


1

You need to remove the curly brackets {} in {sf*1.0} etc., to yield: tone1 = Table[sig[t, sf*1.0], {t, 0, timerange, nyquist}]; tone2 = Table[sig[t, sf*1.2], {t, 0, timerange, nyquist}]; tone3 = Table[sig[t, sf*1.4], {t, 0, timerange, nyquist}]; tone4 = Table[sig[t, sf*1.6], {t, 0, timerange, nyquist}];


Only top voted, non community-wiki answers of a minimum length are eligible