error in ode solver for heaviside function

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');

[Tomorrow]

Has anyone any idea? thanks so much

Answer 1

[Star Strider]

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.

[Tomorrow]

First Really appreciate your help.

1. I don't understand your saying about the FAI1 and FAI2. they are two step functions in "anti-phase". sorry to miss your intention of the first script.

2. agreed with your idea for this discontinuable problem. I have tried your codes, which took me hours to finish, sadly, the results are not so reasonable. in the second segment, matlab reminds that

 Warning: Failure at t=1.101024e+03.  Unable to meet integration tolerances
 without reducing the step size below the smallest value allowed
(3.637979e-12) at time t.

which is just at the beginning of this segment, I don't know why ?

and for the third part, it also reaches the minimum time step adding that

 Error using horzcat
Dimensions of matrices being concatenated are not consistent.

3. the results of the first segment are either consistent with the values when I directly replace all heaviside with the coefficients, which means setting FAI1 to be 1 while FAI2 to be 0.

Again, thanks a lot for your help!

[Star Strider]

My pleasure.

1. The problem is that FAI1 and FAI2 are step functions, and essentially not differentiable. This causes problems for all numeric ODE solvers.

2. The ‘step’ introduces severe discontinuities. You may have to define the ‘tspan’ segments specifically to avoid having to differentiate across the discontinuities. Attempting to integrate across the discontinuities is throwing the Warning.

‘... which means setting FAI1 to be 1 while FAI2 to be 0.’

I suggest you do that, with appropriate limits on the tspan segments. (Those do not have to be in a cell array. I originally was using segments of a time vector rather than vectors of limits to define ‘tspan’.)

The concatenation error may result from some elements of ‘dx’ being empty. You will likely have to troubleshoot that, since I have no idea what you are doing in your code, and what liberties I can take with it.

[Tomorrow]

sorry for the delay to response.

I think the problems may come from the equations itself, or the coefficients. so I am checking this and I will update the results if possible.

Thanks for your inspiring answers.

[Star Strider]

We all have lives outside of MATLAB Answers. No problem with the delay.

I will help as I can.

As always, my pleasure!