I don't know why this code doesn't stops :(

This code doesn't stop when tspan > 0.4. But I want to know more than 10 seconds. Is there any solution?

y0 = [0 0 0 0];
yp0 = [0 0 0 0];
t0 = 0;
tspan = [t0,0.386];
[y0_new,yp0_new] = decic(@f,t0,y0,[1 1 1 1],yp0,[0 0 0 0])

y0_new = 4×1

0 0 0 0

yp0_new = 4×1

0 0 172.0836 289.5405

[T,Y] = ode15i(@f,tspan,y0_new,yp0_new);
y1 = 0.035*Y(:,1);
y2 = 0.035*Y(:,2);
plot(T,[y1 y2])
xlabel('t(s)')
ylabel('x(m)')

%plot(T,[Y(:,1),Y(:,2)])

function res = f(t,y,yp)
  I = 0.000122144; %kgm^2
  m = 0.19744; %kg
  g = 9.80665; %m/s^2
  R = 0.035; %m
  s = 0.090757; %rad
  r = 0.011515; %m
  M = 0.029244; %kg
  k = 15.36;
  A = 0.011065331; %m^2
  res = zeros(4,1);
  res(1) = yp(1) - y(3);
  res(2) = yp(2) - y(4);
  res(3) = I*yp(3) - (-m * sqrt(g^2 + R^2 * yp(4)^2 - 2 * g * R * yp(4) * sin(s)) * r * sin(s + y(1) - atan(R * yp(4) * cos(s) / (g - R * yp(4) * sin(s)))) - k * R * A * (y(3) - y(4)));
  res(4) = 2 * M * R * yp(4) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * yp(3)));
end

Answer 1

[Torsten]

res(3) == 0 and res(4) == 0 are two nonlinear equations in the unknowns yp(3) and yp(4). ode15i seems to have problems following their root path.

You can try to solve

res(3) == 0
res(4) == 0

within the function f using the nonlinear solver "fsolve" and return

res(3) = yp(3) - result(1)
res(4) = yp(4) - result(2) 

to ode15i where "result" is the "fsolve" solution. But maybe there is no or a non-unique solution for these two unknowns in the sequel of the computation - you should test it separately for a set of y-values where ode15i gives up.

[Student]

Can you give me a whole code?? please

[Torsten]

Do you know how the functions should behave for greater values of t ? Can you sketch y and y' ?

[Student]

I can't do such a thing like that... Do you have any solutions to solve these two differetial equations?

[Sam Chak]

@Student, If you're unable to sketch something like that, could you at least describe the nature of the differential equation?

[Student]

these are the two differential equations and the value of constant is mentioned above

[Torsten]

Do you have any solutions to solve these two differetial equations?

Are you kidding ? What else do you think is the code provided ? I explained two ways of solving the problem (one with explicit code, the other explained in detail), and I mentionned the numerical difficulties you could encounter. The rest is trying, trying, trying ...

Octave gives solutions for greater values of t with the original code, but when printing "res", I don't know whether they are trustworthy. It seems that at some time instances, there are branch points where the solver has to decide which branch for yp(3) and yp(4) is the correct to follow (like y = x^2 at x = 0) - and at these branch points, the solution gets inexact. Another possibility is that the denominator in the atan(...) expression becomes zero. That's why I gave you the hint to check the solution for yp(3:4) for the y-vector at the critical t-value.

[Student]

how much time does Ocatave gives solutions?

[Torsten]

how much time does Ocatave gives solutions?

Test it. I stopped it at t = 1.25 s because it became slower and slower.

You should prescribe an explicit array of output times "tspan" - else the output will collide with your available RAM.

As I said: no guarantee that Octave's "ode15i" takes the correct branches. That's why I asked if you have an impression about the expected result.

Better you try to find the reason for the difficulties instead of trying new methods for solution.

