New answers tagged

1

To explain further what I wrote in the comments: Your function has the form F(t,x,co,xo) where x is the state, co the constants for the other objects and xo the dynamical data for the other objects. Variant 1: You will get an order 1 approximation out of RK4 if you use the external dynamic in the RK4 loop as (in python notation) dx1 = dt*F(t[k] , x[k] ...


2

For the first question: You should look into continuation. I don't know python so i can't write the code for you, but basically, you change the algorithm into something like this: search for some solutions for phi = -2. use those solutions as initial guesses when looking for phi = -1.9 repeat step 2 for the rest of the interval Because you are using the ...


3

You may take a look at Peter Gottschling's Discovering Modern C++, especially chapter 7, where Mario Mulansky (one of the authors of odeint) implements a generic ODE Solver (using Runge-Kutta algorithms and C++11/14). To the best of my knowledge, the second version of the book will be published within the next months and covers C++17 and C++20. Another ...


2

Browsing the paper I can confirm that all integrations of the trajectories is done with a fixed time step. The authors compute approximations $A_h(t)$ for all $t \in \Sigma_h$ where $$\Sigma_h = \{ nh \: : \: 0 \leq nh \leq T\}.$$ They use at least two different values of $h$, namely $h_1=0.125$ days and $h_2 = 8$ days. It is important to recognize that $$ \...


2

As the integrand (in the answer you quote) is periodic and fairly smooth, you can evaluate it numerically using the trapezoidal rule. For such integrands, the convergence is exponential, so you shouldn’t need too many points. Here is a nice explanation. I think it is quite feasible to do this in Excel. Here's some matlab code that shows the rapid convergence ...


2

You have not completely vectorized your code. z(1:10) is the first 10 numbers in the data of z, not the first 10 rows. For that you need to write z(1:10,:) like you have already correctly done on the left side. Why the flip in the matrix orientation occurs between the two calls, that is why z(1:10) inherits the column-ness in the first case and the row-ness ...


2

odenumjac calls your function in a vectorized manner it seems, and your function is not vectorized. You can easily change that by changing the second index of f in your function to : instead of 1, for instance: f(10,:) = 2*(x(end-1,:) - x(end,:)); I thought the setting joptions.vectvars=1 would not allow the vectorised call (see one of your other questions). ...


3

The effect is due to quantization noise and/or aliasing, you are not computing with the true minima of the radius along the orbit, but with the closest sample point. This means that the sample point for the half-step integration can lie in-between. Close to the minimum this results in an $O(e·dt^2)$ ($e$ the eccentricity) difference between the computed ...


Top 50 recent answers are included