1

I’m using Mathematica for some conceptually simple algebra, that will ultimately land in Excel or a compiled language.

$Assumptions = Element[{r0,r1,t0,t1,disc0}, Reals] && (t1>t0) && (disc0>0)
(* Also, typically, r0 ≈ r1 && r0 ≈ 0 && r1 ≈ 0 *)

r[t_] = ((t1-t)r0+(t-t0)r1) / (t1-t0)   (* Linear interpolation *)


disc[t_] = disc0 Exp[-Integrate[r[tt], {tt, 0, t}]]


ans = Integrate[disc[t], {t,t0,t1}]


ans // FullSimplify

It works, in that it produces an answer. But the answer has multiple divisions by (r0 - r1), or, worse, by its Sqrt, As said in the first comment in the code, r0 ≈ r1 && r0 ≈ 0 && r1 ≈ 0. Hence dividing by (r0 - r1) is numerically unwise. Please, is there a means of instructing Mathematica of the approximate equalities, so as to discourage division by things near zero?

Perhaps relevant: $Version = “9.0 for Mac OS X x86 (64-bit) (January 24, 2013)”.

Thank you.

jdaw1
  • 499
  • 2
  • 9
  • Function definitions should use := rather than =. – A.G. Apr 23 '21 at 21:49
  • Not true, @A.G. and this poor advice leads to many problems for beginners. Both immediate and delayed assignments have important roles to play (tutorial), and here is a typical example where an immediate assignment makes more sense than a delayed one because we want the Integrate command to be executed at definition-time (i.e., only once), not at function-call time (i.e., at every call). – Roman Apr 24 '21 at 17:26
  • @Roman I only knew about https://reference.wolfram.com/language/tutorial/FunctionsAndPrograms.html ; do you have a reference for using = in defining functions? – A.G. Apr 24 '21 at 17:41
  • Here and here and here and links therein. Keep in mind that Mathematica has no concept of "function", only assignments and pattern-matching; there is no category distinction between a=5 and a[x_]=5. – Roman Apr 24 '21 at 17:58

2 Answers2

1

One way around the issue would be to separate three cases: r0==r1, r0<r1 and r0>r1:

$Assumptions = 
  Element[{r0, r1, t0, t1, disc0}, Reals] && (t1 > t0) && (disc0 > 0);
r[t_] := (r0*(t1-t) + r1*(t-t0))/(t1-t0);  (*Linear interpolation*)
disc[t_] := disc0 Exp[-Integrate[r[tt], {tt, 0, t}]];
ans[t_] := FullSimplify[Integrate[disc[t], {t, t0, t1}]];

Assuming[r0 == r1, ans[t]]
Assuming[r0  < r1, ans[t]]
Assuming[r0  > r1, ans[t]]

Output

A.G.
  • 4,362
  • 13
  • 18
  • Thank you. However, (Erf[…/Sqrt[r0−r1]] − Erf[…/Sqrt[r0−r1]]) … / Sqrt[r0−r1] is still a numerical disaster. – jdaw1 Apr 23 '21 at 22:27
0

A series expansion could help. You say that $r_1\approx r_0$, so let's set $r_1=r_0\cdot\rho$ with $\rho=\frac{r_1}{r_0}$ capturing the relationship between $r_0$ and $r_1$.

Now we can series-expand for small $r_0$ because we've tied the magnitudes of $r_1\approx r_0$ together through $\rho$:

Series[ans /. r1 -> r0*ρ, {r0, 0, 1}] // FullSimplify
(*    disc0 * (t1-t0) +
      disc0/6 * (2*t0*t1*(ρ-1) - t1^2*(ρ+2) + t0^2*(2ρ+1)) * r0 +
      O[r0]^(3/2)                                                    *)

(you can use a higher-oder expansion if you like).

This formula should be stable for small $r_0$. We can back-substitute $r_1$ into it:

ans0 = Normal[%] /. ρ -> r1/r0 // FullSimplify
(*    disc0/6*(t0*(r0*t0+2*r1*t0-6) + 2*(3-r0*t0+r1*t0)*t1 - (2*r0+r1)*t1^2)    *)

Let's try out this approximation ans0:

{ans, ans0} /. {disc0 -> 1, t0 -> 0.3, t1 -> 0.8, r0 -> 0.01, r1 -> 0.0101}
(*    {0.497258 + 0. I, 0.49725}    *)

A higher-order series-expansion may give a more accurate approximation, if needed.

Roman
  • 47,322
  • 2
  • 55
  • 121
  • Ooh, that is a clean way to capture the relationship. Thank you. I’ll explore this, then perhaps mark as answer. But because they are close to zero, I might expand about the difference rather than about the ratio. – jdaw1 Apr 25 '21 at 18:50