function main
  y0 = [0 0 0 0];
  yp0 = [0 0 0 0];
  t0 = 0;
  tspan = [t0,1.25];
  [y0_new,yp0_new] = decic(@f,t0,y0,[1 1 1 1],yp0,[0 0 0 0])
  [T,Y] = ode15i(@f,tspan,y0_new,yp0_new)
  figure(1)
  plot(T,[Y(:,1),Y(:,2)])
  numel(T)

  %I = 0.000122144; %kgm^2
  %m = 0.19744; %kg
  %g = 9.80665; %m/s^2
  %R = 0.035; %m
  %s = 0.090757; %rad
  %r = -0.011515; %m
  %M = 0.029244; %kg
  %k = 15.36;
  %A = 0.011065331; %m^2
  %x0 = [yp0_new(3),yp0_new(4)];
  %for i = 1:numel(T)
  %  y = Y(i,:);
  %  fun = @(x)[I*x(1) - (-m * sqrt(g^2 + R^2 * x(2)^2 - 2 * g * R * x(2) * sin(s)) * r * sin(s + y(1) - atan(R * x(2) * cos(s) / (g - R * x(2) * sin(s)))) - k * R * A * (y(3) - y(4)))...
  %             2 * M * R * x(2) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * x(1)))];
  %  sol = fsolve(fun,x0);
  %  dY(i,1) = y(3);
  %  dY(i,2) = y(4);
  %  dY(i,3) = sol(1);
  %  dY(i,4) = sol(2);
  %  res(i,:) = fun(sol);
  %  x0 = sol;
  %end
  %figure(2)
  %plot(T,res)
end

function res = f(t,y,yp)
  persistent iter
  if isempty(iter)
    iter = 0;
  endif
  iter = iter + 1;
  if mod(iter,10000)==0
    iter = 0;
    t
  endif
  I = 0.000122144; %kgm^2
  m = 0.19744; %kg
  g = 9.80665; %m/s^2
  R = 0.035; %m
  s = 0.090757; %rad
  r = -0.011515; %m
  M = 0.029244; %kg
  k = 15.36;
  A = 0.011065331; %m^2
  res = zeros(4,1);
  res(1) = yp(1) - y(3);
  res(2) = yp(2) - y(4);
  res(3) = I*yp(3) - (-m * sqrt(g^2 + R^2 * yp(4)^2 - 2 * g * R * yp(4) * sin(s)) * r * sin(s + y(1) - atan(R * yp(4) * cos(s) / (g - R * yp(4) * sin(s)))) - k * R * A * (y(3) - y(4)));
  res(4) = 2 * M * R * yp(4) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * yp(3)));
  %fun = @(x)[I*x(1) - (-m * sqrt(g^2 + R^2 * x(2)^2 - 2 * g * R * x(2) * sin(s)) * r * sin(s + y(1) - atan(R * x(2) * cos(s) / (g - R * x(2) * sin(s)))) - k * R * A * (y(3) - y(4)))...
  %           2 * M * R * x(2) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * x(1)))];
  %sol = fsolve(fun,[yp(3),yp(4)]);
  %res(3) = yp(3) - sol(1);
  %res(4) = yp(4) - sol(2);
end

[Sam Chak]

@Student, do the differential equations describe the motion of a Bicycle or a Moon?

[Steven Lord]

If you want to see why the solver's progress slows down, you could add an OutputFcn to your ode15i call and plot each point as the solver progresses in time. I'm betting that one or more of the solutions takes off like a rocket. You can kind of see them already starting their takeoff around t = 0.386 in the plot in the question. See the ballode or orbitode examples for two uses of the OutputFcn for this type of visualization.

[Torsten]

@Steven Lord

In Octave, they grow, but regularily up to about 100 at t = 1.25. But the time step size becomes smaller and smaller - maybe because the nonlinear system for yp(3:4) becomes more and more difficult to solve.

[Torsten]

@Student

You can clearly see from this code that the reason for failure is the division by zero in the atan() term. Again the mathematical problem and not the solver is responsible for failure.

I replaced

atan(R * y(6) * cos(s) / (g - R * y(6) * sin(s)))

by

asin(R * y(6) * cos(s) / sqrt((g - R * y(6) * sin(s))^2 + (R * y(6) * cos(s))^2))

and it works a little better.

