This problem is quite specific and difficult to explain, but I am quite stuck so I wanted to ask to see if anyone has gotten over this problem before.
I have solved a differential equation and am now analysing the solution. For different values of time, I get that my function u(x) is a periodic sin type graph. What I wish to find is the total length for which the function is greater than 0. For example, if I had a perfect sin function from 0 to 2pi, this would be pi. In my case, each function u(x) is defined on [-1,1]. I found a root solver but since u(x) is really just an array, this does not work.
I have not attached any code, as there is not much that I have written. I am looking more for ideas, as I am completely lost as to which direction I should take. I was considering interpolating a function through all the points but this seems cumbersome.
How did you solve it? Probably using DifferentialEquations.jl? In that case, you have interpolation already built-in.
However, this reads to me that you actually solved a partial differential equation. In that case you should specify more how you solved it. Your text reads like you just discretized space. That would explain why you don’t have DiffEq’s interpolation in x. In that case, I think I would just find the first negative value and maybe interpolate linearly between it and the previous value to approximate the root. Ideally your resolution is high enough that this is a good approximation. If it is not high enough, then you should think about more elaborate interpolation schemes and probably need to exploit some “known” features of your solution for a good result.
If you have a numerical solution to a DE all you can do is;
use analytic considerations regarding the structure of the solution in order to get bounds on the oscillation of the solution, or structure of the solution (parametric family of functions in best cases). Then fit the numeric solution to the model, get exact numbers from the model, and use initial set of considerations to get error bounds.
The simplest version of this, provided you dint want error bars, would be to use a linear or pw constant approximation for the function between the solution points. eg pseudo
delta_x * size(filter(x->x>0,u))
or you could dig into the solver for the exact approximation used.
You can use a continuous callback to have DifferentialEquations.jl tell you precisely (to high-order accuracy) when your solution crosses zero. That is, find the roots while you are solving the ODE, rather than as post-processing.
I agree with the others that some custom interpolation is likely what you need. One package that makes this quite straightforward is this
You could throw your x and u(x) arrays into the interpolation function(s) (linear, quadratic, …) and then use a root finder on the resulting interpolation.
Of course, the choice of the interpolation scheme will generally matter, but if your grid along the x direction is fine enough, it should not matter too much.