74 votes
Accepted

What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?

Let me start off with corrections. No, odeint doesn't have any symplectic integrators. No, symplectic integration doesn't mean conservation of energy. What does ...
user avatar
17 votes

What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?

To complement Chris Rackauckas answer, to state some of the mathematical nonsense as well as some stuff you almost certainly know, a dynamical system is Hamiltonian if there is a description with ...
user avatar
  • 2,189
17 votes
Accepted

CUDA & Python for numerical integration and solving differential equations

Julia's DifferentialEquations.jl is all GPU-compatible. If you make your arrays GPU-based arrays, then the solver recompiles to be all on the GPU (no data transfers). For example: ...
user avatar
10 votes
Accepted

4th order Runge-Kutta for $y' = y$

The answer is quite simple. You are already comparing apples and oranges in the first equation. Garbage in, garbage out. The equation $y'=y$ if written properly is $$dy/dx=y.$$ Do you see it now? To ...
user avatar
  • 216
10 votes

How to solve a second order differential equation (diffusion) with boundary conditions using Python

I have found that I must keep the value of dt near dx or the results become unstable This behavior you have noticed is known as the Courant–Friedrichs–Lewy (CFL) condition. Indeed, there are ...
user avatar
9 votes
Accepted

Numerical computation of Lyapunov exponent

It is not the Jacobian of $f$ that you need to use to propagate the local basis, but the Jacobian of the propagator of the numerical method. That is, if $$ v_1=v_0+\Phi_f(v_0,dt)dt $$ then $$ U_1=(I+...
user avatar
  • 3,511
8 votes
Accepted

Visualizing the solutions of the Differential equations by varying different parameters

You can try Geogebra (it is free). With SolveODE command and sliders you can do what yo want. For the usage of SolveODE command see. For example by using following command ...
user avatar
  • 643
8 votes

Finite difference for 1D wave equation: why the spike initial data results in a noisy output?

I think you're probably seeing artifacts that are due to numerical dispersion. In brief, in the discrete case different (spatial) frequencies of a wave function will propagate at different phase/group ...
user avatar
  • 4,286
7 votes
Accepted

What is the case of trade-off in different Runge-Kutta methods

There are so many Runge Kutta methods, including Dormand-Prince 45 Cash-Karp 54 Fehlberge (sic) 78 Is there any comparison between them? Well, sure. Here are some traits to ...
user avatar
7 votes

Specifying ode solver options to speed up compute time

Julia's DifferentialEquations.jl has a lot of tooling for automatically deriving (sparse) matrices. For more information, see the JuliaCon 2020 video on Auto-Optimization and Parallelism in ...
user avatar
6 votes

Some questions on Trace (operators) on the boundary in the context of PDEs

I think the quote from Wikipedia is misleading (and I should put that on my list of things to fix :-) ). You can think of the trace of a function $u: \Omega \rightarrow {\mathbb R}$ that lives on a ...
user avatar
6 votes

