> |
# 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) |
SOLVER = LSODE[ADAMSFULL]
> |
sol := dsolve(MC, numeric, method=lsode[adamsfull]); display( odeplot(sol, [t, X[1](t)], 0..9, color=red, gridlines=true, labels=["", "t"], title="X (red), V (blue)"), odeplot(sol, [t, V[1](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[1](t), X[1](t)], 0..9): ListOfPoints := op(1, op(1, aux)); N := LinearAlgebra:-Dimensions(ListOfPoints)[1];
|
|
(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[1](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[1](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[1](t) is the initial value of X[1](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[1](t)], 0..9, color=red, gridlines=true, labels=["", "t"], title="X (red), V (blue)"), odeplot(sol, [t, V[1](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[1](t), X[1](t)], 0..9): ListOfPoints := op(1, op(1, aux)): N := LinearAlgebra:-Dimensions(ListOfPoints)[1]: # 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 # #
|
|