Error using ODE45

Hello everyone,

I used the ODE function to solve a set of differential equations in which a torque was implied (see in the code line1--->line129) and the simulation was made for different values of the torque. After that, i decide to add another stage at my simulation and turn off my moment/torque at a certain time (700s). To do this, i tried to define the torque as an input function but in the command window, an error appears, the function u seems not to be defined.

Undefined function or variable 'u'.

Error in Ryelandt_288861_Assignment2>@(t,x)mydyn(t,x,myinput(t,u(1)))

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Ryelandt_288861_Assignment2 (line 134)

[~,Xout11] = ode45(@(t,x)mydyn(t,x,myinput(t,u(1))),xx,IC1,opts);

tic; 
%% ME - 221 Dynamical systems, Spring 2020
% Solution to MATLAB Assignment #2
% Student Name: Ryelandt Guillaume
% SCIPER 288861
% Submission date: 10 April 2020
%% Initiate the script 
clc; clear; close all

%% Question 3 a
%initial conditions 
xx = linspace(0,1000,1e3).';
M=[5,15,25,55];         %[Nm]
IC1=[298; 298; 0; 0];
opts= odeset('RelTol',1e-8,'AbsTol',1e-10);
[~,Xout1] = ode45(@(t,x)mydyn(t,x,M(1)),xx,IC1,opts);
[~,Xout2] = ode45(@(t,x)mydyn(t,x,M(2)),xx,IC1,opts);
[~,Xout3] = ode45(@(t,x)mydyn(t,x,M(3)),xx,IC1,opts);
[~,Xout4] = ode45(@(t,x)mydyn(t,x,M(4)),xx,IC1,opts);

%data visualisation
%for m=5
figure(1);
    plot(xx,Xout1(:,1),'linewidth',2,'displayName','M= 5 Nm');
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$T$(K)','interpreter','Latex')
    set(gca,'FontSize',11)
hold on; 

    plot(xx,Xout2(:,1),'linewidth',2,'displayName','M= 15 Nm');

    plot(xx,Xout3(:,1),'linewidth',2,'displayName','M= 25 Nm');

    plot(xx,Xout4(:,1),'linewidth',2,'displayName','M= 55 Nm');
    
    lgd=legend;
% [Xs,YS]=ginput(4);
hold off; 
%% Question 3 b
 figure(2);
Xout1norm=normalize(Xout1(:,1),'range');
Xout2norm=normalize(Xout2(:,1),'range');
Xout3norm=normalize(Xout3(:,1),'range');
Xout4norm=normalize(Xout4(:,1),'range');


    plot(xx,Xout1norm,'linewidth',2,'displayName','M= 5 Nm');
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$T$(K)','interpreter','Latex')
    set(gca,'FontSize',11)
hold on; 

    plot(xx,Xout2norm,'linewidth',2,'displayName','M= 15 Nm');
    
    plot(xx,Xout3norm,'linewidth',2,'displayName','M= 25 Nm');
    
    plot(xx,Xout4norm,'linewidth',2,'displayName','M= 55 Nm');
    
    lgd=legend; 
hold off; 
% yyaxis right
%     plot(tt,Xout1(:,2),'linewidth',2);
%     ylabel('$\dot{x}$(m/s)','interpreter','Latex')
%     set(gca,'FontSize',11)
% set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
%% Question 3 c 
figure(3);


    plot(xx,Xout1(:,4),'linewidth',2,'displayName','M= 5 Nm');
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$v$(rad/s)','interpreter','Latex')
    set(gca,'FontSize',11)
hold on; 

    plot(xx,Xout2(:,4),'linewidth',2,'displayName','M= 15 Nm');
    
    plot(xx,Xout3(:,4),'linewidth',2,'displayName','M= 25 Nm');

    plot(xx,Xout4(:,4),'linewidth',2,'displayName','M= 55 Nm');
    
    lgd=legend;
    
hold off; 
figure(4);
Xout1norm=normalize(Xout1(:,4),'range');
Xout2norm=normalize(Xout2(:,4),'range');
Xout3norm=normalize(Xout3(:,4),'range');
Xout4norm=normalize(Xout4(:,4),'range');

figure (4); 
    plot(xx,Xout1norm,'linewidth',2,'displayName','M= 5 Nm');
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$v$(rad/s)','interpreter','Latex')
    set(gca,'FontSize',11)
hold on;

    plot(xx,Xout2norm,'linewidth',2,'displayName','M= 15 Nm');

    plot(xx,Xout3norm,'linewidth',2,'displayName','M= 25 Nm');

    plot(xx,Xout4norm,'linewidth',2,'displayName','M= 55 Nm');
    
    lgd=legend; 
hold off; 
%% Question 3 d
RelError1=100*abs((Xout1(:,2)-Xout1(:,1))./(Xout1(:,1))); 
RelError2=100*abs((Xout2(:,2)-Xout2(:,1))./(Xout2(:,1)));
RelError3=100*abs((Xout3(:,2)-Xout3(:,1))./(Xout3(:,1)));
RelError4=100*abs((Xout4(:,2)-Xout4(:,1))./(Xout4(:,1)));