Solving coupled differential equations in Python, 2nd order

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{...
user avatar
  • 6,046
5 votes
Accepted

Can Runga-Kutta method be used to solve non-linear differential equations?

Runga-Kutta schemes are multistage recipes for numerical discretizations of temporal derivatives. That is they tell to how to solve equations in the form $\dot{\mathbf{y}} = \mathbf{f}(\mathbf(y,t)$. ...
user avatar
  • 2,189
5 votes
Accepted

Numerical integration in 2D

You can convert this equation into a Poisson equation and use standard numerical methods to solve it like finite difference to obtain $\theta(x,y)$. If you take a divergence from both side: $$\nabla^{...
user avatar
5 votes
Accepted

Numerical Solution to Rayleigh Plesset Equation in Python

I see at least one important problem. On the right hand side you have a term that looks like $$P_o \left( \frac{\dot{R}}{R} \right)^{3 \kappa}$$ This term is dimensionally inconsistent with the ...
user avatar
5 votes
Accepted

Solving the heat diffusion equation with source term

You are starting from a uniform temperature and you have insulated boundary conditions; so there is no heat conduction occurring. Likewise, your initial mole fraction is also constant in $x$ so that ...
user avatar
  • 5,734
5 votes
Accepted

scipy odeint: excess work done on this call and very sensitive to initial value

Essentially, by combining all the constant factors, your ODE is $$ \frac{dL}{dt} = -A + B L + C L^2 $$ With the initial $L(0)=L_0$ large enough, the positive feed-back of the quadratic term drives ...
user avatar
  • 3,511
5 votes
Accepted

2-DOF Robotic Manipulator Trajectory Tracking Simulation

As far as I understand the fist part of your question is regarding the options of the ode15s, and which ones you can use to make the solver more efficient/accurate. ...
user avatar
5 votes

Numerical integrator for $a'(t)=e^{-a(t)}f(t)$

You want a numerical solution, but this might help you check your computed results. If $a$ satisfies the ODE, you know $e^{a(t)}a'(t) = f(t)$. Integrating you get \begin{align} \int_0^t\, f(\tau)\, d\...
user avatar
5 votes

Solving numerically the 1D Kuramoto-Sivashinsky Equation using spectral methods

Indeed the problem was that you were trying to calculate your nonlinear term incorrectly and you forgot that: $$\mathcal{F}[(f(x))^{2}] \neq (\mathcal{F}[f(x)])^{2}$$ I changed your code to this: <...
user avatar
5 votes

solve_ivp from scipy does not integrate the whole range of tspan

You can examine the sol object to see why the integration failed. It provides the message 'Required step size is less than spacing between numbers.' This usually ...
user avatar
4 votes
Accepted

Finding Numerical Stability of Simple System with Integral Term Due to Low Pass Filter

You can introduce an auxiliary variable $$ y(t) = \int_0^t \exp\left(\frac{-(t-\hat t)}{\tau}\right) x(\hat t) \; d\hat t, $$ which you can differentiate to get on ODE for $y(t)$ that depends on $x(...
user avatar
4 votes
Accepted

Small errors accumulate while solving ODE of motion

As noted in the comments, when I ran your code (using LSODE in SciPy) I get ...
user avatar
4 votes

Visualizing the solutions of the Differential equations by varying different parameters

You can use DifferentialEquations.jl Online to visualize solutions to differential equations without a hassle. It's built using the Julia suite DifferentialEquations.jl, and the online interface is a ...
user avatar
4 votes
Accepted

Integrating a nonlinear ordinary differential equation

You don't have just a first-order ODE so you cannot use an explicit Runge-Kutta method. Because of the square term, you cannot bring this into mass matrix form even. Instead, what you have is an ...
user avatar
4 votes
Accepted

Symplectic Algorithms for Hamilton’s Equations as opposed to just Volume-Preserving

Backward error analysis comes to mind. For example, for the Verlet scheme, you can say that the numerical solution turns out to be the exact solution for a perturbed Hamiltonian system (also known as ...
user avatar
4 votes

Calculate forces on atoms from potential energy of system and position of atoms

Generally speaking, the force acting on particle $i$ is just $-\nabla E(\vec{r}_i)$ (note the minus sign), where $\nabla$ is the gradient operator and $\vec{r}_i$ is the position of particle $i$. ...
user avatar
  • 657
4 votes
Accepted

Wrong results for $2$ stage multistep method $y_{n+2} - y_n = h\left[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_n\right]$

The local error in each integration step is composed of 3 parts: the discretization error of the method the floating point error of actually performing the method in precision-limited data types, and ...
user avatar
  • 3,511
4 votes
Accepted

Symplectic linear multistep method?

Yes, but you have to mean symplectic on a higher-dimensional phase space than your original problem that includes previous steps too. As I understand there are also some subtle stability issues too. ...
user avatar
4 votes
Accepted

Using Implicit Euler with second order differential equations

You have to write your second order equation as a system of two first order equations. Let $y' = v$, then your equation $$ y'' + \omega^2 y = 0 $$ becomes $$ \begin{pmatrix} y' \\ v' \end{pmatrix} =...
user avatar
  • 1,198

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