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

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

Answer 1

[Walter Roberson]

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.

[roham farhadi]

thank you for your answer,

so how can I use an ode45 answer into another ode45 without using curve fitting?

[Bruno Luong]

If your question is about : solving the first ode, obtain a discrete solution, use interpolation (not fit) as a continuous function that will be used in the second ode.

Then I think it is better to write both odes as a single ode with both state variables concatenated as a longer but single state variable.

[roham farhadi]

thank you for your answer,

but I can't do this, because the first ode45 has Initial conditions and the second ode45 has final conditions.

[Bruno Luong]

You might take a look at bvp4c instead of ode45

[Walter Roberson]

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.)

Answer 2

[Bruno Luong]

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