I don't know why this error happens...

y = [0 0 0 0];
yp = [0 0 0 0];
tspan = [0,1];
y0_new = [5;0;0;0];
yp0_new = [0;0;5;0];

% y(1) = x , y(2) = y
% y(3) = yp(1) = vx
% y(4) =  yp(2) = vy
% yp(3) = ax
% yp(4) = ay

[T,Y2] = ode15i(@f,tspan,y0_new,yp0_new);

Error using odearguments
The last entry in tspan must be different from the first entry.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in solution>f (line 44)
[t, Y1]= ode45(M1,[0, t2],[a, b / sqrt((2 * k1 * a)^2 + 1)]);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15i (line 116)
odearguments(FcnHandlesUsed, false, solver_name, ode, tspan, y0, ...

plot(y(1), y(2));

function res = f(t2,y,yp)
    
    m2 = 60;
    F2 = 300;
    fs2 = 150;
    fk2 = 0;
    
    syms x1(t);
    
    m1 = 60; 
    k1 = 1;
    F1 = 300;
    fs1 = 150;
    fk1 = 0;
    
    s = ((log(sqrt(1 + (2 * k1 * x1)^2) + 2 * k1 * x1)) / (4 * k1)) + (sqrt(1 + (2 * k1 * x1)^2) * x1) / 2;
    
    Ds1 = diff(s, t);
    D2s1 = diff(s, t, 2);
    
    dnjstlafur1 = 0.5 * (((m1 * Ds1^2) / (((2*k1*x1)^2 + 1)^1.5 / (2 * k1))) + abs(((m1 * Ds1^2) / (((2*k1*x1)^2 + 1)^1.5 / (2 * k1))) - fs1) - fs1); 
    
    ode1 = m1 * D2s1 == sqrt(F1^2 - dnjstlafur1^2);
   
    [V1] = odeToVectorField(ode1);
    M1 = matlabFunction(V1, 'vars', {'t', 'Y'});
    a = 0;
    b = 0;
    [t, Y1]= ode45(M1,[0, t2],[a, b / sqrt((2 * k1 * a)^2 + 1)]);
    
    xx = Y1(:,1);
    xx1 = xx(end);
    yy1 = xx1^2;

    % y(1) = x , y(2) = y
    % y(3) = yp(1) = vx
    % y(4) =  yp(2) = vy
    % yp(3) = ax
    % yp(4) = ay
    
    res = zeros(4,1);
    res(1) = yp(1) - y(3);
    res(2) = yp(2) - y(4);
    r2 = (y(3)^2 + y(4)^2)^1.5 / (y(3)*yp(4) - yp(3)*y(4));
    v2 = sqrt(y(3)^2 + y(4)^2);
    dnjstlafur2 = 0.5 * (((m2 * v2^2) / r2) + abs(((m2 * v2^2) / r2) - fs2) - fs2); 
    a2 = sqrt(yp(3)^2 + yp(4)^2);
    res(3) = m2 * a2 - sqrt(F2^2 - dnjstlafur2^2);
    res(4) = (y(1) - xx1) * yp(2) - (y(2) - yy1) * yp(1);
end

[Star Strider]

The error gets thrown by the ode45 call in ‘f’ the first time ‘f’ executes because ‘t2’ is zero so the ‘tspan’ argument to it is [0 0]. However even adding a small increment to ‘t2’ so it is not zero fails to completely solve the problem because the integration just ‘hangs’, and even changing the solver to ode15s (although the system does not appear to be ‘stiff’) does not change that. (That aside, it would likely be more efficient to solve that system separately if nothing in the differential equation changes between iterations — it does not appear to — rather than doing the derivation inside ‘f’ at each iteration. Save the derived differential equation as a separate function, then call the function from within ‘f’.)

I am not posting this as an answer because it does not completely solve the problem.

[Student]

Thanks for your help!!!

[Star Strider]

My pleasure!

Answer 1

[Walter Roberson]

[T,Y2] = ode15i(@f,tspan,y0_new,yp0_new);

You have an outer ode15i call. That is executing function f each time.

function res = f(t2,y,yp)

Function f is receiving the "current" time as the variable t2

        [t, Y1]= ode45(M1,[0, t2],[a, b / sqrt((2 * k1 * a)^2 + 1)]);

t2 is used as the upper bound of the tspan for an ode45() call.

ode15i is finding f to be very unstable right from the beginning, so ode15i is decreasing the time step again and again trying to find a scale that the function is predictable "enough" over.

ode15i is finding f to be very unstable. It reduces the time right down to eps(0), 4.94065645841247e-324 -- the smallest representable positive time.

ode15s is calling f with that t2 value of 4.94065645841247e-324 . f is calling ode45() passing in a time span of [0 4.94065645841247e-324] .

In the argument processing for ode45, ode45 tries to figure out how large of a step to start with. There are a few calculations involved, but as a simplification they involve taking the final time of 4.94065645841247e-324 and multiplying by 0.1 . But 4.94065645841247e-324 is the smallest representable number, so 0.1 times it underflows to 0.

ode45 then cross-checks and sees that the step size it has calculated is 0, and generates the MATLAB:odearguments:MaxStepLEzero message.

To avoid getting the error message about the step size, use options on the ode15i call to set a minimum step size greater than 10*eps(0) . I would suggest to you that it is not useful to have a minimum step size smaller than realmin, about .

Doing this would avoid the obscure error about minimum step size, and would instead get you a more understandable message from ode15i that it detected a likely singularity at time 0.

[Student]

thanks!

Answer 2

[Student]

y = [0 0 0 0];
yp = [0 0 0 0];
tspan = [0,0.01];
%[y0_new,yp0_new] = decic(@f,t0,y0,[0 0 0 0],yp0,[0 0 0 0]);
y0_new = [5;0;0;0];
yp0_new = [0;0;5;0];

% y(1) = x , y(2) = y
% y(3) = yp(1) = vx
% y(4) =  yp(2) = vy
% yp(3) = ax
% yp(4) = ay

%options = odeset('MaxStep', 1e7);

[T,Y2] = ode15i(@f,tspan,y0_new,yp0_new);
%plot(y(1), y(2));
% y1 = 0.035*Y(:,1)/50;
% y2 = 0.035*Y(:,2)/50;
% plot(T,[y1 y2])
% %plot(T,[y3 y4])
% xlabel('t(s)')
% ylabel('x(m)')
% figure(2)
% y3 = 0.035*Y(:,3);
% y4 = 0.035*Y(:,4);
% plot(T,[y3 y4])

function res = f(t2,y,yp)
    
    t2
    m2 = 60;
    F2 = 300;
    fs2 = 150;
    fk2 = 0;
    
    if t2 == 0
        xx1 = 0;
        yy1 = 0;
    else
        syms x1(t);
          
        m1 = 60;
        k1 = 1;
        F1 = 300;
        fs1 = 150;
        fk1 = 0;
        
        s = ((log(sqrt(1 + (2 * k1 * x1)^2) + 2 * k1 * x1)) / (4 * k1)) + (sqrt(1 + (2 * k1 * x1)^2) * x1) / 2;
        
        Ds1 = diff(s, t);
        D2s1 = diff(s, t, 2);
        
        dnjstlafur1 = 0.5 * (((m1 * Ds1^2) / (((2*k1*x1)^2 + 1)^1.5 / (2 * k1))) + abs(((m1 * Ds1^2) / (((2*k1*x1)^2 + 1)^1.5 / (2 * k1))) - fs1) - fs1); 
        
        ode1 = m1 * D2s1 == sqrt(F1^2 - dnjstlafur1^2);
        
        [V1] = odeToVectorField(ode1);
        M1 = matlabFunction(V1, 'vars', {'t', 'Y'});
        a = 0;
        b = 0;
        [t, Y1]= ode45(M1,[0, t2],[a, b / sqrt((2 * k1 * a)^2 + 1)]);
        
        xx = Y1(:,1);
        xx1 = xx(end);
        yy1 = xx1^2;

    end
    % y(1) = x , y(2) = y
    % y(3) = yp(1) = vx
    % y(4) =  yp(2) = vy
    % yp(3) = ax
    % yp(4) = ay
    
    res = zeros(4,1);
    res(1) = yp(1) - y(3);
    res(2) = yp(2) - y(4);
    res(3) = m2 * (sqrt(yp(3)^2 + yp(4)^2)) - sqrt(F2^2 - (0.5 * (((m2 * (sqrt(y(3)^2 + y(4)^2))^2) / ((y(3)^2 + y(4)^2)^1.5 / (y(3)*yp(4) - yp(3)*y(4)))) + abs(((m2 * (sqrt(y(3)^2 + y(4)^2))^2) / ((y(3)^2 + y(4)^2)^1.5 / (y(3)*yp(4) - yp(3)*y(4)))) - fs2) - fs2))^2);
    res(4) = (y(1) - xx1) * yp(2) - (y(2) - yy1) * yp(1);
end

The question I asked last time couldn't be taken out of a function, so I solved it with an if statement.

And then I ran this code, but this program doesn't ends, and then this error came out.

'matlab option 'MaxStep' must be greater than 0.'

So I looked into why it was like this.

The time gets smaller!!

[Walter Roberson]

You should expect time to get smaller when you use any of the ode*() functions.

All of the ode*() functions are adaptive. At any one stage, they have a position and a current step size. They use the current position and step size to select several different locations "nearby" and calculate the values there (evaluate the ode function with those trial positions.)

They then make a prediction based upon those results, and evaluate the ode function at the location chosen for the prediction: if the predicted values and actual values are "close enough" then the ode*() routine declares the trial to be a success and records the results... and increases the step size.

If, though, the predicted value and the actual value are considered to not be "close enough" then the ode*() routine declares the trial to fail, does not record the results -- and decreases the step size.

If you have several successes in a row, the step size will keep increasing -- the ode is "well behaved", smooth and predictable. Well-behaved functions do not need to be evaluiated as often. If the function can be estimated well, ode15i can rush to the end, taking larger and larger steps. When a prediction fails, it might not necessarily fail by much: all that might be needed is a small "course correction" after which perhaps the step size can get large again. And this implies that even for fairly well-behaved functions, ode15i() is continually trying to operate on the edge of stability.

When a failure occurs, then it might be that you have encounted an area that is significantly different than what has gone before -- so it might be necessary for there to be several failures in a row, each one reducing down to a smaller step size.

With odes presented to functions such as ode15i() as "black boxes", the ode*() routines cannot do analysis of the functions to calculate mathematically what the "right" step size should be at any time. You are near a rocky ledge in the dark; if you take a heck of a lot of very small steps, you will be fine, but that might take a long long time; the larger the step size the higher the risk that you might be eaten by a grue.

So imagine a situation in which ode15i guesses originally that it can use a step size of (say) 10, perhaps because the tspan is [0 10000]. It evaluates several locations "near" time 10, decides that the hills are too weird to be able to risk taking a step that large. So it reduces the step size to (say) 5 and tries again. As far as the ode function is concerned, time has just gone backwards from (near) 10 to (near) 5. Maybe it decides that step of 5 is still too risky and reduces the step down to (say) 2; again as far as the function being evaluated is concerned, time has gone backwards again. And maybe it figures that step 2 is still too risky and reduces down to step 1. Then 1/2 ... and so on. In the case of a very steep function, even the initial step can end up being a time on the order of

And so it is: your function is badly behaved enough that even down at a step size of or so, it is still failing.

[Torsten]

@Student

Hint to make your code simpler:

Y = integral_{t'=t0}^{t'=t} y(t') dt'

is equivalent to

dY/dt = y(t), Y(t0) = 0

Thus the call to ode45 is superfluous - the integral functions for which you solve using ode45 can be included as unknowns in the call to ode15i.

[Student]

thank you for your help!!