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)

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

[Mischa Kim]

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

[saidovic]

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]

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.

[saidovic]

thank you very much sir, i do what you told but the program does not work.

note: I have two heaviside : M=heaviside (t-T) - heaviside (t-Tn),

that is to say : M=1 & =0 & =-1.

for that : <<Those are constants; you could pre-calculate (T=d/v) and (Tn=(L+d))/v) before the function and pass them in to the function>> i can't because i change " v " evry time .

thank you very much sir

[Walter Roberson]

You do not change v every time. You have

v=200*1000/3600;

which is independent of all inputs.

[Walter Roberson]

Are you sure you get -1 ?

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; 
T

T = 1×29

0 0.0540 0.2520 0.3060 0.3694 0.4144 0.7598 0.8102 1.2152 1.2602 1.6652 1.7102 2.1152 2.1602 2.5652 2.6102 3.0152 3.0602 3.4652 3.5102 3.9152 3.9602 4.3652 4.4102 4.8152 4.8602 5.2652 5.3102 5.7152

Tn

Tn = 1×29

0.5400 0.5940 0.7920 0.8460 0.9094 0.9544 1.2998 1.3502 1.7552 1.8002 2.2052 2.2502 2.6552 2.7002 3.1052 3.1502 3.5552 3.6002 4.0052 4.0502 4.4552 4.5002 4.9052 4.9502 5.3552 5.4002 5.8052 5.8502 6.2552

syms t
fplot(heaviside(t-T(1))-heaviside(t-Tn(1)), [0 1])

fplot(heaviside(t-T(2))-heaviside(t-Tn(2)), [0 1])

[saidovic]

i mean, i change v after evry result .

anyway v is not a problem i can calculate T and Tn but the problem is in the loop and in the form of the code.

thank you for your time M.Walter

[Walter Roberson]

I do not know what you mean by "after every result" in this situation. Do you mean that you edit function F between runs? If that is what you mean, then for any one run T and Tn are constant and can be calculated before the run.

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; 

[t,x]=ode23s(@(t,x)F(t,x,T,Tn,L), t, [initial_x initial_dxdt]);

function dx=F(t, x, T, Tn, L)
w=40;
O=b*v/L;
%etc

[saidovic]

Oui, tu as raison, heaviside donne 0 et 1

but the code doesn't work, it says I have errors in the form

anyway thank you so much Sir

[saidovic]

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]

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

However this needs to be extended to handle K1. You show K1 being used in

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

but you do not define it.

Your code used t(i)-Tn(N) but t is scalar (it is the input time) and you were looping over all Tn which does not make sense in context.