I don't know why this code doesn't stops :(
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
Answers (1)
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
Student
on 3 Oct 2023
Torsten
on 3 Oct 2023
Do you know how the functions should behave for greater values of t ? Can you sketch y and y' ?
Student
on 3 Oct 2023
Student
on 3 Oct 2023
these are the two differential equations and the value of constant is mentioned above
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
on 3 Oct 2023
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
Steven Lord
on 3 Oct 2023
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.
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.
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
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.
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
