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

2013 7 29
2 个回答
回答已采纳
4 次查看（30 天）

1 个投票

Sir, I am calling ode45 for integrating a diff equation say x'=(1/(H))*(2x)
Main Program:
x0=0.1; final_time=1; [t,x]=ode45(@myfunc,[0,final_time],x0);% ODE Call plot(t(:,1),x(:,1))% plotting the state variable main ends
Function Program:
function dv=myfunc(t,x) H=1.0; Te=2*x(1);% intermediate calculations/variables for my differential equations %The Differential Equation dv=[(1/(H))*(Te)];
But i want to plot the intermediate variables say Te in this case.But ode returns t,x alone.How to plot the intermediate variables that was used in the integration process.

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

 采纳的回答

Jan
Jan 2013-7-29

0 个投票

You can prepare your function such, that it replies the wanted value according to a special trigger method and accepts vectors as input also:
function dv = myfunc(t, x, flag)
H = 1.0;
Te = 2 * x(:, 1)
dv = Te ./ H; % Btw, No addition square brackets
if nargin == 3
dv = Te;
end
Now run the integration at first:
[t,x] = ode45(@myfunc,[0,final_time],x0);
and obtain the internal variable afterwards:
Te = myFunc(t, x, 'flag');

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Thayumanavan
Thayumanavan 2013-7-30
Sir, Thanks for the immediate response.I tried your code and got the output plot of Te .
Inference:
In workspace i am not getting the array of Te values for all time instants.Getting a value of 4.5 (1x 1 Element)only 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.(1057x1) entries.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.I am getting the plot as, with respect to time (0 to 10) a line at 4.5 .so i guess it is plotting t with respect to the last value of Te alone say Te=4.5 (at t=10).(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 .
Jan
Jan 2013-8-1
I do not understand your comment. Te does not depend on the time, but is simply twice as large as x in your case. Therefore I cannot understand how you get results, which "Te seems to be constant with respect to time".
Please post a copy of the code you use.
Thayumanavan
Thayumanavan 2013-8-3
clear all clc; % Time Tspan =[0 10]; X0=0.7; options=odeset('RelTol',1e-3,'AbsTol',1e-6); %solving the differential Equation [t,x,s]=ode45(@myfunc,Tspan,X0,options);
%Plotting the figure
subplot(2,1,1) plot(t(:,1),x(:,1)) xlabel('t'); ylabel('x'); title('Time Response ') grid on
subplot(2,1,2) IntVar= dvdt(t, x, 'flag'); plot(t(:,1),x(1,:)) xlabel('t'); ylabel('Te') grid on
function dv=myfunc(t,x,flag) H = 1.0; Te = 2 * x(:, 1)
%The Diffenrential Equation
dv=Te./H;
if nargin == 3
dv = Te;
end
This is the code ,what i am expecting is t=0 to 10 ,x is changing means Te should also change.But Te seems to be constant.You try the code sir It is plotting the last value of Te(0.999) wrt time.Hope you understand.
Thayumanavan
Thayumanavan 2013-8-3
Sir ,I got answer for the problem that i asked help.But for simplicity i gave that differential equation.But my code is little bit different with many variables.I am attaching the main program and my function.Please help me to plot Te and Vpcc intermediate variable.
Main Code: 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('Dynamic Wind 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
%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
You copy paste the code execute and suggest me to plot Te,Vpcc.Thanks
Thayumanavan
Thayumanavan 2013-8-3
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('Dynamic Wind 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
Little change in code .i forgot to give 1 parameter .You pls simulate and tell me what change i need to do in the code to get plot of Te,Vpcc

请先登录，再进行评论。

更多回答（1 个）

Shashank Prasanna
Shashank Prasanna 2013-7-29

0 个投票

You can define an output function (outputfcn) that will be called after each iteration.
http://www.mathworks.com/help/matlab/ref/odeset.html#f92-1016858
you can specify that in the options structure. In your custom output function you can use a plot command with a hold on in order to plot all you intermediate values.

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by