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

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

Maple 2015

Hi,

I have an ODE system with a free parameter P, that I can only solve numerically.
To simplify let's suppose this parameter has only 2 values P1 and P2.
Let SOL1 the solution of this system when P=P1 and SOL2 its solution for P=P2:

Again for a sake of simplification let's suppose the independent variable is t and the dependent one is x(t).

How can I plot the difference x1(t)-x2(t) where x1(t) is the solution corresponding to SOL1 and x2(t) the solution corresponding to SOL2
(sorry I screwed something up and changed the config of my keyboard making the interrogation point unaccessible)

Here is a detailed example

 > 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)
 (1)
 > # 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);
 > # 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)
 (2)
 > # 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
 >