Hot answers tagged

9

There are many kinds of RK methods which have extensions to exploit linearity. They all use some form of exponential or Lie Group idea (again exponential) to do so. Thus they generally do some form of integrating factor method and then apply the Runge-Kutta method to the IF transformed equation. If you want to see a listing of some of these, you can check ...


7

As @WolfgangBangerth already commented, this is usually referred to as a differential algebraic equation (DAE). These have their own challenges and there are special numerical methods for them. In general, for a system of DAEs you might have fewer differential equations than unknowns, so the algebraic constraint(s) serve to close the system and identify a ...


6

The first step is to transform the second order equation to a set of two coupled first order equations. Define an auxiliary function $u(T) = \frac{dr(T)}{dT}$. This results in the system $$\begin{align} \frac{du}{dT} &= k-(1-\frac{5}{r})(3+\frac{2}{r^2}) \\ \frac{dr}{dT} &= u\\ \frac{d\phi}{dT} & = \frac{1}{r^2} \end{align} $$ Now you have a ...


6

I got your point: you have the ODE $y''=y'+x$, and we can see that $y(x)=-\frac{x^2}{2} - x$ solves the problem. As you want to integrate it numerically, you want to set two initial conditions, so you imposed $y(0)=0$ (and that's fine) and also another one $y(h)$, which is wrong because it's not the required information you need to solve an ODE, mainly ...


5

The second formula is just the velocity verlet, and it's correct but if you adapt time steps then it's not symplectic. In a separate answer I describe in quite detail that symplecticness is a global property of the integration, it's not a stepwise property. Because of this, local error estimates and local changes of $\Delta t$ do not necessarily preserve the ...


5

Consider the vector $y_i$ for a given time $t_i$. We are looking for the linear operator $\mathcal{B}$ that maps $y_i$ to $y_{i+1}$ by noting that \begin{align} \dot{y}&= \mathcal{A} y\\ \frac{y_{i+1}-y_i}{\Delta t}&=\mathcal{A} y_i\\ \Rightarrow y_{i+1}&=(\mathcal{A}\Delta t+I) y_i\\ y_{i+1}&=\mathcal{B} y_i \end{align} For a set of vectors $...


4

Parareal is fairly effective for parabolic problems, but doesn't work well for hyperbolic problems. So it wouldn't be effective for the advection part, but could be effective for the diffusion part. Making it work well for hyperbolic problems is an area of ongoing research. Note that the parallel efficiency of parareal is quite poor, so if you haven't yet ...


4

You don't actually need co-routines for this -- you can achieve the same by just using regular callbacks into user code. For example, in the C programming language, the ODE solver routine would have an extra argument that denotes a pointer to a function. User code would call the integrator with a pointer-to-function for a user function that the ODE ...


4

You're using a numerical method to approximate the solution to the continuous dynamical system. If you've done this carefully, that approximate solution could be adequate for investigating the properties of the continuous dynamical system. You could also, if you wanted to, analyze the numerical approximation as a discrete dynamical system. The key word ...


4

This is not an answer to your question, but more of an observation: More often than not, an ODE is "stiff" or a linear system is "nearly singular" because of a mistake either in deriving the equation to be solved, or in implementing it. Trying to find a way to solve what you have is then just a way to paper over the problem. If you were to find a way to ...


4

Is there a way to solve for multiple initial conditions without resorting to a slow python for loop? I think the solution to your problem might be parallelization, not vectorization, unfortunately. While Python is going to be slower than C or C++, the for loop in question isn't actually a part of the computation, it's a different problem numerically because ...


4

The correct dynamic equations in the polar coordinates should be $ \dot{v_r} = \omega^2 r - \alpha/r^2 \\ \dot{\omega} = - 2 v_r \omega /r\\ \dot{\theta} = \omega \\ \dot{r} = v_r $ Here is the fixed Python code: from math import * import numpy as np from matplotlib import pyplot as plt from scipy.integrate import odeint def vec(w,t): r,vr,theta,omega=...


4

Many numerical tips and theoretical explanations can be found in this book from Hairer and Wanner: https://www.springer.com/gp/book/9783540566700 In this book, a strategy is described, which uses a time step such that the relative variation of the solution during the first time step is below a certain threshold if you were using explicit Euler (omitting the ...


4

The result you see is due to the long integration interval and the accumulation of truncation errors. Each step with its truncation error will switch to a slightly different exact solution. Now the exact solutions are in general $$ y(x)=-\frac12x^2-x+C+De^x $$ Due to the first order approximation of the first derivative, the coefficients $C$, $D$ will grow ...


3

For debugging the code, there is a set of analytic solutions here for several reduced models corresponding to subsets of terms on the right-hand side. These analytic solutions have to be reproduced by the code. Verification testing of this kind is a standard practice for debugging simulation models. Reduced model 1: $ m \ddot{x} = - \gamma \dot{x} $ Solution:...


3

There are integrators which expose a more fine-grained control to the user, e.g., scipy.integrate.ode: from scipy.integrate import ode def f(t,y): return y[0] initial = [1] ODE = ode(f) ODE.set_integrator("dopri5") ODE.set_initial_value(initial) for time in range(1,10): print(ODE.integrate(time)) Instead of coroutines, this provides the ...


