0

I've obtained a numerical solution of a differential equation in a form of a vector (i.e., M(170,1)) by using ode45 (MATLAB) and tspan=0:500 (time range). How can I now integrate $M(t)$ with respect to time when I am not sure how time sequences are distributed? I know, that I could put tspan=0:0.01:500, but I would like to save some time now.

horchler
  • 3,203
Boki
  • 35
  • How accurate does your integral need to be? Do you have any simple bounds on the derivative(s) of the solution from the differential equation? I would expect that you would get decent accuracy by just applying the trapezoidal rule with your numerical solution, but you may want more accuracy than that will provide. Or are you simply asking about how you could calculate the integral when you don't have the values of time which were used? ode45 can provide those; use [tout,yout]=ode45(...) – Ian Jun 11 '15 at 16:33

2 Answers2

1

Consider augmented ODE of the form $$ M'(t) = F(t, M(t))\\ J'(t) = M(t)\\ M(0) = M_0\\ J(0) = 0 $$ Then $J(T)$ would be the $$ J(T) = \int_0^T M(t) dt. $$

uranix
  • 7,503
0

Since you're using Matlab, you could use cumtrapz, which implements trapezoidal rule quadrature cumulatively:

Mint = cumtrapz(tout,M);

If you think you want something more, you could try combining cubic interpolation (or another interp1 option) with one of Matlab more advanced quadrature methods, e.g., integral:

Mint = zeros(length(tout),1);
for i = 2:length(tout)
    Mint(i) = integral(@(x)interp1(tspan,M,x,'cubic'),tout(1),tout(i))
    % Or maybe try:
    %Mint(i) = Mint(i-1)+integral(@(x)interp1(tout,M,x,'cubic'),tout(i-1),tout(i))
end

Both of these options will return the integral of your ode45 solution at each of the time points in the output time vector tout.

horchler
  • 3,203