I don't know why this code doesn't stops :(
2 views (last 30 days)
Show older comments
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])
[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 Comments
Answers (1)
Torsten
on 3 Oct 2023
Edited: Torsten
on 3 Oct 2023
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 Comments
Torsten
on 3 Oct 2023
Edited: Torsten
on 4 Oct 2023
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))
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));
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
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
on 5 Oct 2023
Edited: Torsten
on 5 Oct 2023
If you download
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!