3

Your issue is that you are modifying y in your function. If you return, [u,v] instead of resetting the values in y, it seems to run just fine. I believe this is because the solver will store and use intermediate values elsewhere, but your function is having the side effect of altering these stored values.


3

It's hard to get around the allocations implicit to SymPy in this case. It wants to allocate the matrix, so the easiest thing to do would be, as you show, build individual scalar functions. But then composing those together can be a bit of a hassle, since you don't want to put them into an array since they are all different types and that would then ruin the ...


3

Well, in answer to the question about the use of higher-order RK methods, I have a few things to add to the discussion: As an example, it is often helpful to plot an "efficiency diagram" for several different methods being applied to a given initial value problem (using a variety of stepsizes). On such a diagram, one can plot (for each method) how the ...


3

I would recommend 3 things to pursue: Use ode15s (a stiff solver). This will allow for better resolution near the jumps Rescale the problem. Your coefficients are very large and if they are changing by several orders of magnitude, this is likely causing undue conditioning problems with your solver so rescaling then undoing the scaling after the ode solve is ...


3

Shortly after reaching $t=0.4$ the system has a pole, in component $z_4$ with component $z_3$ following. Such behavior is to be expected for systems with positive quadratic terms, similar to $y'=y^2$. This leads to overflow and the interruption of the integration. The remaining array elements are/remain filled with the random garbage that they contained at ...


3

So I've tried to compute successive time derivatives for the position $x_1$ of a mass from a coupled mass-spring system (adapted from a previous question). The system is linear: $d_t x = Ax$, hence the analytical $n$-th time derivative is $\frac{d^n x}{dt^n} = A^n x$, which permits an easy estimation of the error. In a "try hard" approach which I ...


3

A few things jump out at me from your code as potential problems. You seem to be using Newton's method separately for the $y_1$ and $y_2$ variables. This is not the same as using Newton's method for a nonlinear system involving both $y_1$ and $y_2$. For the full Newton method, you'll be solving a 2x2 linear system at every step. You'll have to calculate not ...


3

In a one-step method $$ y_{n+1}=y_n+h\Phi_f(x_n,y_n,h) $$ one gets a truncation error for the exact solution $$ y(x_{n+1})=y(x_n)+h\Phi_f(x_n,y(x_n),h)+h^{p+1}\tau(x_n) $$ For the error propagation of $e_n=y_n-y(x_n)$ this gives approximately $$ e_{n+1}=e_n+h\partial_y\Phi_f(x_n,y(x_n),h)e_n-h^{p+1}\tau(x_n)+... $$ The error is then of order $p$, $e_n=c(x_n)...


2

The phenomenon you're dealing with here is called a "stiff system". Check here for a full explanation: http://www.scholarpedia.org/article/Stiff_systems Matlab’s ode45 is a non-stiff solver, thus ill-suited for your problem. You should try to use one stiff solver. Matlab has 4 stiff solvers but other tools like Scipy’s solve_ivp also possess non-stiff ...


2

This is a system of first order differential equations, not second order. It models the geodesics in Schwarzchield geometry. In other words, this system represents the general relativistic motion of a test particle in static spherically symmetric gravitational field. In general, there is a third equation for how coordinate time is related to proper time. ...


2

Another solution is to solve with solve_ivp and use the dense_output option which allows interpolating between solution steps: import numpy as np from scipy.integrate import ode, solve_ivp import matplotlib.pyplot as plt import warnings def f(t,y): l = 1 m = 1 d = 1 g = 9.8 return [y[1], -np.sin(y[0])*g/l-y[1]*d/m] y0, t0 = [np.pi/2, 0], 0 ...


2

Not quite a package but the example at the end of the Revised6 Report on the Algorithmic Language Scheme implements a classical Runge-Kutta integrator as a function that returns 'an infinite stream of states'. This is done in a few lines of code and is easy to port to for example Python using its coroutines, generators, itertools, &c. I think that ...


2

Due to your formulation, I call $X(z) = \begin{bmatrix} x(z) \\ p_{x}(z) \end{bmatrix}$ so your ODE is written in matrix form: $$X^{'}(z) = C(z) X(z)$$ Where: $C(z) = \begin{bmatrix} 0 & A(z) \\ B(z) & 0 \end{bmatrix}$. Your general formula by using backward Euler method is: $$\frac{X(z+\Delta z) - X(z)}{\Delta z} = C(z+\Delta z) X(z+\Delta z)$$ ...


2

With the functions $y_1 := f$, $y_2 := f'$, $\pmb{y} := (y_1,y_2)^{\top}$, we obtain an initial-value problem with an autonomous first-order system: $$ \pmb{y}' = \left( \begin{array}{c} y_2\\ -a \sin(y_1) \end{array} \right) =: \pmb{f}(\pmb{y}), \quad \pmb{y}(0) = \left( \begin{array}{c} b\\ c \end{array} \right) =: \pmb{y}_0. $$ We now choose the Euler ...


Only top voted, non community-wiki answers of a minimum length are eligible