Filter löschen
Filter löschen

How to plot intermediate variables of a function used by ode45 solver:

1 Ansicht (letzte 30 Tage)
Thayumanavan
Thayumanavan am 1 Aug. 2013
Sir, I want to plot the intermediate variables say Te from a function that was used by ode45 solver .But i am not getting the response exactly. The time response of x is not constant(from the plot of (t,x(:,1))). But the The time response of Te seems to be constant with respect to time. The problem is ,in workspace i am not getting the array of Te values for all time instants.I am getting a single last value of Te say 4.5 (1x 1 Element) Workspace is having the last value of Te alone ie. at t=10 if tspan=[0 10]. But the vector of state variable x is coming in workspace.(say 1057x1) entries.(If x is time varying means Te should also vary with respect to time) . I tried giving Te as global variable also.But in vein.Please help me .
Main Pragram:
x0=0.1;
final_time=1;
[t,x]=ode45(@myfunc,[0,final_time],x0);% ODE Call
plot(t(:,1),x(:,1))% plotting the state varaible
Intervar=myfunc(t,x,'flag');
plot(t(:,1),Intervar)
main ends
Function :
function dv=myfunc(t,x)
H=1.0;
Te=2*x(1);% intermediate calculations/variables for my differential equations
%The Diffenrential Equation
dv=[(1/(H))*(Te)];
if nargin == 3
dv=Te; //Returns the intermediate variable Te
end

Antworten (2)

Azzi Abdelmalek
Azzi Abdelmalek am 1 Aug. 2013
Bearbeitet: Azzi Abdelmalek am 1 Aug. 2013
In your code Te represent 2*x, then why you don't just type
Te=2*x(:,2)

Jan
Jan am 1 Aug. 2013
Bearbeitet: Jan am 1 Aug. 2013
Please copy the code I have posted exactly:
function dv = myfunc(t, x, flag)
H = 1.0;
Te = 2 * x(:, 1); % <-- Not "Te=2*x(1);"
dv = Te ./ H; % Btw, No additional square brackets
if nargin == 3
dv = Te;
end
  1 Kommentar
Thayumanavan
Thayumanavan am 18 Aug. 2013
Sir, my code is this:Main Program
clear all clc; % Time Tspan =[0 10]; X0=[0.7; 1.0151; 157 ];
options=odeset('RelTol',1e-3,'AbsTol',1e-6);
%solving the differential Equation [t,x]=ode45(@dvdt,Tspan,X0,options);
%Plotting the figures figure subplot(5,1,1) plot(t(:,1),x(:,1)) xlabel('t'); ylabel('Eprime'); title('Time Response ') grid on
subplot(5,1,2) plot(t(:,1),x(:,2)) xlabel('t'); ylabel('Delta') grid on
wradsec=x(:,3); wmech=(60*wradsec)/(2*pi); wpu=wmech/1500;
subplot(5,1,3) plot(t(:,1),wradsec) xlabel('t'); ylabel('Omega') grid on
% Plotting of Intermediate Variables as per Your Idea subplot(5,1,4) op= dvdt(t, x, 'flag'); plot(t(:,1),op(1,:)) xlabel('t'); ylabel('Vterminal') grid on
subplot(5,1,5) plot(t(:,1),op(2,:)) xlabel('t'); ylabel('ElectromagneticTorque') grid on
Function Code:
function dv=dvdt(t,x,flag)
Xs=0.0754;Xr=0.0326;H=0.20; Rl=0.0083;Rcon=0.0089;Xl=0.033;Xcon=0.01789; Vdinf=1.0;Vqinf=0.0;wspu=157.0; T=0.7855;Xdash=0.1076;D=6.4165; Rt=Rl+Rcon;Xt=Xl+Xcon; xtrans=Xt+(D/Xr); detm=(Rt*Rt)+(xtrans*xtrans); J1=-xtrans;J2=-Rt;J3=-Rt;J4=xtrans; iq=(1/detm)*((J3*((x(1)*cos(x(2)))-Vqinf))+(J4*((-x(1)*sin(x(2)))-Vdinf))); id=(1/detm)*((J2*((-x(1)*sin(x(2)))-Vdinf))+(J1*((x(1)*cos(x(2)))-Vqinf))); Vq=(1/Xr)*((D*id)+(Xr*x(1)*cos(x(2)))); Vd=(-1/Xr)*((D*iq)+(Xr*x(1)*sin(x(2)))); Vdpc=(-Rcon*id)+(Xcon*iq)+Vdinf; Vqpc=(-Rcon*iq)-(Xcon*id)+Vqinf; VPCC=sqrt((Vdpc*Vdpc)+(Vqpc*Vqpc));% Intermediate Variable 1 t1=(Vq*cos(x(2)))-(Vd*sin(x(2))); t2=(Vq*sin(x(2)))+(Vd*cos(x(2))); t3=((Vd*cos(x(2)))+(Vq*sin(x(2)))); Te=(x(1)*t3/Xdash)% Intermediate Variable 2 Tm=1.0;
%The Diffenrential Equation
dv=[ (1/T)*((-Xs*x(1)/Xdash)+(((Xs-Xdash)/Xdash)*t1)); % My 1st Differential Equation
(x(3)-wspu)-(((Xs-Xdash)*t2) /(Xdash*T*x(1))); %My 2nd Differential Equation
(wspu/(2*H))*(Tm+Te)]; % My 3rd Differential Equation
if nargin == 3
dv = [VPCC;
Te];
end
It is giving the last value only.In workspace Te,VPCC is a 1x1 element.iT SHOULD BE AN ARRAY similar to x.I tried giving Te,VPCC Global,but in vein.Simulate the code ,Please help.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Particle & Nuclear Physics finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by