restart:

with(plots):

# This is a 1D problem.
# A body with unit mass is submitted to an acceleration and to a solid-solid friction resistance
# The friction model is a simple Coulomb's

sys := {
         diff(x(t), t) = v(t),
         diff(v(t), t) = t - piecewise(v(t) > 0, f, v(t) < 0, -f, 0),
         x(0) = 0,
         v(0) = 0
       };

sol__0 := dsolve(sys, numeric, parameters=[f], maxfun=10^6)

# This system doesn't seem to be solved formally

sol__exact := dsolve(sys, {v(t), x(t)})

# Attempt of solving it numerically.
#
# Let's observe that the "v" solution must be an increasing function of t
# so the friction model could be replaced by piecewise(v(t) > 0, f, 0) for v(t) >=0.
# But this is a particular case (think of a sine acceleration for instance) and
# I keep using the initial Coulomb's model.
#
# The main problem, as it appears below is the discontinuity of the friction effect at v(t)=0
# (thus at the initial time).

sol := dsolve(sys, numeric, parameters=[f], maxfun=10^6):
sol(parameters=[1]):
sol(0.1);

Error, (in sol) cannot evaluate the solution further right of 0.15837836e-1, maxfun limit exceeded (see ?dsolve,maxfun for details)

# Clearly increasing maxfum will help to simulate on a longer time but it's not the
# most efficient way to handle this problem.
#
# The classical way is to use a regularized Coulomb's model ; for instance

sys := {
             diff(x(t), t) = v(t),
             diff(v(t), t) = t - f *tanh(v(t)/c),
             x(0) = 0,
             v(0) = 0
           };

sol__1 := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8)

# The smaller c the closer the solution to the one invoking the Coulomb's model.
#
# The question is: how can we measure as close to sol__exact we are?
# The idea is to choose a "small" value of c, lets say c__0, and to declare that the solution
# obtained with c__0 is "the exact solution".
#
# Studying the convergence of the solution as c__0 becomes smaller and smaller  will help
# choosing c__0.
# Of course the computational time is an increasing function of 1/c__0 : the smaller c__0
# the higher this time.

c__0 := 1e-8:  # arbitrary choice

sol__ref := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8):
sol__ref(parameters=[1, c__0]):

tmax := 10:

CodeTools:-Usage(sol__ref(tmax)):
numfun__0 := sol__ref(numfun);

odeplot(sol__ref, [t, v(t)], t=0..tmax)

memory used=33.20KiB, alloc change=0 bytes, cpu time=34.76s, real time=36.19s, gc time=0ns

# Assuming c__0 = 1e-8 is a value small enough to consider the solution plotted above
# is "the solution for the true Coulomb's model".
#
# Here are a few solutions for larger values of c:

cs   := 10^~[$-6..-1]:
tmax := 10:

printf("c = 1e-08  number of calls = %d\n", numfun__0):

for i from 1 to nops(cs) do
  sol__1(parameters=[1,cs[i]]):
  sol__1(tmax):
  couleur := ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]):
  Xplot__||i := odeplot(
                         sol__1, [t, x(t)], t=0..tmax,
                         legend=cs[i],
                         color=couleur
                       ):
  printf("c = %1.0e  number of calls = %d\n", cs[i], sol__1(numfun))
end do:


LY  := plottools:-transform((x, y) -> [x, log(y)]):
LXY := plottools:-transform((x, y) -> [log(x), log(y)]):

display(
         seq(Xplot__||i , i=1..nops(cs)),
         view=[0..3, 0..2], labels=[t, x], title="x(t) (lin-lin)"
       );

c = 1e-08  number of calls = 50831166

c = 1e-06  number of calls = 1308504

c = 1e-05  number of calls = 126848
c = 1e-04  number of calls = 12825
c = 1e-03  number of calls = 1513
c = 1e-02  number of calls = 456
c = 1e-01  number of calls = 246

