Skip to main content
Added explanation of answer
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169

Explanation of using WhenEvent

Integrating an ODE long distances along a separatrix is difficult, because the numerical solution can depart rapidly from the true solution due to small errors in the initial condition. One method of improving the accuracy of the initial conditions is to choose initial guesses (bl and bu in the answer above) that bracket the unknown true initial condition, and then systematically reducing the uncertainty in the initial guesses by doing calculations with initial conditions that bifurcate the distance between the guesses. So, it is necessary to stop a calculation when it obviously is departing from the separatrix, and to note whether the trial calculation is departing above or below. In the answer above, the separatrix is expected to be near 1, except at small r. So, {9/10, 11/10} are expected to bracket the separatrix, and WhenEvent is used to stop the calculation, when the solution moves from inside to outside that range. (Merely being outside that range does not stop the calculation, which is why I check for u < 0 to catch cases in which the solution never reaches the desired range in the first place.) For a solution asymptotically approaching 2, use {18/10, 22/10} or something of that sort. Setting these limits may take some experimentation. Ideally, the range selected should bracket the desired solution with only a modest margin of error, because a large margin of error means that more computer time is required to detect when a particular computation is leaving the expected range.

Explanation of using WhenEvent

Integrating an ODE long distances along a separatrix is difficult, because the numerical solution can depart rapidly from the true solution due to small errors in the initial condition. One method of improving the accuracy of the initial conditions is to choose initial guesses (bl and bu in the answer above) that bracket the unknown true initial condition, and then systematically reducing the uncertainty in the initial guesses by doing calculations with initial conditions that bifurcate the distance between the guesses. So, it is necessary to stop a calculation when it obviously is departing from the separatrix, and to note whether the trial calculation is departing above or below. In the answer above, the separatrix is expected to be near 1, except at small r. So, {9/10, 11/10} are expected to bracket the separatrix, and WhenEvent is used to stop the calculation, when the solution moves from inside to outside that range. (Merely being outside that range does not stop the calculation, which is why I check for u < 0 to catch cases in which the solution never reaches the desired range in the first place.) For a solution asymptotically approaching 2, use {18/10, 22/10} or something of that sort. Setting these limits may take some experimentation. Ideally, the range selected should bracket the desired solution with only a modest margin of error, because a large margin of error means that more computer time is required to detect when a particular computation is leaving the expected range.

corrected timing
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169

Convergence is very good for both the solution, u, and the integral, FNB. (The slight irregularity in FNBat large r is due to a slight boundary condition mismatch, which I shall fix as time permits.) The only significant difference in the revised code used here is that (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 replaces FNB[m - 1][r] in eqnNB to improve numerical stability. Note that this computation required 56 hours on my pc. However, mmax was excessively large to assure convergence, and mmax == 1214 could have been run in 34 hours.

Convergence is very good for both the solution, u, and the integral, FNB. (The slight irregularity in FNBat large r is due to a slight boundary condition mismatch, which I shall fix as time permits.) The only significant difference in the revised code used here is that (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 replaces FNB[m - 1][r] in eqnNB to improve numerical stability. Note that this computation required 5 hours on my pc. However, mmax was excessively large to assure convergence, and mmax == 12 could have been run in 3 hours.

Convergence is very good for both the solution, u, and the integral, FNB. (The slight irregularity in FNBat large r is due to a slight boundary condition mismatch, which I shall fix as time permits.) The only significant difference in the revised code used here is that (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 replaces FNB[m - 1][r] in eqnNB to improve numerical stability. Note that this computation required 6 hours on my pc. However, mmax was excessively large to assure convergence, and mmax == 14 could have been run in 4 hours.

improved wording
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169

Addendum: Converged Iterative SolutionAddendum: Converged Iterative Solution

Convergence is very good for both the solution, u, and the integral, FNB. (The slight irregularity in FNBat large r is due to a slight boundary condition mismatch, which I shall fix as time permits.) The only significant difference in the revised code used here is that (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 is used instead ofreplaces FNB[m - 1][r] in eqnNB to improve numerical stability. Note that this computation required 5 hours on my pc. However, mmax was excessively large to assure convergence, and mmax == 12 could have been run in 3 hours.

Addendum: Converged Iterative Solution

Convergence is very good for both the solution, u, and the integral, FNB. (The slight irregularity in FNBat large r is due to a slight boundary condition mismatch, which I shall fix as time permits.) The only significant difference in the revised code used here is that (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 is used instead of FNB[m - 1][r] to improve numerical stability. Note that this computation required 5 hours on my pc. However, mmax was excessively large to assure convergence, and mmax == 12 could have been run in 3 hours.

Addendum: Converged Iterative Solution

Convergence is very good for both the solution, u, and the integral, FNB. (The slight irregularity in FNBat large r is due to a slight boundary condition mismatch, which I shall fix as time permits.) The only significant difference in the revised code used here is that (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 replaces FNB[m - 1][r] in eqnNB to improve numerical stability. Note that this computation required 5 hours on my pc. However, mmax was excessively large to assure convergence, and mmax == 12 could have been run in 3 hours.

added addendum.
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169
Loading
Corrected typo
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169
Loading
fixed typo
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169
Loading
corrected omission of eqnNB, updated sp, added sp1
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169
Loading
More robust starting values for bl, bu; minor other changes
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169
Loading
improved wording
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169
Loading
Source Link
bbgodfrey
  • 63.1k
  • 18
  • 94
  • 169
Loading