Question:Can we trust Maple?

Question:Can we trust Maple?

 > # CAN WE TRUST MAPLE ? # # I want to solve numerically systems of ordinary differential equations (ODE) # # File "MC.m" below contains an example of one of them (a second degree ODE # reduced to a couple of 2 first order ODEs ; details about these equations # are of no importane here). # # I am going to show you a very disturbing behaviour of MAPLE ... # restart: with(plots): read "/Users/marcsancandi/Desktop/MAPLE++SCILAB/BUG-LSODE/MC.m": # A very quick look to the ODE system # # (0.2345... is the mass of the "moving mass" ; some of you will # probably recognize the structure of a mass-spring-damper system) # # Here again, the detailed expression of the 2 RHSs does not matter # map(u -> if is(rhs(u), numeric) then u else lhs(u)=RHS end if, convert(MC, list)); (1)

 > sol := dsolve(MC, numeric, method=lsode[adamsfull]); display(   odeplot(sol, [t, X(t)], 0..9, color=red, gridlines=true, labels=["", "t"], title="X (red), V (blue)"),   odeplot(sol, [t, V(t)], 0..9, color=blue) );  > # Now, for reasons detailled at the end of the worksheet, I do the # following "four stages program" # # # 1/ : capture de list of points odeplot returns aux          := plots:-odeplot(sol, [t, V(t), X(t)], 0..9): ListOfPoints := op(1, op(1, aux)); N            := LinearAlgebra:-Dimensions(ListOfPoints);  (2)
 > # 2/ : Extract from ListOfPoints the ones close to t=8 seconds AllTheTimes := convert(ListOfPoints[..,1], list): GoodRows    := zip((u,v)-> if verify(u, 7.8..8.2, `interval`) then v end if, AllTheTimes, [seq(1..N)]); (3)
 > # 3/ : print the extract of ListOfPoints for GoodRows only printf("\n What odeplot gives\n"); printf("------------------------------------------------\n"); printf("        t               V(t)           X(t)\n"); printf("------------------------------------------------\n"); for k in GoodRows do    printf("%-15.12f  %-15.12f  %-15.12f\n", seq(ListOfPoints[k,n], n=1..3)) end do; print():
 What odeplot gives ------------------------------------------------         t               V(t)           X(t) ------------------------------------------------ 7.830000000000   0.002817384181   0.006597587127 7.875000000000   0.002964569411   0.006727757784 7.920000000000   0.003166704318   0.006865745263 7.965000000000   0.003385619552   0.007013165777 8.010000000000   0.003558854571   0.007170093728 8.055000000000   0.003842118974   0.007336171867 8.100000000000   0.004302831211   0.007517231266 8.145000000000   0.003951552171   0.007708525735 8.190000000000   0.003970924524   0.007886798955 (4)
 > # 4/ Now evaluate "pointwise" sol(t) for the times retained print(); printf("\n What sol(t) gives for the same values of t\n"); printf("------------------------------------------------\n"); printf("        t               V(t)           X(t)\n"); printf("------------------------------------------------\n"); for MyTime in ListOfPoints[GoodRows, 1] do    MySol := map(u -> rhs(u), sol(MyTime)):    printf("%-15.12f  %-15.12f  %-15.12f\n", seq(MySol[n], n=1..3)) end do: print(): What sol(t) gives for the same values of t ------------------------------------------------         t               V(t)           X(t) ------------------------------------------------ 7.830000000000   0.000000000000   0.000987586956 7.875000000000   0.004000359093   0.001166940047 7.920000000000   0.004016468258   0.001347371362 7.965000000000   0.004032379719   0.001528593353 8.010000000000   0.004048109517   0.001710597919 8.055000000000   0.004063691706   0.001893215210 8.100000000000   0.004079090715   0.002076599938 8.145000000000   0.004094338959   0.002260632651 8.190000000000   0.004109918301   0.002445352658 (5)
 > # observe the differences between this table and the previous one !!!
 > # So ... what does the completely reconstructed "pointwise curve" look like ? # ... case of X alone PointwiseCurve := []: for MyTime in ListOfPoints[.., 1] do    MySol          := map(u -> rhs(u), sol(MyTime))[[1,3]]:    PointwiseCurve := [op(PointwiseCurve), MySol]; end do: display(   odeplot(sol, [t, X(t)], 0..9, color=red, gridlines=true),   PLOT(POINTS(PointwiseCurve)) ); > # Astonishing : the odeplot and the pointwise ona are exactly the same ! # # Let's make a zoom around t=8 display(   odeplot(sol, [t, X(t)], 0..9, color=red, gridlines=true, view=[7.5..8.5, 0.006..0.009]),   PLOT(POINTS(PointwiseCurve), VIEW(7.5..8.5, 0.006..0.009)) ); > # Do you see that X(8.1) is close to 0.0075 ? # Yes ? # Are you sure ? # Let's ask for a confirmation ... sol(8.1); (6)
 > # Conclusion : # # MAPLE is MAGIC !!! # # REMARK : Have you seen the answer for X(t) is the initial value of X(0) ? #          (see top of worksheet) # # # You don't believe me ... ABRACABADRA # ... and look to the values of sol(8.1) left to the stars # ListOfTimes := [1.0, 8.1, 10.0, 8.1, 1.0, 8.1, 10.0, 8.1]: printf("------------------------------------------------\n"); printf("   t          V(t)            X(t)\n"); printf("------------------------------------------------\n"); k := 0: for MyTime in ListOfTimes do    k     := k+1:    MySol := map(u -> rhs(u), sol(MyTime)):    if is(k, odd) then       printf("%6.3f  %-15.12f  %-15.12f\n", MyTime, seq(MySol[n], n=2..3))    else       printf("%6.3f  %-15.12f  %-15.12f  ***\n", MyTime, seq(MySol[n], n=2..3))    end if: end do:
 ------------------------------------------------    t          V(t)            X(t) ------------------------------------------------  1.000  0.000000000000   0.000987586956  8.100  0.004283834684   0.007512460823   *** 10.000  0.004613814863   0.015687722648  8.100  0.000000000000   0.000987586956   ***  1.000  0.000000000000   0.000987586956  8.100  0.004283834684   0.007512460823   *** 10.000  0.004613814863   0.015687722648  8.100  0.000000000000   0.000987586956   ***
 > # In fact the value of sol(8.1) depends on the value T had during # the previous evaluation sol(T). # Which suggests that, maybe, some global variable has not been # properly erased when sol(..) is evaluated (?)

