Error using ODE45

1 次查看（过去 30 天）
Guillaume Ryelandt
评论： Guillaume Ryelandt 2020-4-7
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
>>

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

回答（2 个）

darova
darova 2020-4-7
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 个评论
显示 1更早的评论
darova
darova 2020-4-7
That's smart!
[~,Xout1] = ode45(@(t,x)mydyn(t,x,M(1)),xx,IC1,opts);
So just add this line inside your mydyn function

请先登录，再进行评论。

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

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by