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

Question

Student 2023-10-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2028665-i-don-t-know-why-this-code-doesn-t-stops

编辑： Torsten 2023-10-5

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2023-10-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2028665-i-don-t-know-why-this-code-doesn-t-stops#answer_1324560

编辑：Torsten 2023-10-3

在 MATLAB Online 中打开

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.

13 个评论
显示 11更早的评论隐藏 11更早的评论

Torsten 2023-10-3

编辑：Torsten 2023-10-3

在 MATLAB Online 中打开

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

Torsten 2023-10-3

编辑：Torsten 2023-10-4

在 MATLAB Online 中打开

@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.1052

ans = -0.1224

ans = 3.9154

ans = 0.1052

ans = -0.1224

ans = 4.8001

ans = 0.0378

ans = -0.0303

ans = 4.7753

ans = 0.0378

ans = -0.0303

ans = 4.8093

ans = 0.0175

ans = -0.0030

ans = 4.8020

ans = 0.0175

ans = -0.0030

ans = 4.8122

ans = 0.0113

ans = 0.0052

ans = -9.3687e-04

ans = 4.8082

ans = 0.0052

ans = -9.3687e-04

ans = 4.8143

ans = 0.0034

ans = 0.0015

ans = -3.2303e-04

ans = 4.8131

ans = 0.0015

ans = -3.2303e-04

ans = 4.8149

ans = 9.6582e-04

ans = 4.1349e-04

ans = -1.3890e-04

ans = 4.8146

ans = 4.1349e-04

ans = 4.1334e-04

ans = -1.3890e-04

ans = 4.8151

ans = 2.4778e-04

ans = 8.2062e-05

ans = -8.3656e-05

ans = 4.8150

ans = 8.2062e-05

ans = 8.1916e-05

ans = -8.3656e-05

ans = 4.8152

ans = 3.2347e-05

ans = -1.7369e-05

ans = 4.8152

ans = 3.2347e-05

ans = 3.2493e-05

ans = -1.7369e-05

ans = 4.8152

ans = 1.7432e-05

ans = 2.5178e-06

ans = -1.2397e-05

ans = 4.8152

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.6686e-07

ans = 4.8152

ans = 1.1755e-06

ans = 1.0293e-06

ans = -1.6686e-07

ans = 4.8152

ans = 7.7277e-07

ans = 3.7007e-07

ans = -3.2625e-08

ans = 4.8152

ans = 3.7007e-07

ans = 5.1620e-07

ans = -3.2625e-08

ans = 4.8152

ans = 2.4926e-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 = 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.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.4599e-07

ans = -8.3741e-10

ans = 1.4529e-07

ans = 8.7497e-07

ans = 3.0659e-06

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 2023-10-5

编辑：Torsten 2023-10-5

在 MATLAB Online 中打开

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