Question: How to plot the difference between 2 dsolve[numeric] solutions

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