error in ode solver for heaviside function

Question

Tomorrow 2017-9-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/355588-error-in-ode-solver-for-heaviside-function

评论： Star Strider 2017-9-20

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Tomorrow 2017-9-10

Has anyone any idea? thanks so much

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2017-9-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/355588-error-in-ode-solver-for-heaviside-function#answer_280977

在 MATLAB Online 中打开

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Tomorrow 2017-9-11

编辑：Tomorrow 2017-9-11

在 MATLAB Online 中打开

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!

Tomorrow 2017-9-20

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 2017-9-20

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

I will help as I can.

As always, my pleasure!

请先登录，再进行评论。

error in ode solver for heaviside function

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

error in ode solver for heaviside function

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

WeChat

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论