# From this plot it seems that values of c less than 1e-4 give roughly the same curves.
#
# Remembering that the cimputational time is an increasing function of 1/c, we face the
# practical dilema of getting a solution "very close" to the ideal one (Coulomb's model)
# in a reasonable amount of time.
# Thus we need a small value of c nut not too small.
#
# The natural idea is to plot the difference between x(t) obtained with some value of c
# and x(t) obtained with our "reference value" c = c__0.
#
# How can we do this efficiently?
# (unfortunately piecewise and array/Array outputs are not allowed for parameterized solutions)
#
# For instance this is very slow


gap := proc(c, s1, s2, n)
  local s, h, k:
  sol__1(parameters=[1, c]):
  sol__1(s2):
  h := Matrix(n, 2):
  k := 1:
  for s from s1 to s2 by (s2-s1)/(n-1) do
    h[k, 1] := s:
    h[k, 2] := eval(x(t), sol__1(s)) - eval(x(t), sol__ref(s)):
    k       := k+1:
  end do:
  return h
end proc:


GAP := CodeTools:-Usage( gap(1e-5, 0, 1, 11) ):

Statistics:-ScatterPlot(GAP[..,1], GAP[..,2], symbol=solidcircle, color=blue)

memory used=399.35KiB, alloc change=0 bytes, cpu time=34.34s, real time=34.81s, gc time=0ns

restart:

with(Statistics):

# This is the naive coding of the problem
#
# It has the advantage to explain clearly (at least I hope so) what is to be done.

NaiveProc := proc(REP, a, b, N)
  local roll, rep, p, K, N0, N1:
  roll := rand(a .. b);
  N0   := 0:
  N1   := 0:
  for rep from 1 to REP do
    p := roll();
    K := Sample(Binomial(N, p), 1)[1]:
    if K=0 then
      N0 := N0+1:
    elif K=1 then
      N1 := N1+1:
    end if:
  end do:
  N0, N1:
end proc:
    
CodeTools:-Usage(NaiveProc(10^3, 0, 5e-4, 25000));

memory used=78.92MiB, alloc change=68.00MiB, cpu time=1.21s, real time=6.30s, gc time=24.80ms

# As the previous code is very slow I implemented several variants based on the fact
# that I care only on values of K equal to 0 and 1 and that it is then useless to sample
# a Binomial random variable to do this.
#
# Each of them is much faster than the naive coding... but can we even faster?
# For comparison, the same simulation realized in R takes less than 15s to run with REP=10^8

pi := (n, p, k) -> binomial(n, k) * (1-p)^(n-k) * p^k;


f1 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := Select(t -> is(t <= 0), U-P01):
  Got0  := Select(t -> is(t <= 0), U-P0 );
  N01   := numelems(Got01):
  N0    := numelems(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f1(10^5, 0, 5e-4, 25000));

print():

f2 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := select[flatten](`<=`, U-P01, 0):
  Got0  := select[flatten](`<=`, U-P0 , 0):
  N01   := numelems(Got01):
  N0    := numelems(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f2(10^5, 0, 5e-4, 25000));

print():

f3 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := Vector(REP, i -> `if`(U[i] <= P01[i], 1, 0)):
  Got0  := Vector(REP, i -> `if`(U[i] <= P0 [i], 1, 0)):
  N01   := add(Got01):
  N0    := add(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f3(10^5, 0, 5e-4, 25000));

memory used=376.83MiB, alloc change=296.77MiB, cpu time=4.12s, real time=4.03s, gc time=430.38ms

memory used=173.26MiB, alloc change=0 bytes, cpu time=2.30s, real time=1.98s, gc time=461.22ms

memory used=154.19MiB, alloc change=0 bytes, cpu time=2.28s, real time=2.01s, gc time=432.40ms

restart:

interface(version)

KL := (a, b) -> (1/4)*(2*ln(a+b)*a^2+4*ln(a+b)*b*a+2*ln(a+b)*b^2-2*ln(b)*b^2-a^2-2*a*b)/a

evalf(KL(1e-10, 1/2))

evalf(KL(1e-10, 0.5))

evalf(KL(1e-10, convert(0.5, rational)))

limit(KL(a, 1/2), a=0, left);
limit(KL(a, 1/2), a=0, right)

limit(KL(a, 0.5), a=0, left);
limit(KL(a, 0.5), a=0, right)