1
$\begingroup$

I computed the NUMERIC solution to a non-linear system of ODE. Kfn[t] = K(t), t>=0, is an entry in the solution. $K(t)$ is periodic with period inertialPeriod. I want to compute some maxima.

I coded this package:

$Assumptions = t >= 0
eq7 := Kf'[t] + Kf[t]^2 - L1f[t]^2 + f L1f[t] + 2 gp A1f[t] + 
   s Kf[t] == 0
eq8 := L1f'[t] + 2 Kf[t] L1f[t] - f Kf[t] + s L1f[t] == 0
eq9 := A0f'[t] + 2 Kf[t] A0f[t] == 0
eq10 := A1f'[t] + 4 Kf[t] A1f[t] == 0
sys = {eq7, eq8, eq9, eq10}
funs = {Kf[t], L1f[t], A0f[t], A1f[t]}
{gp = 1/100, omega = 2 Pi/(24 3600), fiPos = 19 Pi/90, 
 f = 2 omega Sin[fiPos], s = 10^(-6)}
inertialPeriod = 2 Pi/f
initVals = {(Pi*Sin[(19*Pi)/90])/
   1080000, -1491/1000000000 + (Pi*Sin[(19*Pi)/90])/43200, 500, 
  50*(2223081/1000000000000000000 - (13*Pi^2*Sin[(19*Pi)/90]^2)/
      24300000000)}
initCond = {Kf[0] == initVals[[1]], L1f[0] == initVals[[2]], 
  A0f[0] == initVals[[3]], A1f[0] == initVals[[4]]}
sysInit = Join[sys, initCond]
tFin = 16 inertialPeriod
solN = NDSolve[sysInit, funs, {t, 0, tFin}]
{Kfn[t_], L1fn[t_], A0fn[t_], 
  A1fn[t_]} = {Kf[t], L1f[t], A0f[t], A1f[t]} /. First[solN]
tip = inertialPeriod
t0 = 0
pltKtest = 
 Plot[{Kfn[t], 0}, {t, 0, tFin}, PlotLegends -> "Expressions", 
  FrameLabel -> {"seconds", "K"}, Frame -> True, 
  GridLines -> Automatic]
iFin = 16
(*compute iFin-1 maximum points,add t=0;duplicated points can appear*)
findTmaxVdup = Module[
  {i, t, ti, maxi},
  tMaxVdup = {0}; KMaxVdup = {0};
  For[i = 1, i <= iFin - 1, i++,
    maxi = FindMaximum[{Kfn[t], t >= 0}, {t, (i - 0.1) inertialPeriod}];
    ti = t /. Last[maxi];
    tMaxVdup = Append[tMaxVdup, ti];
  ];(*end for*)
  tKMaxV = tMaxVdup;
  ntKMaxV = Dimensions[tKMaxV][[1]];
  ptsKmaxV = Table[{tKMaxV[[i]], Kfn[tKMaxV[[i]]]}, {i, 1, ntKMaxV}];
  pltKtestV = ListPlot[
    ptsKmaxV, PlotLegends -> "Expressions",
    FrameLabel -> {"seconds", "K"}, Frame -> True,
    PlotStyle -> {Green}, GridLines -> Automatic]
]

Show[pltKtest, pltKtestV]

In the final plot copied below, the GREEN points are {maxi, Kfn[maxi]} points

enter image description here

Green points are NOT maxima! Any hint?

$\endgroup$
4
$\begingroup$

The problem is in the way you set your starting point for the max search. I suggest that you use constrained optimization instead, looking for maxima within ranges of e.g. ($n$ inertial periods ± 0.1).

In other words, replace the following lines of code in your For loop:

...
maxi = FindMaximum[{Kfn[t], t >= 0}, {t, (i - 0.1) inertialPeriod}]; 
ti = t /. Last[maxi]; 
tMaxVdup = Append[tMaxVdup, ti];
...

Your findTamxVdup then becomes:

findTmaxVdup = Module[{i, t, ti, maxi},
  tMaxVdup = {0}; KMaxVdup = {0};
  For[i = 1, i <= iFin - 1, i++,
   (* changed from here... *)
   maxi = NArgMax[{Kfn[t], i inertialPeriod - 0.1 <= t <= (i + 1) inertialPeriod - 0.1 }, t];
   tMaxVdup = Append[tMaxVdup, maxi];
   (* ... to here *)
  ];
  tKMaxV = tMaxVdup;
  ntKMaxV = Dimensions[tKMaxV][[1]];
  ptsKmaxV = Table[{tKMaxV[[i]], Kfn[tKMaxV[[i]]]}, {i, 1, ntKMaxV}];
  pltKtestV = ListPlot[ptsKmaxV, PlotLegends -> "Expressions", 
     FrameLabel -> {"seconds", "K"}, Frame -> True, PlotStyle -> Green,
     GridLines -> Automatic
  ];
 ]

Using the rest of your definitions, then:

Show[pltKtest, pltKtestV]

Mathematica graphics


In fact, I would suggest furher refactoring to avoid the loop altogether and save the value at max as well as the argument at the same time, rather than recalculating the value separately. So you could replace the entire findTmaxVdup = Module[ ... ]; with the following code:

ptsKmaxV = {t /. #2, #1} & @@@
  Map[
   NMaximize[{Kfn[t], # inertialPeriod - 0.1 <= t <= (# + 1) inertialPeriod - 0.1 }, t] &,
   Range[iFin - 1]
  ];

pltKtestV = 
  ListPlot[ptsKmaxV, PlotLegends -> "Expressions", FrameLabel -> {"seconds", "K"}, 
   Frame -> True, PlotStyle -> Green, GridLines -> Automatic];

Show[pltKtest, pltKtestV]
$\endgroup$
  • $\begingroup$ Great! Thank you very much MarcoB: my problem was fixed. When you are in Venice, let me know... $\endgroup$ – Flavio Sartoretto Jun 6 '18 at 19:48
  • $\begingroup$ @FlavioSartoretto Glad it helped. Grazie per l'invito, a presto! $\endgroup$ – MarcoB Jun 8 '18 at 4:18

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.