9

There are many kinds of RK methods which have extensions to exploit linearity. They all use some form of exponential or Lie Group idea (again exponential) to do so. Thus they generally do some form of integrating factor method and then apply the Runge-Kutta method to the IF transformed equation. If you want to see a listing of some of these, you can check ...


7

As @WolfgangBangerth already commented, this is usually referred to as a differential algebraic equation (DAE). These have their own challenges and there are special numerical methods for them. In general, for a system of DAEs you might have fewer differential equations than unknowns, so the algebraic constraint(s) serve to close the system and identify a ...


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

$b_2$ is zero in the Butcher tableau for Dormand-Prince, so you can skip that term. Note that also the code you linked to contains lines b2 = 0.0, b2p = 0.0, which are just there as a reminder because those variables are never used.


4

You are almost there, just put $t$ under the square root, and it will become $ \displaystyle\int_0^{1/1095} {\frac{dt}{\sqrt{(1+644.153t)(4.17 \cdot 10^{-5} + 0.145 t)}}}, $ which eliminates the singularity at $t=0$


3

A few things jump out at me from your code as potential problems. You seem to be using Newton's method separately for the $y_1$ and $y_2$ variables. This is not the same as using Newton's method for a nonlinear system involving both $y_1$ and $y_2$. For the full Newton method, you'll be solving a 2x2 linear system at every step. You'll have to calculate not ...


3

Specific explanation: Consider the ODE $\dot{y}=f(t,y)$. The exact solution $y(t)$ verifies $$y_{n+1}=y_n + \int_{t_n}^{t_{n+1}}f(t, y(t)) dt$$ for each time step. To numerically compute the integral, one of the quadrature rules can be used, which consist of computing the weighted average of some selected function values on the domain of integration (see ...


3

A search for the specific coefficients listed led me to the method ROS3PRL from J. Sieber, Konvergenzanalyse und Numerische Tests für die Prothero–Robinson–Gleichung (Master thesis), TU Darmstadt, 2014. I can't seem to find this thesis online, but the method is mentioned in the following which may be of interest. Rang, Joachim. "Improved traditional ...


3

Well, in answer to the question about the use of higher-order RK methods, I have a few things to add to the discussion: As an example, it is often helpful to plot an "efficiency diagram" for several different methods being applied to a given initial value problem (using a variety of stepsizes). On such a diagram, one can plot (for each method) how the ...


2

There's a closed form solution. For $a,b,c,d,x_0 > 0$, $$ \int_{x_0}^\infty \frac{dx}{x\sqrt{(a+bx)(c+dx)}} = \\ \frac{\log(2 \sqrt{b d (a + b/x_0) (c + d/x_0)} + a d + b c + 2 b d/x_0) - \log(2 \sqrt{b d a c} + a d + b c)}{\sqrt{b d}} $$ a <- 1 b <- 644.153 c <- 4.17e-5 d <- 0.145 x <- 1/1095 (log(2*sqrt(b*d*(a+b*x)*(c+d*x)) + a*d + b*c + ...


2

The reason you don't find something on integrators of higher order ODEs is that you can always convert them to first order ODEs with auxiliary variables. A equation like $$ \frac{d^2x}{dt^2} = f(t, x, \dot{x}) $$ and $\dot{x}(0)=\dot{x}_0, x(0)=x_0$ can be converted to with an auxiliary variable $\xi$ $$ \frac{dx}{dt} = \xi \\ \frac{d\xi}{dt} = f(t, x, \xi) $...


1

I would check what version of scipy you are using. DOP853 was introduced relatively recently in 1.4.0. In 1.4.1, I see DOP853 listed appropriately. >>> import scipy >>> scipy.__version__ '1.4.1' >>> from scipy.integrate import solve_ivp >>> def exponential_decay(t, y): return -0.5 * y >>> sol = solve_ivp(...


1

It is useful, as a first step, to analyze method-of-lines (MOL) treating the time integration as exact. For example, suppose we are solving the 1D advection equation $ \frac{\partial{n}}{\partial{t}} =-c \frac{\partial{n}}{\partial{x}}, $ where $c$ is the advection speed. Then, for example, for central difference, we'd have $ \frac{d}{dt} n_i = - \frac{c}{2 ...


1

First one should verify the numerical solution for $a(t)$ against the analytic solution. $ d_t a = 1/a + 1/a^2 = \frac{a+1}{a^2} $ Use $b=a+1$, then the ODE becomes $ \frac{(b-1)^2}{b} db = dt, $ from which we find the solution $ b^2/2 - 2 b - \ln(b) = t $ This relation (defined up to a constant to satisfy the initial conditions) is not invertible but ...


1

If the constants are positive then $y$ is positive and setting $u=y^2$ you get $$u'=2A\sqrt{B\sqrt{|u|}+C}$$ which is a little more regular and allows the initial value $u(0)=0$ without singularity. This equation is still stiff close to the initial point, which should have consequences for the local error with fixed-step RK4 as long as $u$ is close to zero. ...


1

The integration curves have the overall form as shown in the figure below, they have vertical asymptotes at $y=0$, which certainly creates difficulties for a numerical solver integrating from $y=0$ (x does not matter). This does not necessarily mean the problem is ill posed; that would be the case if the integration curves diverge rapidly as the initial ...


1

Your problem is not well posed. Multiply both sides by $|y|$ to obtain the ODE in the form $$ |y(x)| y'(x) = A \sqrt{By(x)+C}. $$ If your initial condition is $y(0)=0$, then at the initial time you have $$ 0 y'(0) = A\sqrt{C}. $$ This means that if $A\sqrt{C}\neq 0$, no value of $y'(0)$ can satisfy the equation. In other words, your problem is not that ...


1

I think that you don't need a change of variable (yourself) for this problem. Quadpack seems to work just fine for it. It uses a Gauss-Kronrod quadrature. I tried it (in Python) and it seems to work. import numpy as np from scipy.integrate import quad fun = lambda x: 1/(x*np.sqrt((x + 644.153)*(4.15e-5*x + 0.145))) inte, err = quad(fun, 1095, np.inf) ...


1

To give a short summary of the answer of @Gabriele, the main point is that in a numerical solution of a second order system $x''=f(x)$ as you have with a first order solver as RK4, you need to transform the problem into a first order system $$ x'=p,\\ p'=f(x). $$ This is also the point of departure of the question, it is only in the implementation that the ...


1

Runge-Kutta methods solve equations of the form $$\dot y = \frac{\mathrm{d} y}{\mathrm{d} t} = F(t,y(t)),$$ where $y$ can be multidimensional. The first step is reconducting your equation to this form. Starting from $$ \begin{cases} \dot x_n(t) = \frac{p_n(t)}{m} = f(p_{n}(t))& \\ \dot p_n(t) = -k \left[\left(x_n(t)-x_{n-1}(t)\right) + \...


1

I'm not sure why you think your system of ODEs cannot be solved by using RK4, but I think the easiest way to do this is using DifferentialEquations.jl. Basically you have $2n$ unknowns and $2n$ ODEs, which should be good to solve as long as you have their initial conditions. I think you could look at an example for system of ODEs here from their official ...


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