y0 = [0 0 0 0 -172 289];
t0 = 0;
s = @(x)f(t0,[y0(1:4),x(1),x(2)]);
h = @(x) subsref(s(x), struct('type', '()', 'subs', {{[5,6]}}));
sol = fsolve(h,y0(5:6),optimset('TolFun',1e-10,'TolX',1e-10))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

sol = 1×2

172.0836 289.5405

y0 = [y0(1:4),sol];
tspan = [t0,20.0];
M = eye(6);
M(5,5) = 0;
M(6,6) = 0;
[T,Y] = ode15s(@f,tspan,y0,odeset('RelTol',1e-8,'AbsTol',1e-8,'Mass',M));

ans = 0.3248

ans = 0.3248

ans = 0.1052

ans = 0.1052

ans = -0.1224

ans = 3.9154

ans = 0.1052

ans = 0.1052

ans = 0.1052

ans = 0.1052

ans = 0.1052

ans = 0.1052

ans = 0.1052

ans = -0.1224

ans = 4.8001

ans = 0.0378

ans = 0.0378

ans = -0.0303

ans = 4.7753

ans = 0.0378

ans = 0.0378

ans = 0.0378

ans = 0.0378

ans = 0.0378

ans = 0.0378

ans = 0.0378

ans = -0.0303

ans = 4.8093

ans = 0.0175

ans = 0.0175

ans = -0.0030

ans = 4.8020

ans = 0.0175

ans = 0.0175

ans = 0.0175

ans = 0.0175

ans = 0.0175

ans = 0.0175

ans = 0.0175

ans = -0.0030

ans = 4.8122

ans = 0.0113

ans = 0.0113

ans = 0.0052

ans = -9.3687e-04

ans = 4.8082

ans = 0.0052

ans = 0.0052

ans = 0.0052

ans = 0.0052

ans = 0.0052

ans = 0.0052

ans = 0.0052

ans = -9.3687e-04

ans = 4.8143

ans = 0.0034

ans = 0.0034

ans = 0.0015

ans = 0.0015

ans = -3.2303e-04

ans = 4.8131

ans = 0.0015

ans = 0.0015

ans = 0.0015

ans = 0.0015

ans = 0.0015

ans = 0.0015

ans = 0.0015

ans = -3.2303e-04

ans = 4.8149

ans = 9.6582e-04

ans = 9.6582e-04

ans = 4.1349e-04

ans = 4.1349e-04

ans = -1.3890e-04

ans = 4.8146

ans = 4.1349e-04

ans = 4.1349e-04

ans = 4.1349e-04

ans = 4.1349e-04

ans = 4.1349e-04

ans = 4.1349e-04

ans = 4.1334e-04

ans = -1.3890e-04

ans = 4.8151

ans = 2.4778e-04

ans = 2.4778e-04

ans = 8.2062e-05

ans = 8.2062e-05

ans = -8.3656e-05

ans = 4.8150

ans = 8.2062e-05

ans = 8.2062e-05

ans = 8.2062e-05

ans = 8.2062e-05

ans = 8.2062e-05

ans = 8.2062e-05

ans = 8.1916e-05

ans = -8.3656e-05

ans = 4.8152

ans = 3.2347e-05

ans = 3.2347e-05

ans = -1.7369e-05

ans = 4.8152

ans = 3.2347e-05

ans = 3.2347e-05

ans = 3.2347e-05

ans = 3.2347e-05

ans = 3.2347e-05

ans = 3.2347e-05

ans = 3.2493e-05

ans = -1.7369e-05

ans = 4.8152

ans = 1.7432e-05

ans = 1.7432e-05

ans = 2.5178e-06

ans = 2.5178e-06

ans = -1.2397e-05

ans = 4.8152

ans = 2.5178e-06

ans = 2.5178e-06

ans = 2.5178e-06

ans = 2.5178e-06

ans = 2.5178e-06

ans = 2.5178e-06

ans = 2.3717e-06

ans = -1.2397e-05

ans = 4.8152

ans = -1.9566e-06

ans = 4.8152

ans = 1.1755e-06

ans = 1.1755e-06

ans = -1.6686e-07

ans = 4.8152

ans = 1.1755e-06

ans = 1.1755e-06

ans = 1.1755e-06

ans = 1.1755e-06

