55

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 which ...


15

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 ...


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 ...


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 ...


5

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_{...


5

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 other terms in the brackets, which have dimensions of a pressure. This term should actually be $$P_o \left( \frac{R_o}{R} \right)^{3 \kappa}$$ The paper you link ...


5

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 the heat input is uniform along the length. So the fact that the temperature and mole fraction don't change as a function of $x$ is exactly what you should ...


5

I am not sure this is possible with the Python libraries since they are using Fortran under the hood and that can't be easily recompiled, but the Julia DifferentialEquations.jl JIT compile specializes the solvers based on the number types that you give it. Here's a demonstration of some weird types like rational numbers, MPFR BigFloats, and ArbFloats (based ...


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 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} = \begin{pmatrix} v \\ -\omega^2 y \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix} \begin{pmatrix} y \\ v \end{pmatrix} $$ If you ...


4

Is there any out-of-the-shelf method I can use in this case to solve the equation without fully calculating the Jacobian? Using any integrator for stiff ODEs where the implicit equation is solved via Newton-Krylov methods will only require Jacobian-vector products. So something like DASKR, IDA, etc. can be made to do this.


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

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 \...


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 ...


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