A message from our CEO about the future of Stack Overflow and Stack Exchange. Read now.

Hot answers tagged

19

So there is a ton to say about this, and we will actually be putting a paper out that tries to summarize it a bit, but let me narrow it down to something that can be put into a quick StackOverflow post. I will make one statement really early and keep repeating it: you cannot untangle the efficiency of a method from the efficiency of a software. The details ...


16

Note that $\pi/2$ is represented in double precision format in a way that is not exactly equal to $\pi/2$. It's only accurate to about 15 digits. Thus you're starting every so slightly away from the equilibrium position. Since the equilibrium is unstable, it will eventually start moving.


15

I think the two main points have already been made by Brian and Ertxiem: your initial value is an unstable equilibrium and the fact that your numerical computations are never really exact provides the small perturbation that will make the instability kick in. To give a bit more detail how this plays out, consider your problem in the form of a general ...


10

DifferentialEquations.jl library is a library for a high level language (Julia) which has tools for automatically transforming the ODE system to an optimized version for parallel solution on GPUs. There are two forms of parallelism that can be employed: array-based parallelism for large ODE systems and parameter parallelism for parameter studies on ...


8

This is technically still an IVP if you do an appropriate change of variables. Given your time is between $t \in [t^*, 0]$, make a new time variable $\tau = -t$ so that $\tau \in [0, -t^*]$ and you can modify the time derivatives accordingly. This means that you should have the differential equation $\frac{dy}{d\tau} = -f(-\tau, y)$ with $y(\tau = 0) = 0$ as ...


7