ans = 1.1755e-06

ans = 1.1755e-06

ans = 1.0293e-06

ans = -1.6686e-07

ans = 4.8152

ans = 7.7277e-07

ans = 7.7277e-07

ans = 3.7007e-07

ans = 3.7007e-07

ans = -3.2625e-08

ans = 4.8152

ans = 3.7007e-07

ans = 3.7007e-07

ans = 3.7007e-07

ans = 3.7007e-07

ans = 3.7007e-07

ans = 3.7007e-07

ans = 5.1620e-07

ans = -3.2625e-08

ans = 4.8152

ans = 2.4926e-07

ans = 1.2845e-07

ans = 1.2845e-07

ans = 7.6438e-09

ans = 7.6434e-09

ans = -1.1317e-07

ans = 4.8152

ans = 7.6434e-09

ans = 7.6434e-09

ans = 7.6434e-09

ans = 7.6434e-09

ans = 7.6434e-09

ans = 7.6434e-09

ans = 1.5377e-07

ans = -1.1317e-07

ans = 4.8152

ans = -2.8599e-08

ans = 4.8152

ans = -3.2294e-09

ans = 4.8152

ans = 4.3815e-09

ans = 1.1197e-09

ans = -2.1421e-09

ans = 4.8152

ans = 1.1197e-09

ans = 1.1197e-09

ans = 1.1197e-09

ans = 1.1197e-09

ans = 1.1197e-09

ans = 1.1197e-09

ans = 1.4725e-07

ans = -2.1421e-09

ans = 4.8152

ans = 1.4115e-10

ans = -8.3741e-10

ans = 4.8152

ans = 1.4114e-10

ans = 1.4114e-10

ans = 1.4114e-10

ans = 1.4114e-10

ans = 1.4114e-10

ans = 1.4114e-10

ans = -1.4599e-07

ans = -8.3741e-10

ans = 1.4529e-07

ans = 8.7497e-07

ans = 3.0659e-06

ans = 3.0659e-06

ans = 8.7497e-07

ans = 8.7497e-07

ans = 8.7497e-07

ans = 8.7497e-07

ans = 8.7497e-07

ans = 8.7497e-07

ans = 8.6035e-07

ans = 3.0659e-06

ans = -2.7943e-09

ans = 1.3249e-06

ans = -2.1094e-09

ans = 9.9516e-07

ans = -1.9040e-09

ans = 9.0979e-07

ans = -1.8423e-09

ans = 8.8530e-07

ans = -1.8238e-09

ans = 8.7806e-07

ans = -1.8183e-09

ans = 8.7621e-07

ans = -1.8169e-09

Warning: Failure at t=3.853933e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.881784e-16) at time t.

y1 = 0.035*Y(:,1);
y2 = 0.035*Y(:,2);
figure(1)
plot(T,[y1 y2])
xlabel('t(s)')
ylabel('x(m)')

for i = 1:numel(T)
  dydt = f(T(i),Y(i,:));
  res5(i) = dydt(5);
  res6(i) = dydt(6);
end

ans = 0.3248

ans = 0.1052

ans = 0.0378

ans = 0.0175

ans = 0.0113

ans = 0.0052

ans = 0.0034

ans = 0.0015

ans = 9.6582e-04

ans = 4.1349e-04

ans = 2.4778e-04

ans = 8.2062e-05

ans = 3.2347e-05

ans = 1.7432e-05

ans = 2.5178e-06

ans = 1.1755e-06

ans = 7.7277e-07

ans = 3.7007e-07

ans = 2.4926e-07

ans = 1.2845e-07

ans = 7.6434e-09

ans = 4.3815e-09

ans = 1.1197e-09

ans = 1.4114e-10

ans = 1.4529e-07

ans = 8.7497e-07

figure(2)
plot(T,[res5.',res6.'])

%plot(T,[Y(:,1),Y(:,2)])

