# How to Calculate and Plot Fit Residuals from Numerical Integration to High Precision

**URL:** <https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094>\
**Category:** General Usage\
**Created:** [June 9, 2020, 8:38pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094 "2020-06-09T20:38:39Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 9, 2020, 8:38pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/1 "2020-06-09T20:38:39Z")

</div>

Hi, I am learning how to use Numerical Integration, and want to test how well the Num Integration fits the analytical/theoretical model, by subtracting the model from the Num Integration results, and plotting & examining those residuals. But I think the Num Integration output is some complicated interpolation object, not just a simple array of values. I’m going to eventually be doing very high precision Num integrations, any suggestions on how to best preserve this high precision in getting the residuals? Thanks for any info!

An example calculation:

```julia
function Spring1D(ddu,du,u,p,t)
    ddu[1] = -9*u[1]
end
x0 = [5.0]
v0 = [0.0]
tspan = (1.0,10.0)
probSpring1D = SecondOrderODEProblem(Spring1D,v0,x0,tspan)
solSpring1D = solve(probSpring1D,dense=true)

```

(Theoretical model here is 5\*Cos(3t). Any coding suggestions for the best way to calculate & plot a set of high-precision residuals?)

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [June 9, 2020, 8:49pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/2 "2020-06-09T20:49:47Z")

</div>

You can use `BigFloat`s to increase the precision and also use `setprecision()` to tweak it manually. To plot the results, you can use one of the many plotting libraries out there (PGFPlotsX.jl, Plots.jl, Gaston.jl, Makie.jl, UnicodePlots.jl, …).

I am not sure where you struggle to be honest, given that the code snippet you provide is neither working, nor is it self-descriptive. Anyways, I think both of your questions are answered now 😉

Btw. you can quote your code to make it more readable (here is a quick guide how to make it easier to help you: [Please read: make it easier to help you](https://discourse.julialang.org/t/psa-make-it-easier-to-help-you/14757))

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [June 9, 2020, 8:58pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/3 "2020-06-09T20:58:12Z")

</div>

If I understand your question correctly, you’re trying to access the raw solution values rather than interpolating. You seem to be using DifferentialEquations.jl, which describes how to access those values on the [Solution Handling](https://docs.sciml.ai/dev/basics/solution/#Accessing-the-Values-1) page:

> The solution type has a lot of built in functionality to help analysis. For example, it has an array interface for accessing the values. Internally, the solution type has two important fields:
> 
> 1. `u` which holds the Vector of values at each timestep
> 2. `t` which holds the times of each timestep.

…so you’d use something like

```julia
u_DiffEq = solSpring1D.u
u_analytical = your_theoretical_model(solSpring1D.t, params...)

```

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 9, 2020, 9:04pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/4 "2020-06-09T21:04:12Z")

</div>

Actually tamasgal, none of my questions seem to be answered above. (Although BigFloats and setprecision should be useful.) First, the snippet of code I gave does actually work fine. And, I am still asking how to write the code to define the function and subtract the residuals.

`Println(solSpring1D.u)` gives output like this:

`ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}[[0.0][5.0], [-2.22133e-8][5.0], . . . `

How to define the theoretical function `F(t) = 5 Cos(3t)`, and subtract it from the numerical results contained in the above mess? The function is continuous, the calculated array is some kind of interpolation object. What code would be be best to calculate and plot the residuals?  
Thanks…

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [June 9, 2020, 9:09pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/5 "2020-06-09T21:09:33Z")

</div>

Yeah, now that I see you are `using DifferentialEquations`, the code will work (I have not tried, but I trust you). I have not used that package yet but @stillyslalom’s hint looks very useful.

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 9, 2020, 9:09pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/6 "2020-06-09T21:09:51Z")

</div>

Ok stillyslalom, I’ll reads through the Solution Handling page. I know that solSpring.u has the f(t)and f’(t) data in it, just not sure how to extract it properly yet. Thanks.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [June 9, 2020, 9:28pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/7 "2020-06-09T21:28:06Z")

</div>

The reason you’re getting a complicated-looking `u` is because each timestep’s solution contains both `u'` and `u`. For example, `solSpring1D.u[1] = ([0.0], [5.0])`, which matches your initial conditions. You can use broadcasted `getindex` to access the individual components:

```julia
julia> v_DiffEq = getindex.(solSpring1D.u, 1)
52-element Array{Float64,1}:
   0.0
  -2.0478248738186e-8
  -8.618467381075347e-8
   ...

julia> x_DiffEq = getindex.(solSpring1D.u, 2)
52-element Array{Float64,1}:
  5.0
  5.0
  5.0
  4.999999999999994
  ...

```

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [June 9, 2020, 9:47pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/8 "2020-06-09T21:47:19Z")

</div>

A full MWE:

```julia
using DifferentialEquations, Plots

function Spring1D(ddu,du,u,p,t)
  ddu[1] = -9*u[1]
end
x0 = [5.0]
v0 = [0.0]
tspan = (0.0,10.0)
probSpring1D = SecondOrderODEProblem(Spring1D,v0,x0,tspan)
solSpring1D = solve(probSpring1D,dense=true)

u_DiffEq = getindex.(solSpring1D.u, 2)
u_analytical = @. 5cos(3solSpring1D.t)
err = abs.(u_analytical - u_DiffEq)
@. err[err == 0] = NaN # set zeros to NaN to avoid problems caused by log(0) = -Inf 

plot(solSpring1D.t, err, legend=false, 
     xlabel = "time", ylabel="Absolute error", yaxis=:log10)

```

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 9, 2020, 11:00pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/9 "2020-06-09T23:00:15Z")

</div>

Thanks stillyslalom, that works to get a discrete array & plot of residuals. I’m working through the syntax to figure it out (particularly the @ command).

Since the solution, `solSpring1D(time)[2]` , is interpolatable as a continuous function – and the analytical solution is obviously continuous – I was actually thinking of getting & plotting a “continuous” residuals error plot. But I guess it will look continuous enough if I make dt very small, with many more steps.

Thanks!

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [June 9, 2020, 11:44pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/10 "2020-06-09T23:44:25Z")

</div>

The `@.` command is a macro that converts all the subsequent commands to broadcasted form:

```plaintext
help?> @.
  @. expr

  Convert every function call or operator in expr into a "dot call" (e.g.
  convert f(x) to f.(x)), and convert every assignment in expr to a "dot
  assignment" (e.g. convert += to .+=).

```

You could equivalently write that line as `err[err .== 0] .= NaN`, but I prefer to use the dot macro for conditional replacement, since the conditions can get pretty convoluted. You can read more about the `@.` macro in this blog post: [More Dots: Syntactic Loop Fusion in Julia](https://julialang.org/blog/2017/01/moredots/)

The ‘problem’ is that your example is too easy for DiffEq.jl to integrate, so it needs very few timesteps to find an acceptable solution. Artificially clamping the timestep won’t really give you a clearer picture of the solver’s accuracy, since you wouldn’t choose such a fine discretization in day-to-day use. I would instead use a (smooth) line plot for the interpolated solution’s residual, overlaid with a scatter plot of the raw residual. Continuing on from my previous code,

```julia
t_fine = LinRange(tspan..., 1000)
u_DiffEq_fine = [solSpring1D(t)[2] for t in t_fine]
u_analytical_fine = @. 5cos(3t_fine)
err_fine = abs.(u_analytical_fine - u_DiffEq_fine)
@. err_fine[err_fine == 0] = NaN

plot(t_fine, err_fine, label="Interpolated solution")
scatter!(solSpring1D.t, err, label="Raw solution", 
     xlabel = "time", ylabel="Absolute error", yaxis=:log10, legend=:bottomright)

```

![image](https://global.discourse-cdn.com/julialang/original/3X/6/2/624deda59228883688ad2148c7db933e3e781e80.png)

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 11, 2020, 4:12am UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/11 "2020-06-11T04:12:18Z")

</div>

That is very helpful, I’m learning a lot from this, thanks again!

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 11, 2020, 10:37pm UTC](https://discourse.julialang.org/t/how-to-calculate-and-plot-fit-residuals-from-numerical-integration-to-high-precision/41094/12 "2020-06-11T22:37:22Z")

</div>

I’m kind of curious about why it spends the first ~5-7 steps at pretty much the same time value. I suppose the Diff-Eq solver needs a few points to get started.