This is a way of writing the flow of a vector field. In other words, if you have an ordinary differential equation (ODE) given by $\dot{x}(t) = v(x(t))$ with initial value $x(0) = x_0$, you could formally write its solution as $x(t) = \mathrm{e}^{t v} x_0$. The formal solution coincides with the actual solution when $v$ is a linear mapping (think of the case ...


7

The three best sources of Butcher tableaus are, according to me, the reference books Numerical Methods for Ordinary Differential Equations, 3rd ed., J.C. Butcher, Solving Ordinary Differential Equations I, 2nd ed., Hairer, Norsett and Wanner and Solving Ordinary Differential Equations II, 2nd ed., Hairer and Wanner. You will have to browse through the ...


7

But clearly, this is not the case as my programs do come up with (an approximate) solution though. I believe you did not continue the integration until you see that your integration is not convergent and is not bounded. I could rewrite your system of ODEs as: $$\dot{x_{1}} = x_{2}$$ $$\dot{x_{2}} = -kx_{1}$$ Or in matrix form: $$\dot{X} = AX$$ Where: ...


7

The 14th order methods due to Feagin can be found in DifferentialEquations.jl. An example using them with 128-bit floating point arithmetic is as follows: using OrdinaryDiffEq, DoubleFloats function lorenz(du,u,p,t) du[1] = 10.0(u[2]-u[1]) du[2] = u[1]*(28.0-u[3]) - u[2] du[3] = u[1]*u[2] - (8/3)*u[3] end u0 = [1.0;0.0;0.0] tspan = (0.0,100.0) prob = ...


6

Two systematic ways of smoothing a function $h$ would be: 1. Join the piecewise smooth parts of your function using Hermite interpolation so that the derivatives are matched to your satisfaction. 2. Convolve your function $h(x)$ with a heat kernel of the form $f(x) = \frac{\exp\left\{-\frac{x^2}{2 \sigma^2}\right\}}{\sqrt{2 \pi \sigma^2}}$ so that instead ...


6

Do you have more suggestions to improve the precision and accuracy of the method ? What are the cause of such errors from the mathematical point of view ? For instance if we have different time scale in the dynamics ... You need to use a reversible ODE solver method if you want to do this. I actually recently showed in a blog post that there are many cases ...


6

There are explicit, symplectic methods for certain types of Hamiltonian problems. For example, the symplectic Euler method \begin{align} p_{n+1} &= p_n - h H_q(p_{n+1}, q_n) \\ q_{n+1} &= q_{n} + h H_p(p_{n+1}, q_n) \end{align} is symplectic, see e.g. Theorem 3.3 on p. 189 in the book by Hairer, Wanner and Lubich (see full reference below). For ...


5

What you are looking for is called "bootstrapping". It is a common problem of all multistep ODE integrators and is discussed in many books on the topic. Among your options are to use a lower-order method with smaller time step, or to use a one-step method of higher order for the first few steps (e.g., a Runge-Kutta method).


5

You can't. Nonlinear systems of equations are in general not solvable exactly. What you need to do is to use a method to solve nonlinear systems, of which there are of course quite a lot: A simple approach would be to use $D(U)\approx D(U^{n-1})$, where $D^{n-1}$ is the solution of the previous time step. A possibly smarter approach would be to use $D(U)\...


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


4

The initial assumption was that the initial position was at a stable equilibrium (i.e., a minimum of the potential energy) with zero kinetic energy and the system started moving away from the equilibrium. Since physically it can't happen (if we consider classical mechanics), two things came to my mind: The first one is that maybe the initial position is: ...


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

since you're dealing with a BVP, it's not a good choice to reduce to a first order system. That's because you should use a finite difference approach. Given a uniform grid from $0$ to $1$, with $N$ equally-spaced knots, i.e. $x_i=x_{i-1}+nh$, where $h$ is the discretization step, you can discretize the second and the first derivative (assuming knowledge of ...


4

I do not believe that a direct solve is possible in your situation, there are simply too many applications with too many users which would benefit too much if it was possible. In particular, the problem of computing, say, the transfer function $T(s) = C(A - sI)^{-1}B+D$ of a linear time invariant system $$ \begin{align} x'(t) &= Ax(t) + Bu(t) \\ y (t) &...


4

Look at the components of the forces calculated in your functions. You will probably find they are never exactly zero, because as other answers have said, you can't represent the value of $\pi$ exactly in computer arithmetic, and the routines that calculate trig functions are not exact either. Eventually, the tiny forces (probably of order $10^{-16}$ at ...


3

$$G(u)=\delta u''+u(u'-1) =0 \\ u(0)=a, u(1)=b$$ If you wish to calculate implicitly the Jacobian before applying any discretization scheme on your PDE first (e.g. so you can reuse the code of a linear solver already available to you), you need the linearization of the nonlinear operator $G(u)$. This is the Frechet derivative of $G$, defined as $$G_u \Delta\...


3

I did not yet succeed to use the sde.sim() function of the SDE package, however, I succeeded to solve the system (with and without noise) using the suggestions of Chris and the diffeqr package in R. library(plotly) # Lotka-Volterra Model (SDE without noise) f <- function(u,p,t) { du1 = p[1]*u[1]-p[2]*u[1]*u[2] du2 = p[3]*u[1]*u[2]-p[4]*u[2] return(...


3

The common way of tracking such properties is indeed to add an additional ODE to the system, here: $$\frac{df}{dt}=a(t)+b(t)+x(t)+z(t)$$ with initial condition $$f(0)=0.$$ If this is not possible, then a smart choice of output times e.g. quadrature points in t for an integration afterwards is recommend.


3

The basic idea is that You use the estimated error given to you (cheaply) by the embedded methods; You use a metric to define acceptance using a user-defined relative and absolute tolerance; Based on the order properties of the code and this metric, a new step size is computed; You avoid large differences in step sizes from one step to the next by limiting ...


3

You can do this computationally for both one-sides (first order similar to Euler) and center (second order) differences easily using the code below. Basically the error due to floating point takes over much sooner than machine epsilon would suggest even for double precision. The case is even worse for single precision of course (you can change the below ...


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

You will need to write your problem such that the unknowns are a single vector, not a matrix. In your example with $N=2$, you will have an unknown vector $x(t)$ of size $20\times 1$ (not a matrix of $10\times 2$). You will solve a problem of the following shape $$\dot{x}(t) = A(t)\, x(t), \mathrm{ with }\ x \in \mathbb{R}^{n},\ A(t) \in \mathbb{R}^{n\times n}...


3

You should be able to set this up as some sort of generalized eigenvalue problem like $$Au=\lambda Bu,$$ where $Au = (du')'+fu$ and $Bu = -(cu')'-eu$. I simplified the derivatives a bit since this formulation is more numerically stable. You should now be able to represent both as tridiagonal matrices using standard finite differences and any generalized ...


3

This problem has an invariant which is the total energy $$ E(t) = \frac{1}{2}(\dot{x}^2 + k x^2) = \textrm{constant} $$ As done by AloneProgrammer, write as first order system $$ \dot{x}_1 = x_2, \qquad \dot{x}_2 = - k x_1 $$ In the phase space $(x_1,x_2)$, the solution must stay on an ellipse whose size is determined by the initial energy. Applying forward ...


3

I think you have a problem due to a misunderstanding of how ODE45 works. If I'm understanding your problem correctly, you want to plot the behavior at certain times, but you are concerned that you can't do this with ODE45. This is a misunderstanding of how ODE45 works. ODE45 uses adaptive time-stepping, and will never use the time array you give it. It will ...


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