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


6

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


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}


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


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

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

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

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

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

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


3

This problem is too small to actually be sparse. Sparse handling has a big overhead because the indexing is not "direct", i.e. you don't necessarily know where the next value will be without branch checking. So you need it to be "sparse enough" that the O(n^3) dense LU-factorization cost shrinking to the purely non-zero terms overcomes ...


3

You are implementing the system $$ \pmatrix{\ddot x_1\\\ddot x_2} = \pmatrix{1&1\\-4&2} \pmatrix{x_1\\x_2} $$ in its first-order version with $v_k=\dot x_k$. It is not surprising that this different system results in different solutions. For your stated system you would have to use def V(u,t): x1, x2 = u return np.array([ (x1+x2), (4*x1-2*x2)...


3

As far as I understand the fist part of your question is regarding the options of the ode15s, and which ones you can use to make the solver more efficient/accurate. Providing the Jacobian to the solver, especially if its reasonably easy to calculate, is a good idea to make things faster and more accurate. Using sparse systems for small systems of equations ...


3

You could use a discretization method such as the finite element or finite difference method with a linearization technique such as Picard or Newton method to solve this problem. A similar question is this. The main idea of these solution techniques is as follows: Pick an initial guess for the solution. Linearize your equation and write an updated solution ...


2

You're printing out the time steps. If you print out the solution, you should find the IC is initialized properly.


2

You have not completely vectorized your code. z(1:10) is the first 10 numbers in the data of z, not the first 10 rows. For that you need to write z(1:10,:) like you have already correctly done on the left side. Why the flip in the matrix orientation occurs between the two calls, that is why z(1:10) inherits the column-ness in the first case and the row-ness ...


2

There is no higher magic necessary, just transcribe into the canonical first-order system, encode the boundary conditions, make a reasonable initial guess of the solution shape and call the BVP solver Pr = 5 def odesys(t,u): F,dF,ddF,θ,dθ = u return [dF, ddF, θ-0.25/Pr*(2*dF*dF-3*F*ddF), dθ, 0.75*F*dθ] def bcs(u0,u1): return [u0[0], u0[1], u1[2]-1, ...


2

There is really no such thing as a good finite difference equivalent to an operator. In the earliest days of scientific computing, the thought was that each differential operator would be replaced by some finite difference expression, and that finite difference operator would be the most accurate one available, usually a central difference. I believe that ...


2

This use of the numerical solver is completely wrong. The numerical ODE solvers are for problems that have a smooth right side. As long as the existence of an exact solution is ensured, they can also be applied to problems where the right side is piecewise smooth, but that will tend to slow down the integration with very small step sizes at the ...


2

You were trying to be too helpful. While it is true that for $u(t)=y(T-t)$ you get the equation $\dot u(t)=-\dot y(T-t)=-F(u(t))$ for an autonomous system, this means that to get the correct solution you now need to integrate the modified bw system forward in time, with increasing time points. The second variant that you used is simply to call the solver ...


2

Are you sure that you have the correct boundary conditions? I think the equation you state for the analytical case, $$H=x/\tanh(x)-1$$ is incompatible with the boundary condition $H(\bar{\epsilon})=\bar{\epsilon}$ where $\bar{\epsilon}$ is large. Maybe this is why you are getting incorrect behaviour for large $x$? I tried solving your equation for the $\...


2

Indeed the problem was that you were trying to calculate your nonlinear term incorrectly and you forgot that: $$\mathcal{F}[(f(x))^{2}] \neq (\mathcal{F}[f(x)])^{2}$$ I changed your code to this: import numpy as np import matplotlib.pyplot as plt from matplotlib.pyplot import cm nu = 1 L = 100 nx = 1024 t0 = 0 tN = 200 dt = 0.05 nt = int((tN - t0) / 0....


1

If you want to know the why and the how, you probably want to watch lecture 21.65 here: https://www.math.colostate.edu/~bangerth/videos.html Which you probably want to do after you've watched lecture 21.6. Disclaimer: I'm the one in these videos.


1

The way I computed the solution in the linked answer is the classical one: I just took a reference solution (assuming the code was correct, i.e. the numerical solution I found is the right one) with a small enough step size $h$, say $h=10^{-9}$. Then I computed the solution with smaller $h$s and for each one of those $h$ I computed the $|| e_h ||$ in a ...


1

I have seen that people use CuPy; and scikit-cuda. But, I don't see that any of those provide differential equations capabilities. Looking around, I found CudaPyInt and it uses PyCuda.


1

The simple example as posted wasn't stiff, but I put it together in Julia anyways for show. I modified it to be a system of 2 PDEs with N=1200 and get in the 10's of ms using ModelingToolkit, LinearAlgebra, BenchmarkTools, DifferentialEquations # Setup matrices N = 1200 mat1 = hcat(zeros(N),Tridiagonal(ones(N-1),-2*ones(N),ones(N-1)),zeros(N)) mat1[1,1] = 1 ...


1

The exact solution is $$x(t)=\frac{3}{1-3t}.$$ This means that the domain of the solution is $(-\infty,\frac13)$. odeint with the integration interval $[0,5]$ stalls at this singularity, as it can, correctly, not move beyond this point.


1

Problem solved: I just have to include the solution of the homogeneous equation into the solution of the differential equation. The solution to the differential equation is $$V_H(r) = -\frac{K}{r} - \frac{r+1}{r} \exp(-2r).$$ By requiring that $\lim_{r \to 0} V_H(r)$ is real we get $K=-1$ and $$V_H(r) = \frac{1}{r} - \frac{r+1}{r} \exp(-2r).$$


1

So the size of the problem doesn't necessarily stop you from using the full Jacobian if you store it in a compressed/sparse storage format, you also should be able to know which residuals are dependent on which unknowns which should make it cheaper as well. That said, just use Jacobian free Newton-Krylov (JFNK) solvers. FGMRES is a really strong linear ...


1

There's a really nice book by Li and Ito on the immersed interface method, which was designed to solve problems like yours. Chapter 2 describes 1D problems. Basically, the finite difference method works well when the coefficients are smooth, but when they're rapidly varying or have discontinuities things can go to hell very quickly. You can make it work by ...


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