MATLAB: Error in ode solver for heaviside function

heavisideode'sstiff

how to solve a ODEs with heaviside function?
I am solving a ODEs with a heaviside funstion in it. since it will bring stiffness to the problems, so i use ODE15s and ODE23s. the former took a lot of time and the result is not what I want to some unknown reasons. and the latter can' t solve due to low accuracy inherently. below is my function, How can I solve this?
%function
function dx=Integrated_ode(t,x)
t1=1100; t2=1200;
FAI1=1+heaviside(t-t2)-heaviside(t-t1);
FAI2=heaviside(t-t1)-heaviside(t-t2);
dx=zeros(4,1); 
dx(1)=x(2); % x(1)=displacement;c
dx(2)=-2.1100e-04*cos(0.4438*t)-2*0.0045*x(2)-(1-1.0429)*x(1)-167.1087*x(1)^3-0.1606*x(3)+0.1014*x(4);
dx(3)=x(2)-7.4957*x(3)+(-2.4162e+03*FAI2)*x(4);
dx(4)=-x(2)+(-3.8247e+03*FAI2)*x(3)+(591.7674*FAI1)*x(4);
by the way, the ODEs without the heaviside can be solved by the ODE23s and ODE15s. so the equation is ok I think.
%
heaviside(sym(0));
oldparam=sympref('HeavisideAtOrigin', 0);
[T,z]=ode15s(@Integrated_ode,[0 3000],[0 0 0 0]);
plot3(T,z(:,1),z(:,2),'b');

Best Answer

Numeric ODE solvers do not handle discontinuities well, so it is necessary to integrate it for each side of the discontinuities, using the previous ‘end’ results of the integration for the initial conditions for the subsequent integration.
Also, ‘FAI1’ and ‘FAI2’ seem to me to conflict, creating problems for the ODE solver.
Consider:
t = 0:3000;
t1=1100; t2=1200;
FAI1=1+heaviside(t-t2)-heaviside(t-t1);
FAI2=heaviside(t-t1)-heaviside(t-t2);
figure(1)
subplot(3,1,1)
plot(t, FAI1)
axis([0  3000    -0.1  1.1])
subplot(3,1,2)
plot(t, FAI2)
axis([0  3000    -0.1  1.1])
subplot(3,1,3)
plot(t, FAI1+FAI2)
axis([0  3000    -0.1  1.1])
You can deal with the discontinuity problem by breaking up the integration to solve it for ‘before’ and ‘after’ each discontinuity:
t1=1100; t2=1200;
td = {[0 t1-1], [t1+1 t2-1], [t2+1 3000]};
ic = [0 0 0 0];
for k1 = 1:length(td)
    Iteration = [k1 td{k1} ic]
    tspan = td{k1};
    [T,z]=ode15s(@Integrated_ode, tspan, ic);
    Tv{k1} = T;
    zv{k1} = z;
    ic = z(end,:);
end
figure(1)
plot3([Tv{:}],[zv{:}(:,1)],[zv{:}(:,2)],'b');
grid on
NOTE — I tested this partially. The integration for the centre segment takes forever, even with ode15s, so I leave you to test it beyond the first segment.
Related Question