# Numerics

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

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

**Page:** 4

---

## [Calling C Code with Automatic Differentiation in Julia](https://discourse.julialang.org/t/calling-c-code-with-automatic-differentiation-in-julia/129503)

<div class="topic-metadata">

**Author:** [@fgittins](https://discourse.julialang.org/u/fgittins)\
**Replies:** 6\
**Last updated:** [June 2, 2025, 7:39am UTC](https://discourse.julialang.org/t/calling-c-code-with-automatic-differentiation-in-julia/129503 "2025-06-02T07:39:48Z")

</div>

Hi all, I’m reaching out to see if anyone in the Julia community has experience or advice on calling C code from Julia in a way that still supports automatic differentiation (AD). My motivation I’m working on a project…

---

## [Nonlinear boundary conditions in ApproxFun](https://discourse.julialang.org/t/nonlinear-boundary-conditions-in-approxfun/129504)

<div class="topic-metadata">

**Author:** [@MingjiaYan](https://discourse.julialang.org/u/MingjiaYan)\
**Replies:** 0\
**Last updated:** [May 31, 2025, 3:43pm UTC](https://discourse.julialang.org/t/nonlinear-boundary-conditions-in-approxfun/129504 "2025-05-31T15:43:07Z")

</div>

Hi, I’m trying to solve for the steady state solution of a one-dimensional heat equation describing a nonlinear heat conduction in a finite domain \[0,1\] using ApproxFun.jl: T\_{xx} = 0 subject to a Dirichlet boundary c…

---

## [Transient Stokes - Mass Matrix DAE - Dirichlet Boundary Conditions - Jacobian](https://discourse.julialang.org/t/transient-stokes-mass-matrix-dae-dirichlet-boundary-conditions-jacobian/129459)

<div class="topic-metadata">

**Author:** [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Replies:** 1\
**Last updated:** [May 30, 2025, 1:24pm UTC](https://discourse.julialang.org/t/transient-stokes-mass-matrix-dae-dirichlet-boundary-conditions-jacobian/129459 "2025-05-30T13:24:44Z")

</div>

I am attempting to solve transient Stokes flow in a canonical channel (parabolic inlet, no-slip bottom and top wall, fixed pressure and inlet). The Stokes problem leads to a mass-matrix DAE. A naive hand code backward E…

---

## [Finding Parameters for Boundary-Value Problem using DifferentialEquations.jl](https://discourse.julialang.org/t/finding-parameters-for-boundary-value-problem-using-differentialequations-jl/129386)

<div class="topic-metadata">

**Author:** [@parf](https://discourse.julialang.org/u/parf)\
**Replies:** 2\
**Last updated:** [May 28, 2025, 9:50am UTC](https://discourse.julialang.org/t/finding-parameters-for-boundary-value-problem-using-differentialequations-jl/129386 "2025-05-28T09:50:21Z")

</div>

Hello everyone! I have a system of equations, the solution of which defines a trajectory on a plane: dx/dt = Vx(x,y,t,p), dy/dt = Vy(x,y,t,p), where p is some unknown parameter. I know the starting point of the trajec…

---

## [UpperTriangular arrays in RecursiveArrayTools](https://discourse.julialang.org/t/uppertriangular-arrays-in-recursivearraytools/129352)

<div class="topic-metadata">

**Author:** [@Menno](https://discourse.julialang.org/u/Menno)\
**Replies:** 1\
**Last updated:** [May 26, 2025, 6:32pm UTC](https://discourse.julialang.org/t/uppertriangular-arrays-in-recursivearraytools/129352 "2025-05-26T18:32:15Z")

</div>

Hi all, In the (otherwise absolutely wonderful) RecursiveArrayTools package, I’m fairly sure it must be possible to use UpperTriangular or LowerTriangular arrays as parts of an ArrayPartition. The source code (Recursive…

---

## [Enzyme.jl Critical Issues with Basic Arithmetic Operations in Automatic Differentiation](https://discourse.julialang.org/t/enzyme-jl-critical-issues-with-basic-arithmetic-operations-in-automatic-differentiation/129216)

<div class="topic-metadata">

**Author:** [@mra](https://discourse.julialang.org/u/mra)\
**Replies:** 7\
**Last updated:** [May 26, 2025, 12:59pm UTC](https://discourse.julialang.org/t/enzyme-jl-critical-issues-with-basic-arithmetic-operations-in-automatic-differentiation/129216 "2025-05-26T12:59:59Z")

</div>

According to the Enzyme.jl documentation the “Enzyme differentiates arbitrary multivariate vector functions as the most general case in automatic differentiation”, however I encountered critical errors (including the ker…

---

## [ODE model with preallocated MVectors and still compatible with AD/stiff solvers?](https://discourse.julialang.org/t/ode-model-with-preallocated-mvectors-and-still-compatible-with-ad-stiff-solvers/129253)

<div class="topic-metadata">

**Author:** [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)\
**Replies:** 2\
**Last updated:** [May 22, 2025, 8:48pm UTC](https://discourse.julialang.org/t/ode-model-with-preallocated-mvectors-and-still-compatible-with-ad-stiff-solvers/129253 "2025-05-22T20:48:50Z")

</div>

I am writing an ODE model that describes the positions of a system of N interacting particles living in a 2D box. My ODE state vector u is an N x 2 array of numbers. Each call to the model updates a pre-allocated N x 2 a…

---

## [Please explain LinearAlgebra.\_\_normalize!](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636)

<div class="topic-metadata">

**Author:** [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Replies:** 7\
**Last updated:** [May 21, 2025, 3:10pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636 "2025-05-21T15:10:47Z")

</div>

I am trying to understand LinearAlgebra.\_\_normalize!, quoted in full here: @inline function \_\_normalize!(a::AbstractArray, nrm) # The largest positive floating point number whose inverse is less than infinity δ …

---

## [Multi-argument Jacobian and gradient](https://discourse.julialang.org/t/multi-argument-jacobian-and-gradient/129029)

<div class="topic-metadata">

**Author:** [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Replies:** 7\
**Last updated:** [May 21, 2025, 12:04pm UTC](https://discourse.julialang.org/t/multi-argument-jacobian-and-gradient/129029 "2025-05-21T12:04:31Z")

</div>

I have a function that takes multiple arguments (all of them vectors), and I want to obtain the Jacobian in each parameter (or, similarly, the gradient if the function maps to \\mathbb{R}). Is it best to unpack/pack like…

---

## [Eigen fails to find eigenvectors](https://discourse.julialang.org/t/eigen-fails-to-find-eigenvectors/129072)

<div class="topic-metadata">

**Author:** [@martinkjlarsson](https://discourse.julialang.org/u/martinkjlarsson)\
**Replies:** 28\
**Last updated:** [May 20, 2025, 3:14pm UTC](https://discourse.julialang.org/t/eigen-fails-to-find-eigenvectors/129072 "2025-05-20T15:14:41Z")

</div>

I am confused. I’m finding the eigenvalues and eigenvectors of a 3x3 matrix, but it seems like eigen does not find all the vectors correctly. julia\> using LinearAlgebra julia\> A = \[1 1 1; 1 1 1; -1 -1 -1\] 3×3 Matrix{In…

---

## [Linear interpolation given a triangulation in arbitrary dimensions](https://discourse.julialang.org/t/linear-interpolation-given-a-triangulation-in-arbitrary-dimensions/128998)

<div class="topic-metadata">

**Author:** [@Pablo\_Marchant](https://discourse.julialang.org/u/Pablo_Marchant)\
**Replies:** 9\
**Last updated:** [May 16, 2025, 2:15pm UTC](https://discourse.julialang.org/t/linear-interpolation-given-a-triangulation-in-arbitrary-dimensions/128998 "2025-05-16T14:15:05Z")

</div>

Hi all, I’m in need of the equivalent of LinearNDInterpolator in scipy for julia, in order to do interpolation of irregularly sampled points (not on a grid). There are existing packages such as NaturalNeighbours.jl, but…

---

## [How to think about the cost of vJ and Jv (pullback and pushforward)](https://discourse.julialang.org/t/how-to-think-about-the-cost-of-vj-and-jv-pullback-and-pushforward/129021)

<div class="topic-metadata">

**Author:** [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Replies:** 6\
**Last updated:** [May 15, 2025, 8:58am UTC](https://discourse.julialang.org/t/how-to-think-about-the-cost-of-vj-and-jv-pullback-and-pushforward/129021 "2025-05-15T08:58:49Z")

</div>

In a calculation involving an f: R^n \\to R^m function, where n ranges from 50 to a few hundred and m \\ge n with a similar scale, I have the option of either calculating the Jacobian J once and saving it, or calcula…

---

## [Reusing NoiseProcess in StochasticDiffEq](https://discourse.julialang.org/t/reusing-noiseprocess-in-stochasticdiffeq/128947)

<div class="topic-metadata">

**Author:** [@ymuv1](https://discourse.julialang.org/u/ymuv1)\
**Replies:** 2\
**Last updated:** [May 14, 2025, 9:33am UTC](https://discourse.julialang.org/t/reusing-noiseprocess-in-stochasticdiffeq/128947 "2025-05-14T09:33:13Z")

</div>

I am trying to use the same NoiseProcess in two different SDE integrations. I don’t mean to just use the same type of noise, but the same values as well (and potentially values in intermediate times obtained via the brid…

---

## [Optimizing the computation of jacobians for multiple shooting implementation](https://discourse.julialang.org/t/optimizing-the-computation-of-jacobians-for-multiple-shooting-implementation/128976)

<div class="topic-metadata">

**Author:** [@jackyarndley](https://discourse.julialang.org/u/jackyarndley)\
**Replies:** 10\
**Last updated:** [May 13, 2025, 9:38pm UTC](https://discourse.julialang.org/t/optimizing-the-computation-of-jacobians-for-multiple-shooting-implementation/128976 "2025-05-13T21:38:36Z")

</div>

Hi, I’ve been working with multiple shooting for trajectory optimization. I am currently using the sparse backend from DifferentiationInterface.jl and using DifferentialEquations.jl to perform the integration. A snippet …

---

## [Counter-intuitive difference in matrix-vector multiplication when the matrix is transposed](https://discourse.julialang.org/t/counter-intuitive-difference-in-matrix-vector-multiplication-when-the-matrix-is-transposed/127479)

<div class="topic-metadata">

**Author:** [@mrkn](https://discourse.julialang.org/u/mrkn)\
**Replies:** 4\
**Last updated:** [May 8, 2025, 6:41am UTC](https://discourse.julialang.org/t/counter-intuitive-difference-in-matrix-vector-multiplication-when-the-matrix-is-transposed/127479 "2025-05-08T06:41:11Z")

</div>

I found that matrix-vector multiplication in Julia uses a different algorithm when the matrix is wrapped by Transpose or Adjoint, compared to when it is a regular Matrix. As a result, the multiplication can sometimes pro…

---

## [Questions about Cholesky and QR factorizations](https://discourse.julialang.org/t/questions-about-cholesky-and-qr-factorizations/128728)

<div class="topic-metadata">

**Author:** [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Replies:** 1\
**Last updated:** [May 5, 2025, 3:58pm UTC](https://discourse.julialang.org/t/questions-about-cholesky-and-qr-factorizations/128728 "2025-05-05T15:58:14Z")

</div>

I’ve learned that. (But the title of the topic is very general, which is tempting for people to discuss things.) But I find that I cannot perform U \\ UQ (see the 🟠 line) due to ERROR: MethodError: no method matching /…

---

## [Slighly different results for conduction-advection with Neuman boundary conditions: Crank-Nicolson VS. Method of Lines](https://discourse.julialang.org/t/slighly-different-results-for-conduction-advection-with-neuman-boundary-conditions-crank-nicolson-vs-method-of-lines/128488)

<div class="topic-metadata">

**Author:** [@David\_F](https://discourse.julialang.org/u/David_F)\
**Replies:** 1\
**Last updated:** [May 4, 2025, 7:31am UTC](https://discourse.julialang.org/t/slighly-different-results-for-conduction-advection-with-neuman-boundary-conditions-crank-nicolson-vs-method-of-lines/128488 "2025-05-04T07:31:57Z")

</div>

INTRO: Heat conduction + advection problem… 1D spherical coordinates. For Dirichlet BC (bc\_left\_type = “Dirichlet” and bc\_right\_type = “Dirichlet”), the two codes (CN and MOL) return identical results up to machine prec…

---

## [Stokes Flow with Pressure Inlet Conditions](https://discourse.julialang.org/t/stokes-flow-with-pressure-inlet-conditions/128577)

<div class="topic-metadata">

**Author:** [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Replies:** 6\
**Last updated:** [May 1, 2025, 10:03am UTC](https://discourse.julialang.org/t/stokes-flow-with-pressure-inlet-conditions/128577 "2025-05-01T10:03:24Z")

</div>

Are there examples of simulating (Navier-)Stokes flow by imposing pressure inlet conditions using VoronoiFVM, WaterLilly, Oceananigans, Trixi or similar? Thx.

---

## [VoronoiFVM.BlockStrategy: Results documented?](https://discourse.julialang.org/t/voronoifvm-blockstrategy-results-documented/128576)

<div class="topic-metadata">

**Author:** [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Replies:** 5\
**Last updated:** [April 30, 2025, 6:46pm UTC](https://discourse.julialang.org/t/voronoifvm-blockstrategy-results-documented/128576 "2025-04-30T18:46:07Z")

</div>

Congrats to the VoroniFVM developers to making their BlockStrategy available! Have results obtained using this strategy been reported somewhere? Thx!

---

## [Mean‐Field Games “Crowd Control” Problem in Julia: HJB OK, but Struggling With Fokker‐Planck Positivity in Gridap](https://discourse.julialang.org/t/mean-field-games-crowd-control-problem-in-julia-hjb-ok-but-struggling-with-fokker-planck-positivity-in-gridap/128134)

<div class="topic-metadata">

**Author:** [@Ciccisio](https://discourse.julialang.org/u/Ciccisio)\
**Replies:** 2\
**Last updated:** [April 29, 2025, 10:21am UTC](https://discourse.julialang.org/t/mean-field-games-crowd-control-problem-in-julia-hjb-ok-but-struggling-with-fokker-planck-positivity-in-gridap/128134 "2025-04-29T10:21:30Z")

</div>

Hi everyone, I’m relatively new to Julia and I’m trying to solve a crowd control Mean Field Game (MFG) problem (inspired by Yves Achdou’s work " Mean Field Games for Modeling Crowd Motion") consisting of: A Hamilton–J…

---

## [Custom dot product for GMRES](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313)

<div class="topic-metadata">

**Author:** [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Replies:** 14\
**Last updated:** [April 28, 2025, 3:39am UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313 "2025-04-28T03:39:58Z")

</div>

Is it possible to define a custom inner product for gmres in Krylov.jl ?

---

## [Could \`svd!\` support option \`check=false\` like \`lu!\` and other factorizations?](https://discourse.julialang.org/t/could-svd-support-option-check-false-like-lu-and-other-factorizations/128386)

<div class="topic-metadata">

**Author:** [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Replies:** 5\
**Last updated:** [April 25, 2025, 12:17pm UTC](https://discourse.julialang.org/t/could-svd-support-option-check-false-like-lu-and-other-factorizations/128386 "2025-04-25T12:17:38Z")

</div>

The check=false option present in lu!, bunchkaufman! and other factorization routines is super useful to handle exceptions efficiently. Unfortunately, the svd! doesn’t have this option and I am getting LinearAlgebra.LAPA…

---

## [Error with jac\_prototype in Mass Matrix Form of ODEFunction (svd! on Sparse Matrix)](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513)

<div class="topic-metadata">

**Author:** [@sanjeev](https://discourse.julialang.org/u/sanjeev)\
**Replies:** 5\
**Last updated:** [April 18, 2025, 3:31am UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513 "2025-04-18T03:31:04Z")

</div>

Hi all, I’m solving a DAE system using DifferentialEquations.jl in the mass matrix form. The issue arises when I include jac\_prototype in the ODEFunction. If I define ODEFunction without jac\_prototype, everything works …

---

## [IncompleteLU.jl efficiency much worse than Matlab ilu](https://discourse.julialang.org/t/incompletelu-jl-efficiency-much-worse-than-matlab-ilu/106062)

<div class="topic-metadata">

**Author:** [@rli1996](https://discourse.julialang.org/u/rli1996)\
**Replies:** 3\
**Last updated:** [April 17, 2025, 6:01pm UTC](https://discourse.julialang.org/t/incompletelu-jl-efficiency-much-worse-than-matlab-ilu/106062 "2025-04-17T18:01:11Z")

</div>

Hello Julia user, I’m developing a model with large matrix and systems of linear equation solver. I plan to change the original matlab code in Julia as I heard Julia is much faster than MatLab. However, I found ilu in …

---

## [Automatic Unit Conversion Function](https://discourse.julialang.org/t/automatic-unit-conversion-function/128121)

<div class="topic-metadata">

**Author:** [@HSnowden](https://discourse.julialang.org/u/HSnowden)\
**Replies:** 0\
**Last updated:** [April 16, 2025, 10:57am UTC](https://discourse.julialang.org/t/automatic-unit-conversion-function/128121 "2025-04-16T10:57:33Z")

</div>

Hi Is it possible in Unitful.jl or DynamicQuantities.jl given a set of preferred units e.g. \[u"ms",u"mm",u"mJ",u"g"\] to have a function that can convert any given value e.g. 10u"W" to the correct units. I have tried upr…

---

## [Reusing random numbers in EnsembleProblem](https://discourse.julialang.org/t/reusing-random-numbers-in-ensembleproblem/127514)

<div class="topic-metadata">

**Author:** [@aleCombi](https://discourse.julialang.org/u/aleCombi)\
**Replies:** 6\
**Last updated:** [April 14, 2025, 2:13pm UTC](https://discourse.julialang.org/t/reusing-random-numbers-in-ensembleproblem/127514 "2025-04-14T14:13:46Z")

</div>

I am trying to take derivatives of Montecarlo expectations against parameters using DifferentialEquations.jl. I am doing it using geometric brownian motion, and I want to bump the sigma and reevaluate (later I will try A…

---

## [Debugging non-finite Jacobian elements (w/ ForwardDiff)](https://discourse.julialang.org/t/debugging-non-finite-jacobian-elements-w-forwarddiff/127448)

<div class="topic-metadata">

**Author:** [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Replies:** 4\
**Last updated:** [April 8, 2025, 12:02pm UTC](https://discourse.julialang.org/t/debugging-non-finite-jacobian-elements-w-forwarddiff/127448 "2025-04-08T12:02:04Z")

</div>

I have a rather complex calculation that maps a vector of floats to a longer vector of floats. At some inputs, ForwardDiff gives me a Jacobian where some elements are NaN. Inputs are of course finite, and so is the func…

---

## [Solving First Order Differential Equations with Julia](https://discourse.julialang.org/t/solving-first-order-differential-equations-with-julia/127485)

<div class="topic-metadata">

**Author:** [@ghoshr](https://discourse.julialang.org/u/ghoshr)\
**Replies:** 2\
**Last updated:** [March 29, 2025, 5:11pm UTC](https://discourse.julialang.org/t/solving-first-order-differential-equations-with-julia/127485 "2025-03-29T17:11:46Z")

</div>

Solving First Order Differential Equations with Julia I wrote a full tutorial on solving first order ODEs in Julia using the DifferentialEquations.jl package. It provides a gentle introduction to Differential Equations i…

---

## [Heston SDE exact methods](https://discourse.julialang.org/t/heston-sde-exact-methods/127196)

<div class="topic-metadata">

**Author:** [@aleCombi](https://discourse.julialang.org/u/aleCombi)\
**Replies:** 5\
**Last updated:** [March 28, 2025, 10:31am UTC](https://discourse.julialang.org/t/heston-sde-exact-methods/127196 "2025-03-28T10:31:04Z")

</div>

I am new to DifferentialEquations.jl and I tried out some SDE simulations with EM. I implemented an exact simulation of Heston SV model using the paper of Broadie and Kaya (https://www.researchgate.net/publication/22024…

---

## [Integrating over piecewise approxfun](https://discourse.julialang.org/t/integrating-over-piecewise-approxfun/127309)

<div class="topic-metadata">

**Author:** [@nilsbecker](https://discourse.julialang.org/u/nilsbecker)\
**Replies:** 8\
**Last updated:** [March 25, 2025, 4:08pm UTC](https://discourse.julialang.org/t/integrating-over-piecewise-approxfun/127309 "2025-03-25T16:08:26Z")

</div>

Following up on my earlier post, I now have an ApproxFun.jl Fun defined over a piecewise domain: using ApproxFun f\_pw = Fun(x -\> (x \< 1.0 ? exp(x) : 1.), union(Segment(0.,1.), Segment(1.,2.))) f\_pw\_i = integrate(f\_pw) T…

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

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