# Numerics

**URL:** https://discourse.julialang.org/c/domain/numerics/20.md?page=11

[Latest](https://discourse.julialang.org/latest.md) · [Categories](https://discourse.julialang.org/categories.md) · [Tags](https://discourse.julialang.org/tags.md)

**Page:** 12

---

## [Package for matrix polynomials](https://discourse.julialang.org/t/package-for-matrix-polynomials/107687)

<div class="topic-metadata">

**Author:** [@fph](https://discourse.julialang.org/u/fph)\
**Replies:** 2\
**Last updated:** [December 15, 2023, 10:23pm UTC](https://discourse.julialang.org/t/package-for-matrix-polynomials/107687 "2023-12-15T22:23:22Z")

</div>

What package do you recommend to work with matrix polynomials, i.e., expressions of the form A(x) = A\_0 + A\_1 x + \\dots + A\_d x^d? I will use mainly the coefficients to do linear algebra, so I am interested in a package…

---

## ["Warning: Instability detected. Aborting" while using lyapunovspectrum](https://discourse.julialang.org/t/warning-instability-detected-aborting-while-using-lyapunovspectrum/107536)

<div class="topic-metadata">

**Author:** [@Daniel\_Montesinos\_Ca](https://discourse.julialang.org/u/Daniel_Montesinos_Ca)\
**Replies:** 7\
**Last updated:** [December 13, 2023, 7:15pm UTC](https://discourse.julialang.org/t/warning-instability-detected-aborting-while-using-lyapunovspectrum/107536 "2023-12-13T19:15:30Z")

</div>

Hey, I’m trying to study the lyapunov spectrum for a 4D dynamical system using the lyapunovspectrum function from DynamicalSystem.jl but depending on the parameters I get “Warning: Instability detected. Aborting” and I …

---

## [Accurate derivative for tanh](https://discourse.julialang.org/t/accurate-derivative-for-tanh/107568)

<div class="topic-metadata">

**Author:** [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Replies:** 13\
**Last updated:** [December 13, 2023, 2:43pm UTC](https://discourse.julialang.org/t/accurate-derivative-for-tanh/107568 "2023-12-13T14:43:45Z")

</div>

When debugging a computation, I found that the derivative of tanh as defined in DiffRules.jl suffers from catastrophic cancellation for pretty mild values (eg outside 20 in absolute value, examples below). I propose a s…

---

## [All positive eigen values but false in isposdef?](https://discourse.julialang.org/t/all-positive-eigen-values-but-false-in-isposdef/107393)

<div class="topic-metadata">

**Author:** [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Replies:** 2\
**Last updated:** [December 11, 2023, 1:25am UTC](https://discourse.julialang.org/t/all-positive-eigen-values-but-false-in-isposdef/107393 "2023-12-11T01:25:28Z")

</div>

Hello, I was curious about my Julia code spitting out negative variance for some variations in sandwich form vcov matrices (in Fixed Effects, clustering etc.). Let’s say my matrix of control variables, X. I found out t…

---

## [Performance of lu(A)\\B versus inv(A)\*B to solve AX=B](https://discourse.julialang.org/t/performance-of-lu-a-b-versus-inv-a-b-to-solve-ax-b/107331)

<div class="topic-metadata">

**Author:** [@Pablo\_Marchant](https://discourse.julialang.org/u/Pablo_Marchant)\
**Replies:** 3\
**Last updated:** [December 8, 2023, 4:48pm UTC](https://discourse.julialang.org/t/performance-of-lu-a-b-versus-inv-a-b-to-solve-ax-b/107331 "2023-12-08T16:48:56Z")

</div>

Hi all, I’m working in a problem where I need to solve various Ax=b systems, where A is a dense (but usually not that large) matrix. Initially I was doing this by directly computing inverses and multiplying. Awful to do…

---

## [Problems saving/opening ODESolution from SciML with JLD2](https://discourse.julialang.org/t/problems-saving-opening-odesolution-from-sciml-with-jld2/106958)

<div class="topic-metadata">

**Author:** [@yhchang96](https://discourse.julialang.org/u/yhchang96)\
**Replies:** 12\
**Last updated:** [December 2, 2023, 12:49am UTC](https://discourse.julialang.org/t/problems-saving-opening-odesolution-from-sciml-with-jld2/106958 "2023-12-02T00:49:40Z")

</div>

I define a custom struct to hold the parameters for the coupled ODEs I’m solving with OrdinaryDiffEq.jl, which typically stores various numbers and arrays (1D or 2D), as well as vectors of dualcache vectors in par.dualca…

---

## [Huge computation error in Float32?](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746)

<div class="topic-metadata">

**Author:** [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Replies:** 12\
**Last updated:** [November 27, 2023, 2:37pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746 "2023-11-27T14:37:39Z")

</div>

the following shows that s1 was hugely mis-computed. It happens when doing calculations in Float32. Everything is fine in Float64. Random.seed!(1); n = 1000000; a = rand(Float32, n); b = rand(Float32, 100, n); s1 = 0f0;…

---

## [Waveform relaxation and other parallel in time integrators](https://discourse.julialang.org/t/waveform-relaxation-and-other-parallel-in-time-integrators/106754)

<div class="topic-metadata">

**Author:** [@RahulManavalan](https://discourse.julialang.org/u/RahulManavalan)\
**Replies:** 2\
**Last updated:** [November 26, 2023, 9:10pm UTC](https://discourse.julialang.org/t/waveform-relaxation-and-other-parallel-in-time-integrators/106754 "2023-11-26T21:10:05Z")

</div>

I am looking into waveform relaxation for my project. Based on a cursory search, it seems to me that there is no registered package specifically for waveform relaxation. I did find an old issue Parallel ODE Solvers · Iss…

---

## [Solving Stochastic Differential Equations in Different Time Windows](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807)

<div class="topic-metadata">

**Author:** [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)\
**Replies:** 19\
**Last updated:** [November 26, 2023, 5:59pm UTC](https://discourse.julialang.org/t/solving-stochastic-differential-equations-in-different-time-windows/14807 "2023-11-26T17:59:23Z")

</div>

I want to solve a stochastic differential equation in different time windows under Wiener process, e.g., instead of solving the system for a time interval of (0, 2) seconds, I want to solve the system for a time duration…

---

## [Solving coupled ODE systems](https://discourse.julialang.org/t/solving-coupled-ode-systems/106685)

<div class="topic-metadata">

**Author:** [@JunmingDuan](https://discourse.julialang.org/u/JunmingDuan)\
**Replies:** 2\
**Last updated:** [November 24, 2023, 5:24pm UTC](https://discourse.julialang.org/t/solving-coupled-ode-systems/106685 "2023-11-24T17:24:26Z")

</div>

I would like to solve two coupled ode systems. The first is integrated by any method as usual, while the second is integrated once every 10 steps and its solution is used to modify the parameters (also the rhs) in the fi…

---

## [Taking Quaternions Seriously](https://discourse.julialang.org/t/taking-quaternions-seriously/44834)

<div class="topic-metadata">

**Author:** [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Replies:** 69\
**Last updated:** [November 24, 2023, 4:00pm UTC](https://discourse.julialang.org/t/taking-quaternions-seriously/44834 "2023-11-24T16:00:16Z")

</div>

There are at least 4 separate implementations of quaternions in Julia. Some of them are \<:AbstractVector, some are \<:AbstractMatrix (:neutral\_face:), some are \<:Number, some follow sijk convention, some follow ijks conv…

---

## [Strategy for finding eigenvalues of an infinite banded matrix](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513)

<div class="topic-metadata">

**Author:** [@Mason](https://discourse.julialang.org/u/Mason)\
**Replies:** 9\
**Last updated:** [November 23, 2023, 10:42am UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513 "2023-11-23T10:42:54Z")

</div>

So I’ve got a (Hermitian) like the following, and I want to find its first few eigenvalues: \\begin{pmatrix} \\ddots & \\ddots &\\ddots\\\\ \\ddots & T\_{i,i} &T\_{i+1,i} &T\_{i+2,i} & 0 & 0 & 0 &\\\\ \\ddots &…

---

## [Gridap.jl: Help with a Coupled PDE](https://discourse.julialang.org/t/gridap-jl-help-with-a-coupled-pde/84204)

<div class="topic-metadata">

**Author:** [@James\_Gorman](https://discourse.julialang.org/u/James_Gorman)\
**Replies:** 8\
**Last updated:** [November 23, 2023, 3:26am UTC](https://discourse.julialang.org/t/gridap-jl-help-with-a-coupled-pde/84204 "2023-11-23T03:26:26Z")

</div>

Looking over the various FEM PDE solvers (those written in julia and other languages), few of them have tutorials that deal with fracture mechanics. While there are different formulations that could be used, I want to im…

---

## [ReverseDiff + SciMLSensitivity method ambiguity](https://discourse.julialang.org/t/reversediff-scimlsensitivity-method-ambiguity/105718)

<div class="topic-metadata">

**Author:** [@Larbino1](https://discourse.julialang.org/u/Larbino1)\
**Replies:** 7\
**Last updated:** [November 18, 2023, 5:47pm UTC](https://discourse.julialang.org/t/reversediff-scimlsensitivity-method-ambiguity/105718 "2023-11-18T17:47:38Z")

</div>

Hi, I’m getting a method ambiguity when trying to use ReverseDiff + SciMLSensitivity. ERROR: MethodError: kwcall( ::NamedTuple{(:save\_everystep,), Tuple{Bool}}, ::typeof(DiffEqBase.solve\_up), ::ODEProblem{SVector{…

---

## [Implement JFNK with NonlinearSolve.jl](https://discourse.julialang.org/t/implement-jfnk-with-nonlinearsolve-jl/106339)

<div class="topic-metadata">

**Author:** [@gjparker](https://discourse.julialang.org/u/gjparker)\
**Replies:** 2\
**Last updated:** [November 17, 2023, 2:20pm UTC](https://discourse.julialang.org/t/implement-jfnk-with-nonlinearsolve-jl/106339 "2023-11-17T14:20:02Z")

</div>

I am trying to implement a Jacobian-Free Newton Krylov (JFNK) solver using existing functions in NonlinearSolve.jl and related packages. These methods should only require the evaluation of the directional derivatives J\*…

---

## [Why my diffeq solver idea doesn't work?](https://discourse.julialang.org/t/why-my-diffeq-solver-idea-doesnt-work/106350)

<div class="topic-metadata">

**Author:** [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)\
**Replies:** 4\
**Last updated:** [November 17, 2023, 1:32pm UTC](https://discourse.julialang.org/t/why-my-diffeq-solver-idea-doesnt-work/106350 "2023-11-17T13:32:09Z")

</div>

I was testing a silly idea on how to solve a difeq. #f(x+Δx) ≈ a + bΔx + cΔx^2 + dΔx^3 #f'(x+Δx) ≈ b + 2cΔx + 3dΔx^2 #a = f(x) #b = f'(x) # cΔx^2 + dΔx^3 = f(x+Δx) - a - bΔx # 2cΔx + 3dΔx^2 = f'(x+Δx) - b # (c,d) = \[Δx…

---

## [Gridap.jl: Solving non-linear coupled PDEs](https://discourse.julialang.org/t/gridap-jl-solving-non-linear-coupled-pdes/106326)

<div class="topic-metadata">

**Author:** [@charvey](https://discourse.julialang.org/u/charvey)\
**Replies:** 0\
**Last updated:** [November 16, 2023, 4:35pm UTC](https://discourse.julialang.org/t/gridap-jl-solving-non-linear-coupled-pdes/106326 "2023-11-16T16:35:56Z")

</div>

We are having some difficulties solving two coupled multi-variate non-linear PDEs using Gridap. The various tutorials haven’t quite filled the gaps in our knowledge. We have two PDEs that we need to solve together. Thei…

---

## [CuArrays with KrylovKit](https://discourse.julialang.org/t/cuarrays-with-krylovkit/106116)

<div class="topic-metadata">

**Author:** [@Gattu\_Mytraya](https://discourse.julialang.org/u/Gattu_Mytraya)\
**Replies:** 8\
**Last updated:** [November 14, 2023, 10:21am UTC](https://discourse.julialang.org/t/cuarrays-with-krylovkit/106116 "2023-11-14T10:21:51Z")

</div>

What is the recommended way of diagonalizing linear operators on CuArrays using KrylovKit? I am targeting the first 100 eigenvalues of a large matrix - Operating directly on CuArrays leads to a very large allocation on t…

---

## [Memory allocations with Zygote](https://discourse.julialang.org/t/memory-allocations-with-zygote/106236)

<div class="topic-metadata">

**Author:** [@Nikos\_Gianniotis](https://discourse.julialang.org/u/Nikos_Gianniotis)\
**Replies:** 8\
**Last updated:** [November 15, 2023, 7:40pm UTC](https://discourse.julialang.org/t/memory-allocations-with-zygote/106236 "2023-11-15T19:40:46Z")

</div>

I have noticed that when I use Zygote to do automatic differentiation in my code (typically I use Zygote to obtain the gradient of my loss function and pass it to Opim.optimize to minimise the loss), I often get high num…

---

## [Wrong gradient from Zygote?](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204)

<div class="topic-metadata">

**Author:** [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)\
**Replies:** 13\
**Last updated:** [November 14, 2023, 11:39pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204 "2023-11-14T23:39:58Z")

</div>

I have the following example where I get different gradient values for forward and backward: using Random, ChainRulesCore, Zygote, Statistics function simulate(N, drift, sigma, rng,i) ChainRulesCore.@ignore\_derivat…

---

## [Why would A \\ b grow with O(n^2)?](https://discourse.julialang.org/t/why-would-a-b-grow-with-o-n-2/106083)

<div class="topic-metadata">

**Author:** [@tawheeler](https://discourse.julialang.org/u/tawheeler)\
**Replies:** 5\
**Last updated:** [November 12, 2023, 8:03am UTC](https://discourse.julialang.org/t/why-would-a-b-grow-with-o-n-2/106083 "2023-11-12T08:03:00Z")

</div>

Hello, I was curious how long it would take my laptop to solve Ax = b for square A as a function of the problem size. I tried to profile it using BenchmarkTools: ns = \[1,2,3,5,7,10,20,30,50,70,100,200,300,500,700,1000,…

---

## [AD-able cubature](https://discourse.julialang.org/t/ad-able-cubature/106036)

<div class="topic-metadata">

**Author:** [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Replies:** 13\
**Last updated:** [November 11, 2023, 6:07pm UTC](https://discourse.julialang.org/t/ad-able-cubature/106036 "2023-11-11T18:07:24Z")

</div>

I have functions defined on \[0,1\]^4 \\times \[0,\\infty) that I want to integrate numerically, ideally in a way that I can use with an automatic differentiation framework. I know very little about these functions ex ante ot…

---

## [OrdinaryDiffEq solve: "TypeError in setfield!" when using MVector](https://discourse.julialang.org/t/ordinarydiffeq-solve-typeerror-in-setfield-when-using-mvector/105960)

<div class="topic-metadata">

**Author:** [@laikq](https://discourse.julialang.org/u/laikq)\
**Replies:** 0\
**Last updated:** [November 8, 2023, 4:20pm UTC](https://discourse.julialang.org/t/ordinarydiffeq-solve-typeerror-in-setfield-when-using-mvector/105960 "2023-11-08T16:20:53Z")

</div>

I get an error TypeError: in setfield!, expected Tuple{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, Base.RefValue{Int64}}, got a value of type Tuple{LinearAlgebra.LU{Float64, StaticArraysCore.MMatrix{4, 4…

---

## [Smith Normal Form?](https://discourse.julialang.org/t/smith-normal-form/105477)

<div class="topic-metadata">

**Author:** [@BLI](https://discourse.julialang.org/u/BLI)\
**Replies:** 3\
**Last updated:** [November 8, 2023, 5:44pm UTC](https://discourse.julialang.org/t/smith-normal-form/105477 "2023-11-08T17:44:40Z")

</div>

Does anyone have knowledge of packages for Smith Normal Form (e.g., SmithNormalForm.jl)? I did a simple test of SmithNormalForm.jl: Apart from the example in the GitHub page,… function smith seems to work: rectangula…

---

## [Anyone developing anything for Julia regarding SPH?](https://discourse.julialang.org/t/anyone-developing-anything-for-julia-regarding-sph/23708)

<div class="topic-metadata">

**Author:** [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Replies:** 7\
**Last updated:** [November 4, 2023, 6:03pm UTC](https://discourse.julialang.org/t/anyone-developing-anything-for-julia-regarding-sph/23708 "2023-11-04T18:03:42Z")

</div>

Where SPH, stands for Smoothed Particle Hydrodynamics. Does anyone happen to have a special interest in the field or done something in Julia regarding this methodology? Kind regards

---

## [Missing precision in the solution of combinations of ODE and jump processes using DifferentialEquations.jl and JumpProcesses.jl](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370)

<div class="topic-metadata">

**Author:** [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Replies:** 11\
**Last updated:** [October 31, 2023, 8:58am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370 "2023-10-31T08:58:48Z")

</div>

Hey everyone, For a toy model, I tried to simulate an ODE together with a jump process, which jump rate at a certain time is dependent on the ODE at that time. (Motivated by the growth of a tumor and a death process, wh…

---

## [Block diagonal and sparse matrices](https://discourse.julialang.org/t/block-diagonal-and-sparse-matrices/105551)

<div class="topic-metadata">

**Author:** [@loisel](https://discourse.julialang.org/u/loisel)\
**Replies:** 13\
**Last updated:** [October 30, 2023, 5:07pm UTC](https://discourse.julialang.org/t/block-diagonal-and-sparse-matrices/105551 "2023-10-30T17:07:23Z")

</div>

I have a code that does something like this: Z = inv(ABCDE) Here, A is a diagonal matrix. B is block diagonal, its diagonal blocks are 2x2. C is block diagonal, its diagonal blocks are 4x4. Subsequent matrices have lar…

---

## [Unexpectedly Large Compile Time Using MethodOfLines.jl](https://discourse.julialang.org/t/unexpectedly-large-compile-time-using-methodoflines-jl/105212)

<div class="topic-metadata">

**Author:** [@aharbick](https://discourse.julialang.org/u/aharbick)\
**Replies:** 3\
**Last updated:** [October 28, 2023, 11:07am UTC](https://discourse.julialang.org/t/unexpectedly-large-compile-time-using-methodoflines-jl/105212 "2023-10-28T11:07:10Z")

</div>

So I am trying to use the MethodOfLines package to solve some equations which are somewhat similar to 2D diffusion equations. As a starting point, I have modified the heat equation example from the documentation to work …

---

## [ApproxFun.jl for Helmholtz equation on cylinder](https://discourse.julialang.org/t/approxfun-jl-for-helmholtz-equation-on-cylinder/14864)

<div class="topic-metadata">

**Author:** [@mitkoge](https://discourse.julialang.org/u/mitkoge)\
**Replies:** 3\
**Last updated:** [October 27, 2023, 4:05pm UTC](https://discourse.julialang.org/t/approxfun-jl-for-helmholtz-equation-on-cylinder/14864 "2023-10-27T16:05:30Z")

</div>

I am trying to figure out if ApproxFun can be used for solving Helmholtz equation on cylindrical domain. Naively i hoped for something like Bessel spaces readily available. May be i have to construct cylinder domain us…

---

## [QuadGK returns NaNs](https://discourse.julialang.org/t/quadgk-returns-nans/105308)

<div class="topic-metadata">

**Author:** [@aner-sanchez](https://discourse.julialang.org/u/aner-sanchez)\
**Replies:** 4\
**Last updated:** [October 24, 2023, 1:24pm UTC](https://discourse.julialang.org/t/quadgk-returns-nans/105308 "2023-10-24T13:24:46Z")

</div>

Hello everyone :slight\_smile: Would someone know why using QuadGK n = 538 \_, w, \_ = kronrod(n,1e-11,40); @show w returns a vector of NaN weights at n=538 and not before?

[Previous page](https://discourse.julialang.org/c/domain/numerics/20.md?page=10)

[Next page](https://discourse.julialang.org/c/domain/numerics/20.md?page=12)