Numerical integration methods

[Albert Chan]

(07-21-2022 03:12 PM)Albert Chan Wrote:  Rewrite expression using expm1() does not help. It get "simplified" away.

It would be nice if user requested expm1(), expm1() is what you get.
As a patch, we can make our own version, that does not get simplify away.

(04-02-2014 05:18 PM)Thomas Klemm Wrote:  You could use \(e^x-1=2 \cdot sinh(\frac{x}{2}) \cdot e^{\frac{x}{2}}\) instead.

Code:

#cas
exp1(z) := 2*sinh(z/2)*exp(z/2);

log1(z) :=
BEGIN
  local y, z2;
  y := ln(1+z);
  z2 := exp1(y);
  RETURN y + (z-z2)/(z2+1); // taylor correction
END;
#end

log1(x) = log(1+x); exp1(x) = exp(x)-1
log1(x) is coded using result of exp1(), again to avoid it simplified away.

Bonus, this version work for complex numbers

Quote:Turns out the issue is less related to integration routine.
It is more related to inability to evaluate integrand accurately.
Without this fuzziness, Gaussian Quadrature should not have trouble.

Confirmed.

CAS> int(x*x*exp1(1/x^3), x, 1., 100.)                              → 4.80504346032
CAS> int(exp(3y)*exp1(exp(-3y)), y, 0., ln(100))                 → 4.80504346032
CAS> ln(100/1) - int((exp1(x)-x)/(x*x),x,1.,1/100.^3)/3      → 4.80504346032