Could anyone please explain why MATLAB doesn't return answer to this code?

1 view (last 30 days)
clc, clear
% Inputs:
g0 = 1.62;
R = 1;
M = 2;
th = 1;
D = 2;
u0 = [-418.9536 -34.9682 123.2243 119.9791 50.7267 6.1709 3.4923 2.4725 4.0626 3.8980 2.8555];
w0 = 817.5053;
x0 = [R 0 0 0]' ;
[~ ,r] = size(u0);
tt = linspace(0,1,r);
uu = fit(tt' , u0' , 'linearinterp')
[T , X] = ode45(@(taw,x) state_2(taw, x, uu, w0,g0,R,M,th), [0 1], x0);
%-------------------------------------------------------------------------------------------------------------------
% and this is the function "state_2" :
function xdot = state_2(taw,x,uu,w,g0,R,M,th)
xdot = w*[ x(3)
x(4)/x(1)
(x(4)^2)/x(1) - (g0 * R^2)/(x(1)^2) + (th/M)*sin(uu(taw))
-(x(3) * x(4))/x(1) + (th/M)*cos(uu(taw))];
end

Accepted Answer

Walter Roberson
Walter Roberson on 27 Jul 2022
uu = fit(tt' , u0' , 'linearinterp')
linear interpolants do not have continuous derivatives, so using linear interpolant violates the mathematics behind ode45 .
The code does finish running, in my tests, but when you look at the plots you will see a lot of apparent noise for X(:,3) and X(:,4) below t = 0.2 -- where t = 0.2 is the location of the most significant breakpoint in uu.
If you change to 'smoothingspline' then the calculation will eventually finish, but it is slow.
If you change to 'spline' then the calculation is quite fast, but possibly not accurate.
  5 Comments
Walter Roberson
Walter Roberson on 27 Jul 2022
If you have discrete data from evaluating ode45 at various time steps, then you could consider using cubic spline interpolation or smoothing spline, if that would be accurate enough for your purposes. cubic spline can cause a bit of "ringing" near change points -- for example
B-C
A D
then between B and C, cubic spline would predict something that rises above B and C; cubic spline will not handle sharp corners without some projection beyond the available data. (Polynomial fits are often worse than cubic spline for this purpose.)

Sign in to comment.

More Answers (1)

Bruno Luong
Bruno Luong on 27 Jul 2022
Edited: Bruno Luong on 27 Jul 2022
Walter is correct, you should break the interval and do integration on each sequentially
clc, clear
% Inputs:
g0 = 1.62;
R = 1;
M = 2;
th = 1;
D = 2;
u0 = [-418.9536 -34.9682 123.2243 119.9791 50.7267 6.1709 3.4923 2.4725 4.0626 3.8980 2.8555];
w0 = 817.5053;
x0 = [R 0 0 0]' ;
r = size(u0,2);
tt = linspace(0,1,r);
T = [];
X = [];
for k = 1:r-1
tk = tt(k:k+1);
uk = u0(k:k+1);
[Tk , Xk] = ode45(@(taw,x) state_3(taw, x, [tk; uk], ...
w0,g0,R,M,th), tk, x0);
x0 = Xk(end,:);
T = [T; Tk];
X = [X; Xk];
end
for j = 1:size(X,2)
subplot(2,2,j);
plot(T,X(:,j));
title(sprintf('X(:,%d)', j))
end
%-------------------------------------------------------------------------------------------------------------------
% and this is the function "state_3" :
function xdot = state_3(taw,x,tu,w,g0,R,M,th)
% linear interpolation between t(1) and t(2)
t = tu(1,:);
u = tu(2,:);
p = (taw-t(1))./(t(2)-t(1));
uu = (1-p).*u(1) + p.*u(2);
xdot = w*[ x(3)
x(4)/x(1)
(x(4)^2)/x(1) - (g0 * R^2)/(x(1)^2) + (th/M)*sin(uu)
-(x(3) * x(4))/x(1) + (th/M)*cos(uu)];
end

Categories

Find more on Mathematics and Optimization in Help Center and File Exchange

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!