I have a problem that when I use the Heaviside function in the For loop with ODE23 or ODE45, the code does not give exact results (always give 0)

Question

saidovic 2021-1-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/724183-i-have-a-problem-that-when-i-use-the-heaviside-function-in-the-for-loop-with-ode23-or-ode45-the-cod

评论： Walter Roberson 2021-1-23

function tesode23
t=[0 11];
initial_x=0;
initial_dxdt=0;
[t,x]=ode23s(@F,t,[initial_x initial_dxdt]);
plot(t,x(:,1));
xlabel('t'); ylabel('x');
figure 
plot(t,x(:,2));
xlabel('t'); ylabel('x');
function dx=F(t,x)
L=30;
d=[0 3 14 17 20.525 23.025 42.2125 45.0125 67.5125 70.0125 92.5125 95.0125 117.5125 120.0125 142.5125 145.0125 167.5125 170.0125 192.5125 195.0125     217.5125 220.0125 242.5125 245.0125 267.5125 270.0125 292.5125 295.0125 317.5125 ];
v=200*1000/3600;
T=d/v;
Tn=(L+d)/v; 
w=40;
O=b*v/L;
    for i=1:length(t)
         for N=1:length(T);
                 K1=t(i)-T(N);
                 FAI1=heaviside(t(i)-T(N))-heaviside(t(i)-Tn(N));
                dx=zeros(2,1);
                dx(1)=x(2);
                dx(2)=-w^2*x(1)-sum((cosh(O*K1)+sinh(O*K1)+cos(O*K1)+sin(O*K1))*FAI1);   
    end
       end
end
end

2 个评论
显示无隐藏无

Mischa Kim 2021-1-22

Hi Said, please share the differential equation you are trying to solve, e.g. as a screen shot.

saidovic 2021-1-23

编辑：saidovic 2021-1-23

the differential equation : x"-w^2*x'= sum((cosh(O*K1)+sinh(O*K1)+cos(O*K1)+sin(O*K1))*M)

with M=heaviside(t-tj)-heaviside(t-tn);

i just want to know if the code (heaviside function in For loop in ODE) is correct.

thank you

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2021-1-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/724183-i-have-a-problem-that-when-i-use-the-heaviside-function-in-the-for-loop-with-ode23-or-ode45-the-cod#answer_604362

在 MATLAB Online 中打开

The ode*() functions are all designed using mathematical theory that requires that the functions being used and their first two derivatives are continuous. When that condition is violated, if you are lucky the ode*() functions will explicitly tell you that it was not able to find a solution; if you are not fortunate then the ode*() function will produce an answer anyhow, but that answer will be incorrect (unless the discontinuity is very small.)

To write that another way: it is not valid to use heaviside() inside a function called by ode23s() and you had the bad fortunate that it gave you a number anyhow instead of giving you and error message.

Though to be more correct, it is not that it is always invalid to have heaviside() there, but in any one call to ode23s(), the heaviside() always has to give the same answer.

L=30;
d=[0 3 14 17 20.525 23.025 42.2125 45.0125 67.5125 70.0125 92.5125 95.0125 117.5125 120.0125 142.5125 145.0125 167.5125 170.0125 192.5125 195.0125     217.5125 220.0125 242.5125 245.0125 267.5125 270.0125 292.5125 295.0125 317.5125 ];
v=200*1000/3600;
T=d/v;
Tn=(L+d)/v; 

Those are constants; you could pre-calculate T and Tn before the function and pass them in to the function if you needed tohttp://www.mathworks.com/help/matlab/math/parameterizing-functions.html

heaviside(t-tj)-heaviside(t-tn)

The T, Tn pairs give you a series of intervals over which the difference of the heaviside() calls is non-zero. Instead of using ode options for event functions, you can loop using different timespan entries.

for K = 1 : length(Tn)
    tspan = [T(K), Tn(K)];
    [T{K,1}, outputs{K,1}] = ode23s(function_with_M_1, tspan, current_boundary);
    current_boundary = outputs{K,1}(end,:);
    if K < length(Tn)
        tspan = [Tn(K), T(K+1)];
        [T{K,2}, outputs{K,2}] = ode23s(function_with_M_0, tspan, current_boundary);
        current_boundary = outputs{K,2}(end,:);
    end
end

Here, function_with_M_1 is the function for the case where M is 1, and function_with_M_0 is the function for the case where M is 0.

8 个评论
显示 6更早的评论隐藏 6更早的评论

saidovic 2021-1-23

在 MATLAB Online 中打开

matlab say the code has errors

for K = 1 : length(Tn)
    tspan = [T(K), Tn(K)];
    [T{K,1}, outputs{K,1}] = ode23s(function_with_M_1, tspan, current_boundary);
    current_boundary = outputs{K,1}(end,:);
    if K < length(Tn)
        tspan = [Tn(K), T(K+1)];
        [T{K,2}, outputs{K,2}] = ode23s(function_with_M_0, tspan, current_boundary);
        current_boundary = outputs{K,2}(end,:);
    end
end

Walter Roberson 2021-1-23

在 MATLAB Online 中打开

L = 30;
d = [0 3 14 17 20.525 23.025 42.2125 45.0125 67.5125 70.0125 92.5125 95.0125 117.5125 120.0125 142.5125 145.0125 167.5125 170.0125 192.5125 195.0125     217.5125 220.0125 242.5125 245.0125 267.5125 270.0125 292.5125 295.0125 317.5125 ];
v = 200*1000/3600;
T = d/v;
Tn = (L+d)/v; 
P.w = 40;
P.O = b*v/L;
    
current_boundary = [initial_x initial_dxdt];
for K = 1 : length(Tn)
    tspan = [T(K), Tn(K)];
    [T{K,1}, outputs{K,1}] = ode23s(@(t,x)function_with_M(t,x,P,1), tspan, current_boundary);
    current_boundary = outputs{K,1}(end,:);
    if K < length(Tn)
        tspan = [Tn(K), T(K+1)];
        [T{K,2}, outputs{K,2}] = ode23s(@(t,x)function_with_M(t,x,P,0), tspan, current_boundary);
        current_boundary = outputs{K,2}(end,:);
    end
end
function dx = function_with_M_1(t, x, P, FAI1)
    dx = zeros(2,1);
    dx(1) = x(2);
    dx(2) = -P.w^2*x(1)-sum((cosh(P.O*K1)+sinh(P.O*K1)+cos(P.O*K1)+sin(P.O*K1))*FAI1);   
end