1
$\begingroup$

In Excel I compiled some spreadsheets implementing the Runge Kutta method of the 4th order. Initially I tested with some systems of differential equations that are known to the analytical solution and the error was always the order of 10^-7 - 10^-5.

So I passed to the resolution of some of differentiable system of equations can not be solved analytically and the graphs obtained were perfectly coincident with those Wolfram Mathematica 11.0.1. It seemed to work wonders when all I had to face in the following system of differential equations:

ode1 = D[(a[s]^2 + 4a[s]^4t[s]^2)t'[s], s] - 4a[s]^4t[s]t'[s]^2 - 
         (a[s] + 8a[s]^3t[s]^2)a'[s]^2 == 0;

ode2 = D[4a[s]^2t[s]^4a'[s], s] - 4a[s]t[s]^4a'[s]^2 - 8a[s]^2t[s]^3t'[s]^2 == 0;

ic = {a[-10] == 3.03721, a'[-10] == -0.0555245, t[-10] == 1.81013, t'[-10] == -0.0608313};

sol = NDSolve[{ode1, ode2, ic}, {a, t}, {s, -10, 10}];

Plot[Evaluate[{a[s], t[s]} /. sol],
     {s, -10, 10},
     AxesLabel -> {"s", "a, t"},
     PlotLegends -> {"a[s]", "t[s]"},
     PlotRange -> {{-10, 10}, {0, 5}}]

or equivalently:

sol = NDSolve[
   {a'[s] == b[s],
    b'[s] == (-b[s]^2 t[s] - 4 a[s] b[s] u[s] + 2 a[s] u[s]^2)/(a[s] t[s]),
    t'[s] == u[s],
    u'[s] == (b[s]^2 + 8 a[s]^2 b[s]^2 t[s]^2 - 2 b[s] u[s] - 16 a[s]^2 b[s] t[s]^2 u[s] - 
              4 a[s]^3 t[s] u[s]^2)/(a[s] (1 + 4 a[s]^2 t[s]^2)),
    a[-10] == 3.03721,
    b[-10] == -0.0555245,
    t[-10] == 1.81013,
    u[-10] == -0.0608313},
   {a, b, t, u}, {s, -10, 10}];

Plot[Evaluate[{a[s], t[s]} /. sol],
     {s, -10, 10},
     AxesLabel -> {"s", "a, t"},
     PlotLegends -> {"a[s]", "t[s]"},
     PlotRange -> {{-10, 10}, {0, 5}}]

enter image description here

while in Microsoft Excel 2016 I get the following graph:

enter image description here

You might point out some strategies to check the accuracy of such graphs? How do you be certain of the graphs obtained from MMA without the possibility of an analytical feedback? It may be that there are numeric problems?

Thanks so much!

$\endgroup$
  • 1
    $\begingroup$ I would try to change the method and various options, and see if the solution changes. One of the easiest things to try is to see if the solution is sensitive to the MaxStepSize setting. If possible, try to do the same in Excel (does it have a similar setting)? $\endgroup$ – Szabolcs Mar 29 '17 at 9:18
  • 1
    $\begingroup$ It is also good to double check that you are really solving the same equations in the two systems. Little mistakes are too easy to overlook with such a long equation. $\endgroup$ – Szabolcs Mar 29 '17 at 9:19
  • 2
    $\begingroup$ It might interest you to know that with NDSolve[], there's often no need to convert high-order DEs into a bunch of first-order DEs. So, NDSolve[{y''[x] == -y[x], y[0] == 1, y'[0] == 0}, y, {x, 0, Pi}] works directly, without the fuss of having to do the conversion to a pair of DEs. $\endgroup$ – J. M. will be back soon Mar 29 '17 at 9:21
  • $\begingroup$ @Szabolcs Of course I had already tried changing the integration step both in MMA and Excel, but does not change anything substantial. I'm aware that just a small mistake to unpack everything, but after a week I think I wrote everything correctly, so I think the mistake is to look just in numbers. $\endgroup$ – TeM Mar 29 '17 at 9:35
  • 1
    $\begingroup$ It might be helpful to attach an Excel sheet as well, so people can compare themselves. (Personally, I have no idea how to solve DEs in Excel.) $\endgroup$ – Szabolcs Mar 29 '17 at 9:37