47

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 symplectic mean and when should you use it? First of all, what does symplectic mean? Symplectic means that the solution exists on a symplectic manifold. A symplectic manifold is a solution set ...


14

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 coordinates $\mathbf{p}$ and $\mathbf{q}$ and a functional, $\mathcal{H(\mathbf{p},\mathbf{q})}$ such that $$\frac{d\mathbf{q}}{dt}=+\frac{\partial \mathcal{H}}{\...


10

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 correct it, simply write: $dy/dx=ay,$ where $a$ is a constant and in our example, $a=1$ in units of $1/x$.


8

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 SolveODE[ <f'(x, y)>, <Start x>, <Start y>, <End x>, <Step> ] with SolveODE[A + B y + C sin(y), l, m, 10, 0.1] I got the solution curve below. You can vary ...


8

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 velocities. This stands in contrast to the continuous case, where the all frequency components travel with exactly the same speed (c). In the first simulation, ...


6

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 domain $\Omega$ as that function $g: \partial\Omega \rightarrow {\mathbb R}$ that lives on the boundary of $\Omega$ and so that $u(\mathbf x)=g(\mathbf x)$ for ...


5

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)$. In order to achieve the advertised accuracy the right hand side should be a sufficiently smooth function of $y$, and nonlinearity can complicate the question ...


4

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(t)$, but not any of the times before (i.e., no longer contains an integral). In other words, you end up with a system of two ODEs that you can analyze as you ...


4

As noted in the comments, when I ran your code (using LSODE in SciPy) I get >>> result = odeint(f,y0,t) lsoda-- warning..internal t (=r1) and h (=r2) are such that in the machine, t + h = t on the next step (h = step size). solver will continue anyway in above, r1 = 0.2036188800000D+10 r2 = 0.1292689215500D-07 lsoda-- ...


4

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 subset of features which includes explicit parameters and visualization. Here's an example of your equation, assuming that l was the initial time point and ...


4

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 implicit ODE. This falls into the class of problems known as Differential-Algebraic Equations, and that still works even though you don't have any pure algebraic ...


4

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 the shadow Hamiltonian, the resulting system of differential equations are called the modified equations). See also: S. Reich, “Backward error analysis for ...


4

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$. However, note that this expression assumes that we are dealing with conservative forces. But if that's the case, the energy is conserved by definition, which means ...


4

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. Rather than try to summarize, I'll just refer you to chapter 15, section 4 of Hairer, Lubich, Wanner, Geometric Numerical Integration. That book is a must-have ...


4

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^{2} \theta = \frac{1}{k}\nabla \cdot \vec{v}$$ Since, you know the $\vec{v}$ you can approximate its divergence as: $$\nabla \cdot \vec{v} = \frac{\partial v_{...


3

Your model seems to have a single stable critical point that attracts all orbits. You're plotting $(S(t), I(t), R(t))$ as functions of time for specific initial conditions, but actually the easier way to see this directly is to plot the vector field $(\dot S, \dot I)$ in the $S$-$I$ plane like so: Immediately you can see that for this particular choice of ...


3

What you're looking for is a discussion of absolute stability of multistep methods for ODEs. See for example chapter 7 of LeVeque's book on finite difference methods . Once you know the region of absolute stability for your time discretization, it's just a matter of checking that the eigenvalues of D (scaled by the time step) are inside that region. This ...


3

This seems like a difficult problem to solve in general, but here's one idea. Suppose you look for all constants of motion of the form $$ \sum_{i=1}^m \lambda_i \varphi_i(t), $$ where $\lambda_i$ are some unknown coefficients, and $\{\varphi_i\}_i$ is some set of basis functions that can be computed by numerical integration of the ODE. For example, $\...


3

I found that the solution of the PML problem by the method proposed here should not work. Next is a detailed analysis of the problem. I start with the von Neumann stability analysis for the FTCS system corresponding to (eq 3.) on the original question. \begin{eqnarray*} \frac{ w_{i j+1} - w_{ij}}{\Delta t} = \frac{c}{2 \Delta x} (w_{i+1 j} - w_{i-1 ...


3

I would call your problem a eigenvector nonlinear two-parameter eigenvalue problem. It is essentially a combination of two difficult problems: Eigenvector nonlinear eigenproblem $$A(v)v=\lambda v.$$ Two-parameter eigenvalue problem $$(A_1+A_2\lambda +A_3\mu)u=0$$ $$(B_1+B_2\lambda +B_3\mu)v=0$$ These problems have received separate attention in ...


3

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 specifically for implicit methods, the error in solving the implicit equation. The first error is of size $O(h^{p+1})$, the second of size $\sim \mu$, $\mu$ ...


3

The situation seems to be: You have some input function $x$ which more or less follows a model $\dot x=-ax+bu$. There may be noise involved, so the values of $x$ are not exact, and simply computing difference quotients will in general not be close to the right side of the differential equation. To filter out the noise some averaging is required. This you do ...


3

The parameter can be any type, so here I pass in a time-dependent function for p and use it in the differential equation: # Packages library(tidyverse) library(diffeqr) library(JuliaCall) diffeq_setup() # Drift function f <- function(u,p,t){ du1 = p(t) return(c(du1)) } # Diffusion function g <- function(u,p,t){ du1 = 0 # note that there is ...


2

First of all, you have no exact mirroring, which can be clearly seen when superimposing your plots: I suggest to consider the following to get a better idea of the possible mistake: Consider different sets of parameters, in particular ones with a less symmetric dynamics, e.g., $κ=0.2$ and $τ=0.1$. If there is a disagreement between numerical and analytical ...


2

Thanks to TroyHaskin, I've realized that this is in reality a linear problem, since the u(s) values are known from the initial conditions, and is therefore easy to solve in Matlab. Thanks to everyone else for help as well. EDIT: To be more specific, the system can be solved linearly by separating the u(s+1) terms and their coefficients from everything else. ...


2

If you can write down the problem using generic constants then this becomes a parameter estimation problem. The DifferentialEquations.jl docs explains a bunch of methods in more detail, but I'll focus on a simple and standard method exemplified here. That example is about finding biological constants (birth and death rates in Lotka-Volterra) but the idea is ...


2

You cannot take deterministic methods and slap a kind of Milstein correction on them to get higher order. Instead, you should be looking at methods which are specifically designed for integrating SDEs at high order if you want good accuracy. If pure efficiency is what matters and the noise is additive/diagonal, then the adaptive Rossler methods mentioned by @...


2

First of all, it is impossible to intertwine a multi-step Runge–Kutta method and the Milstein–Itō methods for a multitude of reasons that go beyond the scope of this question¹. So the best you can possibly do in is: Make a deterministic Runge–Kutta step, ignoring the noise term. Apply noise (in a Milstein–Itō fashion). Go to 1. This has the problem that ...


2

You just add the diffusion along the other dimensions. This superposition from orthogonal directions makes some sense, as they are independent. So, going by wikipedia for Fick's second law of diffusion in 1D: $$ \frac{\partial \psi}{\partial t} = D \frac{\partial^2 \psi}{\partial x^2} $$ We extend it to 2d as: $$ \frac{\partial \psi}{\partial t} = D \...


2

This is a tricky question because mathematicians usually ignore units. I actually ran into this myself when making DifferentialEquations.jl compatible with units. In this notebook, I use Unitful.jl's unitful-arithmetic. Essentially, you just use Julia's string macros to make u"s" put units in seconds and implicit multiplication to make 2u"s" mean "two ...


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