MATLAB Answers

Error using ODE45

4 views (last 30 days)
Guillaume Ryelandt
Guillaume Ryelandt on 7 Apr 2020
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
>>

  0 Comments

Sign in to comment.

Answers (2)

darova
darova on 7 Apr 2020
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

  3 Comments

Guillaume Ryelandt
Guillaume Ryelandt on 7 Apr 2020
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
darova on 7 Apr 2020
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
Guillaume Ryelandt on 7 Apr 2020
Ok thanks!

Sign in to comment.


Steven Lord
Steven Lord on 7 Apr 2020
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.

  1 Comment

Guillaume Ryelandt
Guillaume Ryelandt on 7 Apr 2020
Thank you for your help, i will try this!

Sign in to comment.

Sign in to answer this question.


Translated by