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:
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?)
You can use BigFloats 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)
1 Like
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 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:
u which holds the Vector of values at each timestept which holds the times of each timestep.
…so you’d use something like
u_DiffEq = solSpring1D.u
u_analytical = your_theoretical_model(solSpring1D.t, params...)
1 Like
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…
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.
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.
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> 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
...
1 Like
A full MWE:
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)
1 Like
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!
The @. command is a macro that converts all the subsequent commands to broadcasted form:
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
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,
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)
That is very helpful, I’m learning a lot from this, thanks again!
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.