58

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


12

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: using OrdinaryDiffEq, CUDA, LinearAlgebra u0 = cu(rand(1000)) A = cu(randn(1000,1000)) f(du,u,p,t) = mul!(du,A,u) prob = ODEProblem(f,u0,(0.0f0,1.0f0)) # Float32 is better on ...


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


7

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+D_v\Phi_f(v_0,dt)dt)U_0 $$ In a pinch you can of course replace the derivation of the method step with the derivation of the Euler step, as the Lyapunov ...


7

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 DifferentialEquations.jl. Combined with the orders of magnitude acceleration commonly seen over the MATLAB solvers, this might be a good option for you and is a quick ...


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

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\tau &= \int_0^t e^{a(\tau)}a'(\tau)\, d\tau \\ &= \int_{a(0)}^{a(t)} e^a da \\ &= e^{a(t)} - e^{a(0)}. \end{align}


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

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


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

You're approaching this the wrong way by using the product rule of differentiation. Rather, making use of the fact that $$ \frac{\partial u(x_i)}{\partial x} \approx \frac{u_{i+1/2}-u_{i-1/2}}{\Delta x}, $$ the first step in finding what you need to do is to recognize that $$ z(x_i)\frac{\partial u(x_i)}{\partial x} \approx z_i \frac{u_{i+1/2}-u_{i-1/2}...


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.


4

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

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 the equation towards a dynamic blow-up in finite time, as you observed with the values getting very large. This is a property of the exact solution, so it is ...


4

Some things I can think of: use sparse matrices for Matrix1 and Matrix2 to speed up the computations of dZand dY use larger integration tolerances reltoland abstol, especially if your are searching for steady-state solutions and/or do not need a precise resolution of the transient dynamics of your system. you are solving with ode15s, which is an implicit ...


4

The go-to reference on this topic is the extensive book by Hairer, Lubich, & Wanner (HLW). The kind of property you are dealing with is known as a quadratic invariant, since $\|x(t)\|^2$ is constant. These are covered in Section IV.2 of the book. Any quadratic invariant will be preserved by a Runge-Kutta method if the method coefficients satisfy $$ ...


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