SOLVER = RKF45   (slower)

 > # Let's do the same operations after replacing lsode by rkf45 # restart: with(plots): read "/Users/marcsancandi/Desktop/MAPLE++SCILAB/BUG-LSODE/MC.m": sol := dsolve(MC, numeric, method=rkf45, maxfun=500000); display(   odeplot(sol, [t, X(t)], 0..9, color=red, gridlines=true, labels=["", "t"], title="X (red), V (blue)"),   odeplot(sol, [t, V(t)], 0..9, color=blue) );  > # (Visual) comparison with the previous curves provides a very good agreement
 > # 1/ : capture de list of points odeplot returns aux          := plots:-odeplot(sol, [t, V(t), X(t)], 0..9): ListOfPoints := op(1, op(1, aux)): N            := LinearAlgebra:-Dimensions(ListOfPoints): # 2/ : Extract from LisrOfPoints the ones close to t=8 seconds AllTheTimes := convert(ListOfPoints[..,1], list): GoodRows    := zip((u,v)-> if verify(u, 7.8..8.2, `interval`) then v end if, AllTheTimes, [seq(1..N)]): # 3/ : print the extract of ListOfPoints for GoodRows only printf("\n What odeplot gives\n"); printf("------------------------------------------------\n"); printf("        t               V(t)           X(t)\n"); printf("------------------------------------------------\n"); for k in GoodRows do    printf("%-15.12f  %-15.12f  %-15.12f\n", seq(ListOfPoints[k,n], n=1..3)) end do; print():
 What odeplot gives ------------------------------------------------         t               V(t)           X(t) ------------------------------------------------ 7.830000000000   0.002815865032   0.006594922503 7.875000000000   0.002960849581   0.006724509859 7.920000000000   0.003162952386   0.006862303055 7.965000000000   0.003381243400   0.007009462432 8.010000000000   0.003556187082   0.007165994538 8.055000000000   0.003835074137   0.007331726709 8.100000000000   0.004281464768   0.007511890647 8.145000000000   0.003952037325   0.007703617600 8.190000000000   0.003971371977   0.007881895231 (7)
 > # Observe these results are very similar to those obtained with lsode # # But what are the results I obtaine by pointwise evaluations of # sol(t) for t in 7.8..8.2 ??? # # 4/ Now evaluate "pointwise" sol(t) for the times retained print(); printf("\n What sol(t) gives for the same values of t\n"); printf("------------------------------------------------\n"); printf("        t               V(t)           X(t)\n"); printf("------------------------------------------------\n"); for MyTime in ListOfPoints[GoodRows, 1] do    MySol := map(u -> rhs(u), sol(MyTime)):    printf("%-15.12f  %-15.12f  %-15.12f\n", seq(MySol[n], n=1..3)) end do: print(): What sol(t) gives for the same values of t ------------------------------------------------         t               V(t)           X(t) ------------------------------------------------ 7.830000000000   0.002815850376   0.006594922415 7.875000000000   0.002960860137   0.006724509763 7.920000000000   0.003162950433   0.006862302956 7.965000000000   0.003381213356   0.007009462330 8.010000000000   0.003556219238   0.007165994431 8.055000000000   0.003835077964   0.007331726592 8.100000000000   0.004281482695   0.007511890548 8.145000000000   0.003952038054   0.007703617628 8.190000000000   0.003971371542   0.007881895259 (8)
 > ListOfTimes := [1.0, 8.1, 10.0, 8.1, 1.0, 8.1, 10.0, 8.1]: printf("------------------------------------------------\n"); printf("   t          V(t)            X(t)\n"); printf("------------------------------------------------\n"); k := 0: for MyTime in ListOfTimes do    k     := k+1:    MySol := map(u -> rhs(u), sol(MyTime)):    if is(k, odd) then       printf("%6.3f  %-15.12f  %-15.12f\n", MyTime, seq(MySol[n], n=2..3))    else       printf("%6.3f  %-15.12f  %-15.12f  ***\n", MyTime, seq(MySol[n], n=2..3))    end if: end do:
 ------------------------------------------------    t          V(t)            X(t) ------------------------------------------------  1.000  0.001763825587   0.001540242179  8.100  0.004281482695   0.007511890548   *** 10.000  0.004614016677   0.015687242808  8.100  0.004281482695   0.007511890548   ***  1.000  0.001763825587   0.001540242179  8.100  0.004281482695   0.007511890548   *** 10.000  0.004614016677   0.015687242808  8.100  0.004281482695   0.007511890548   ***
 > # Wow, these two tables give very close results (at least if you don't look # farther than the tenth decimal position ... which could probably change # by fixing Digits to 15 or 20 (?) ... even if this hypothesis is not # fully satisfactory ?) # # Here pointwise evaluations of sol(t) always give the correct answer (up to # the decimal representation induced by the default value of Digits). # # #------------------------------------------------------------------------ # # Question : why the "magics" no longer operate with rkf45 ? #            (of course it is humor at the second degree !) # #------------------------------------------------------------------------ # # # It could be funny to play with "Magic MAPLE". # # But I have to solve serious problems, and among them, what interests me # in the solution of the "MC" system, is "events capturing". # Example of such an event "event" E_n is "the moving mass has reached the # position X = S_n where S_n is a given stroke (for exemple S_N = 0.006). # And "capturing E_n" then means "find the time T_n such that X(T_n) = S_n". # # I thus wrote a fixed-point method to do this, that is a procédure that # repeatedly asks for evaluations of procedure "sol". # But, for the behaviour depicted above, it doesn't work correctly with lsode. # So, despite the funny magical behaviour of MAPLE, I have a serious problem ! # # A problem I can sum up while saying # #     "I CANNOT TRUST MAPLE WHEN IT SAYS ME sol(t) = xxx" # # # Has anyone already been faced with these kind of behaviour ? # # Maplesoft support seems to beat around the bush with me, asking for the # ODEs equations, or criticizing the expressions of their RHSs ..., wanting # to see where these equations come from (they come from a more than 10.000 # lines application I have developped, but it's not the point here given # rkf45 proceeds correctly). # # If anyone is interested to investigate this problem, I can provide # her the "MC.m" file that contains the ODE system used here # # Thanks a lot for the time you spent to reading me # #
 > 