figure (5);
plot(xx,RelError1,'linewidth',2,'displayName','M= 5 Nm');
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$Relative Error(\%)$','interpreter','Latex')
    set(gca,'FontSize',11)
    
hold on;

    plot(xx,RelError2,'linewidth',2,'displayName','M= 15 Nm');
 
    plot(xx,RelError3,'linewidth',2,'displayName','M= 25 Nm');

    plot(xx,RelError4,'linewidth',2,'displayName','M= 55 Nm');
   
    
    lgd=legend;
toc;
%%
xx = linspace(0,1000,1e3).';
                                    
IC1=[298; 298; 0; 0];
opts= odeset('RelTol',1e-8,'AbsTol',1e-10);
[~,Xout11] = ode45(@(t,x)mydyn(t,x,myinput(t,u(1))),xx,IC1,opts);
[~,Xout21] = ode45(@(t,x)mydyn(t,x,myinput(t,u(2))),xx,IC1,opts);
[~,Xout31] = ode45(@(t,x)mydyn(t,x,myinput(t,u(3))),xx,IC1,opts);
[~,Xout41] = ode45(@(t,x)mydyn(t,x,myinput(t,u(4))),xx,IC1,opts);
%% Function definition
function dx = mydyn(t,x,M)
% Define the parameters 
Ct= 200;                %[J/K]
ht=140;                 %[W/m^2K]
At=0.02;                %[m^2]
he=200;                 %[W/m^2K]
Ae=pi;                  %[m^2]
Te= 298;                %[K]
V=6*pi;                 %[m^3]
rho=1;                  %[kg/m^3]
cp=4200;                %[J/kgK]
J=300;                  %[kgm^2]
c=2;                    %[Nms]
c3=0.1;                 %[Nms^3]
cq=0.09;                %[Ws^4]


T=x(1);
Tt=x(2);
theta=x(3);
thetadot=x(4);

u=M; 
%state model
dT=(V*rho*cp)^-1*(he.*Ae.*(Te-T)+ht*At*(Tt-T)+cq*thetadot^4);
dTt=(ht*At)*(Ct)^-1.*(T-Tt);
dtheta=thetadot;
dthetadot=(M-c*thetadot-c3*thetadot^3)/J;

dx=[dT,dTt,dtheta,dthetadot]';
end
function u=myinput(t,M)

if t<700 
    u=[5,15,25,55];         %[Nm]
else t>700
    u=[0,0,0,0];         %[Nm]
end 
end 

>>

Answer 1

[darova]

The name of function is myinput

function u=myinput(t,M)

But you are trying to pass as argument some u

[~,Xout11] = ode45(@(t,x)mydyn(t,x,myinput(t,u(1))),xx,IC1,opts);

You don't call function inside mydyn

u=M;            % what is this?

You don't use M argument inside

function u=myinput(t,M)
if t<700 
    u=[5,15,25,55];         %[Nm]
else t>700
    u=[0,0,0,0];         %[Nm]
end 
end 

here is the solution

[Guillaume Ryelandt]

thank you for your help, i really appreciate this!! However, i do not understand what you are doing with the a(M) and why you just replaced u(1) by 1. I actually tried your solution, but it still doesn't work. In the first part of the code, i just set the value of my input as M (the torque) like this :

M=[5,15,25,55];         %[Nm]
IC1=[298; 298; 0; 0];
opts= odeset('RelTol',1e-8,'AbsTol',1e-10);
[~,Xout1] = ode45(@(t,x)mydyn(t,x,M(1)),xx,IC1,opts);
[~,Xout2] = ode45(@(t,x)mydyn(t,x,M(2)),xx,IC1,opts);
[~,Xout3] = ode45(@(t,x)mydyn(t,x,M(3)),xx,IC1,opts);
[~,Xout4] = ode45(@(t,x)mydyn(t,x,M(4)),xx,IC1,opts);

and then i used the ode45 function. Now the problem is that i need to keep these value for the torque until i reach 700s, and at this moment, i just turn my torque off and continue to solve the differential equation with ode45. That is actually why i tried this input function but it doesn't seem to work...

[darova]

That's smart!

[~,Xout1] = ode45(@(t,x)mydyn(t,x,M(1)),xx,IC1,opts);

So just add this line inside your mydyn function

[Guillaume Ryelandt]

Ok thanks!

Answer 2

[Steven Lord]

A cleaner (IMO) approach would be to solve your system of ODEs with torque applied from time t = 0 to time t = 700. Use the result of that first ode45 call at time t = 700 as the initial condition for a second call to ode45 that solves a system of ODEs without torque from time t = 700 to time t = 1000.

[Guillaume Ryelandt]

Thank you for your help, i will try this!