Issue integrating inside ode45 loop
When I try to use this Matlab code it goes in an infinite loop. I am trying to perform integration inside ode45:
clear clc options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); [T,Y] = ode45(@rigid,[0 12],[0 1 1],options); plot(T,Y(:,1),'+',T,Y(:,2),'*',T,Y(:,3),'.') function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) ; dy(2) = -y(1) * y(3); fun = @(t) exp(-t.^2).*log(t).^2+y(1); q = integral(fun,0,Inf); dy(3) = y(2) * y(3) + q;
There is no "infinite loop." Your function just takes a very long time to integrate. Try setting tspan to [0 1e-7]. It appears to be a high frequency oscillation, but I don't know if your equations are correct (that's a math question rather than a programming one). Such systems are hard to integrate accurately (ode15 might be a better choice), let alone quickly.
You also didn't bother to mention the important fact that the call to integral is generating a warning message:
Warning: Minimum step size reached near x = 1.75484e+22. There may be a singularity, or the tolerances may be too tight for this problem. > In funfun/private/integralCalc>checkSpacing at 457 In funfun/private/integralCalc>iterateScalarValued at 320 In funfun/private/integralCalc>vadapt at 133 In funfun/private/integralCalc at 84 In integral at 88 In rtest1>rigid at 17 In ode15s at 580 In rtest1 at 5
Printing out warning messages on each iteration greatly slows down integration. There's a good reason for this warning. You do realize that the function that you're evaluating with integral from 0 to Inf is equivalent to the following, right?
sqrt(pi)*((eulergamma + log(4))^2/8 + pi^2/16) + Inf*y(1)
where eulergamma is psi(1) or double(sym('eulergamma')). Your integrand doesn't converge.
If you like, can try to avoid the warning message in one of two ways.
1. Turn off the the warning (being sure to re-enable it afterwards). You can do that with the following code:
... warning('OFF','MATLAB:integral:MinStepSize'); [T,Y] = ode45(@rigid,[0 12],[0 1 1],options); warning('ON','MATLAB:integral:MinStepSize'); ...
You can obtain the the ID for a waring via the lastwarn function.
2. The other option might be to change your integration bounds and avoid the warning altogether, e.g.:
... q = integral(fun,0,1e20); ...
This may or may not be acceptable, but integral is not be returning the correct solution either because the result doesn't converge.