You can use interp1 to create an interpolated lookup table function: fx = [0.5644 0.6473 0.7258 0.7999 0.8697 0.9353 0.9967 1.0540 1.1072 1.1564 … 1.2016 1.2429 1.2803 1.3138 1.3435 1.3695 1.3917 1.4102 1.4250 1.4362 … 1.4438 1.4477 1.4482 1.4450 1.4384 1.4283 1.4147 1.3977 1.3773 1.3535 … 1.3263 1.2957 1.2618 1.2246 1.1841 1.1403 1.0932 1.0429 0.9893 … Read more
Your RK4 function is taking fixed steps that are much smaller than those that ode45 is taking. What you’re really seeing is the error due to polynomial interpolation that is used to produce the points in between the true steps that ode45 takes. This is often referred to as “dense output” (see Hairer & Ostermann … Read more
Let me start with the code and then step through to explain it. InnerFunc = function(x) { x + 0.805 } InnerIntegral = function(y) { sapply(y, function(z) { integrate(InnerFunc, 15, z)$value }) } integrate(InnerIntegral , 15, 50) 16826.4 with absolute error < 1.9e-10 The first line is very easy. We just need the function f(x) … Read more
I’m not sure if you can do exactly what you want, but it is possible to do quite a lot with events. First, this sounds like some sort of numerical calculation of first passage time (a.k.a first hitting time). If these “particles” are stochastic, stop and don’t use ode45 but rather a a method appropriate … Read more
Inserting large discontinuities in your ODE in the way you suggest (and the way illustrated by @macduff) can lead to less precision and longer computation times (especially with ode45 – ode15s might be a better option or at least make sure that your absolute and relative tolerances are suitable). You’ve effectively produced a very stiff … Read more
The AUC is approximated pretty easily by looking at a lot of trapezium figures, each time bound between x_i, x_{i+1}, y{i+1} and y_i. Using the rollmean of the zoo package, you can do: library(zoo) x <- 1:10 y <- 3*x+25 id <- order(x) AUC <- sum(diff(x[id])*rollmean(y[id],2)) Make sure you order the x values, or your … Read more