function dydt = f(t,y)
  I = 0.000122144; %kgm^2
  m = 0.19744; %kg
  g = 9.80665; %m/s^2
  R = 0.035; %m
  s = 0.090757; %rad
  r = 0.011515; %m
  M = 0.029244; %kg
  k = 15.36;
  A = 0.011065331; %m^2
  
  if t > 3.85e-1
    g - R * y(6) * sin(s)
  end
  
  dydt = zeros(6,1);
  dydt(1) = y(3);
  dydt(2) = y(4);
  dydt(3) = y(5);
  dydt(4) = y(6);  
  dydt(5) = I*y(5) - (-m * sqrt(g^2 + R^2 * y(6)^2 - 2 * g * R * y(6) * sin(s)) * r * sin(s + y(1) - atan(R * y(6) * cos(s) / (g - R * y(6) * sin(s)))) - k * R * A * (y(3) - y(4)));
  %dydt(5) = I*y(5) - (-m * sqrt(g^2 + R^2 * y(6)^2 - 2 * g * R * y(6) * sin(s)) * r * sin(s + y(1) - asin(R * y(6) * cos(s) / sqrt((g - R * y(6) * sin(s))^2 + (R * y(6) * cos(s))^2))) - k * R * A * (y(3) - y(4)));
  dydt(6) = 2 * M * R * y(6) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * y(5)));
end

[Torsten]

If you download

www.unige.ch/~hairer/prog/MatlabStiff.7z

and use it together with the driver file

warning off
y0 = [0 0 0 0 172 289];
t0 = 0;
s = @(x)OdeFcn(t0,[y0(1:4),x(1),x(2)]);
h = @(x) subsref(s(x), struct('type', '()', 'subs', {{[5,6]}}));
sol = fsolve(h,y0(5:6),optimset('TolFun',1e-10,'TolX',1e-10));
y0 = [y0(1:4),sol];
tspan = [t0,4.07];
OdeFcn = @OdeFcn;
MassFcn = @MassFcn;
options = rdpset('RelTol',1e-8,'AbsTol',1e-8,'MassFcn',MassFcn,'InitialStep',1e-10);
[T,Y] = radau(OdeFcn,tspan,y0,options);
figure(1)
plot(T,[Y(:,1),Y(:,2)])
for i = 1:numel(T)
  dydt = OdeFcn(T(i),Y(i,:));
  res5(i) = dydt(5);
  res6(i) = dydt(6);
end
figure(2)
plot(T,[res5.',res6.'])
[m5,idx5] = max(res5);
[m6,idx6] = max(res6);
m5
T(idx5)
m6
T(idx6)

function dydt = OdeFcn(t,y)
  persistent iter
  if isempty(iter)
    iter = 0;
  endif
  iter = iter + 1;
  if mod(iter,10000)==0
    iter = 0;
    t
  endif
  I = 0.000122144; %kgm^2
  m = 0.19744; %kg
  g = 9.80665; %m/s^2
  R = 0.035; %m
  s = 0.090757; %rad
  r = 0.011515; %m
  M = 0.029244; %kg
  k = 15.36;
  A = 0.011065331; %m^2

  %if t > 3.85e-1
  %  g - R * y(6) * sin(s)
  %end

  dydt = zeros(6,1);
  dydt(1) = y(3);
  dydt(2) = y(4);
  dydt(3) = y(5);
  dydt(4) = y(6);
  %dydt(5) = I*y(5) - (-m * sqrt(g^2 + R^2 * y(6)^2 - 2 * g * R * y(6) * sin(s)) * r * sin(s + y(1) - atan(R * y(6) * cos(s) / (g - R * y(6) * sin(s)))) - k * R * A * (y(3) - y(4)));
  dydt(5) = I*y(5) - (-m * sqrt(g^2 + R^2 * y(6)^2 - 2 * g * R * y(6) * sin(s)) * r * sin(s + y(1) - asin(R * y(6) * cos(s) / sqrt((g - R * y(6) * sin(s))^2 + (R * y(6) * cos(s))^2))) - k * R * A * (y(3) - y(4)));
  dydt(6) = 2 * M * R * y(6) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * y(5)));
end

function M = MassFcn()
  M = eye(6);
  M(5,5) = 0;
  M(6,6) = 0;
end

you will get the best performance I could achieve until now.

But you will see that the functions blow up at this time, and continuing integration seems no longer possible.