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

1 次查看（过去 30 天）
roham farhadi
roham farhadi 2022-7-27
评论： Walter Roberson 2022-7-27
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

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

采纳的回答

Walter Roberson
Walter Roberson 2022-7-27
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 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2022-7-27
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.)

请先登录，再进行评论。

更多回答（1 个）

Bruno Luong
Bruno Luong 2022-7-27
编辑：Bruno Luong 2022-7-27
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Splines 的更多信息

标签

产品


版本

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by