4
$\begingroup$

I have eqs about the NDSolve, I know this code given the solving automatically.

How can I find out what method is used behind the scenes? How can I gauge the reliability level, find how many iterations have been used, the order of method. How can I estimate the error?

I found hints on this site, but I still do not fully understand.

It is impossible to say NDSolve has automatically solution for publishing paper?

I used this code related to my system:

r = 0.431201; β = 2.99 *10^-6; σ = 0.7; δ = 0.57;
{m = 0.3, η = 0.1, μ = 0.1, ρ = 0.3};


S = {N1'[t] == r N1[t] (1 - β N1[t]) - η  N1[t] I1[t],
     I1'[t] == σ + (ρ  N1[t]  I1[t])/( m + N1[t]) - δ I1[t] - μ  N1[t] I1[t]};

c = {N1[0] == 1, I1[0] == 1.22};

Select[Flatten[
  Trace[
    NDSolve[{S, c}, {N1, I1}, {t, 0, 30}], 
    TraceInternal -> True]], 
  !FreeQ[#, Method | NDSolve`MethodData] &]

but I don't understand the output.

$\endgroup$
  • 2
    $\begingroup$ Partial duplicate: mathematica.stackexchange.com/questions/145/… $\endgroup$ – Michael E2 Mar 24 at 1:17
  • 1
    $\begingroup$ Another partial duplicate: mathematica.stackexchange.com/questions/102704/… $\endgroup$ – Michael E2 Mar 24 at 1:37
  • 1
    $\begingroup$ You say you don't understand some technique or other, nor the output of your Trace[] command. But the first is a very general statement about things already explained and the second is about a command that no one else can reproduce $\endgroup$ – Michael E2 Mar 24 at 1:44
  • 2
    $\begingroup$ "It is impossible to say NDSolve has automatically solution for publishing paper. " Simply saying "I've used NDSolve function of software Mathematica" is enough in many cases, AFAIK. $\endgroup$ – xzczd Mar 24 at 3:39
  • 3
    $\begingroup$ Well, if the reviewer insists on such stuff, given that your system isn't that difficult, a possible workaround at this point is to choose a primary method like classical RK4 to solve the problem. The way to choose classical RK4 in NDSolve can be found in tutorial/NDSolveExplicitRungeKutta#1456351317, then you just need to set Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, "Coefficients" -> ClassicalRungeKuttaCoefficients}, StartingStepSize -> 1/20000, MaxSteps -> Infinity in NDSolve. The solving process is slower but gives the same result as given by default. $\endgroup$ – xzczd Mar 24 at 3:59
6
$\begingroup$

Comment

In response to your question, you already got very valuable comments. I will just try to comment on

How can I estimate the error?

For this I am going to plot residual error at steps and time, which will show the reliability and accuracy of NDSolve,

r = 0.431201; \[Beta] = 2.99*10^-6; \[Sigma] = 0.7; \[Delta] = 0.57;
m = 0.3; \[Eta] = 0.1; \[Mu] = 0.1; \[Rho] = 0.3;   

ode = {N1'[t] == r N1[t] (1 - \[Beta] N1[t]) - \[Eta] N1[t] I1[t], 
   I1'[t] == \[Sigma] + (\[Rho] N1[t] I1[t])/(m + N1[t]) - \[Delta] I1[t] - \[Mu] N1[t] I1[t]};

bcs = {N1[0] == 1, I1[0] == 1.22};

residuals = ode /. Equal -> Subtract;

{s} = NDSolve[{ode, bcs}, {N1, I1}, {t, 20}, InterpolationOrder -> All];

N1["Coordinates"] /. s;

residuals /. t -> N1["Coordinates"] /. s;

ListPlot[Abs[Flatten /@ (residuals /. t -> N1["Coordinates"] /. s)], Frame -> True]

enter image description here

With[{data = {Table[{t, Abs@residuals[[1]]} /. s, {t, N1["Coordinates"] /. s // Flatten}]}}, 
 ListLogPlot[data, Frame -> True, PlotRange -> All]]

enter image description here

Note: I adopted the above from this website but unable to find the link.

$\endgroup$
  • $\begingroup$ Thank you so much @zhk but how can defend the axes for both figures? the first as x represented steps and y residual error. the second one x represents the t time and y residual error. sorry if my question is trivial but it is first time to see the code $\endgroup$ – sana alharbi Mar 24 at 9:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.