# Accuracy of ModelingToolkit linearizations

**URL:** <https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902>\
**Category:** Modelling & Simulations\
**Tags:** question, package, modelingtoolkit\
**Created:** [October 12, 2023, 5:29pm UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902 "2023-10-12T17:29:51Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 12, 2023, 5:29pm UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902/1 "2023-10-12T17:29:51Z")

</div>

I have a non-linear model with about 16 equations. I linearize it with different inputs and outputs around different operating point. I use both the linearize() function and the linerize\_symbolic() function. This works already much better than when using Matlab/Simulink. But I often get quite different results depending on my approach.

Question: Is it possible to increase the accuracy of the linearization by using a higher precision than Float64? If yes, how?

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [October 13, 2023, 2:05am UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902/2 "2023-10-13T02:05:10Z")

</div>

I didn’t run any tests, but it seems to be that the linearization itself probably introduces the main approximation error, not the round-off effects.

Depending on your nonlinear equations, the linearization around x\_0 introduces an error in the right-hand side which grows at the scale of \Vert Df(x\_0)(x - x\_0) - f(x) \Vert \leq C \Vert x - x\_0 \Vert^2. Which can cause errors at the scale \frac{\varepsilon}{L} ( e^{L \Vert x - x\_0 \Vert} - 1) where \varepsilon is the error in the RHS and L Lipschitz constant of the system (Hairer, Norsett, Wanner; Solving ODE I, Thm 10.2). As a result, the round-off error while computing the linearisation (which is maybe at most \varepsilon \approx 10^{-10}) is probably less relevant than the linearization error itself  
(\varepsilon \approx \sup\_x \Vert Df(x\_0)(x - x\_0) - f(x) \Vert \approx C \Vert x - x\_0 \Vert^2 ).

Maybe you could try to get a quadratic approximation instead of a linear one?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [October 13, 2023, 5:18am UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902/3 "2023-10-13T05:18:04Z")

</div>

> [@ufechner7](#):
>
> Question: Is it possible to increase the accuracy of the linearization by using a higher precision than Float64? If yes, how?

The linearization uses ForwardDiff.jl, so it’s already accurate to about machine epsilon. Unless your model is exceptionally poorly scaled numerically, I’d guess that there is something else going on.

> [@ufechner7](#):
>
> But I often get quite different results depending on my approach.

What are the approaches and the results you are talking about, `linearize` vs `linearize_symbolic`? Do you get any warnings from the linearization function?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 13, 2023, 8:41am UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902/4 "2023-10-13T08:41:51Z")

</div>

> [@baggepinnen](#):
>
> Unless your model is exceptionally poorly scaled numerically

Is there a way to figure out how well or poor a model is scaled?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [October 13, 2023, 1:55pm UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902/5 "2023-10-13T13:55:43Z")

</div>

You haven’t shared any details about the system or what results are unexpected, so we’re kind-of in the dark here

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [October 13, 2023, 2:55pm UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902/6 "2023-10-13T14:55:21Z")

</div>

A way to check for scaling problems is to look at the gradient at the approximation point. If some element of the gradient is very close to zero, then there could be a scaling problem. Generally speaking, it is nice if the elements of the gradient are of roughly the same order of magnitude.

Different results at a given point, say x, are entirely expected when the approximations are done at different points, say, x0 and x1.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2023, 9:32am UTC](https://discourse.julialang.org/t/accuracy-of-modelingtoolkit-linearizations/104902/7 "2023-10-16T09:32:34Z")

</div>

Automatic differentiation and symbolic differentiation should give fairly good accuracy in “most” contexts. On mathematical models for computing Jacobians, it should be extremely close (people normally say “within floating point error”, but that’s not exactly true, but I digress). It may not be exactly true because different optimizations in the AD library, particularly SIMD, will change the summation order, but that should be like 2 ulp in most instances.

In order to know more, this will need an example.