Could anyone please explain why MATLAB doesn't return answer to this code?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
roham farhadi
am 27 Jul. 2022
Kommentiert: Walter Roberson
am 27 Jul. 2022
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
0 Kommentare
Akzeptierte Antwort
Walter Roberson
am 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 Kommentare
Walter Roberson
am 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.)
Weitere Antworten (1)
Bruno Luong
am 27 Jul. 2022
Bearbeitet: Bruno Luong
am 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
0 Kommentare
Siehe auch
Kategorien
Mehr zu Mathematics and Optimization finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!