From: Can on
Hi,
I need help on this fallowing code. I need a change so that there is damping on the power angle. The code is written for oscillatory stable. I would appreciate if you help me out.

delta0=31*pi/180;
X0=[delta0;0];
% Simulation time from 0 to 4 cycles.
TCLEAR=0.2;
TSPAN =0:0.01:TCLEAR;
%call ODE solver to integrate from 0 to TCLEAR
[t, x]=ODE45('go1',TSPAN,X0);
last=length(t);
X0=[x(1); x(2)];
xold=x(1:last-1,:);
told=t(1:last-1);
% Simulation time from 4 cycles to 5 seconds
TSPAN =TCLEAR:0.01:20;
%call ODE solver to integrate from TCLEAR to 5 seconds.
[t, x]=ODE45('go2',TSPAN,X0);
% Now join xhold with x.
delta=[xold(:,1)' x(:,1)']';
delta=delta*180/pi;
omega=[xold(:,2)' x(:,2)']';
ttot=[told' t']';

plot(ttot,delta)
grid
ylabel('Angle, delta (degrees)')
xlabel('Time (seconds)')

function xp1 = go1(t,x)
%
% In: t (time) is a scalar.
% x is a 2-element vector whose elements are:
% x(1)= delta(t)
% x(2)= omega(t)
%
% Out: xp1 is a 2-element vector whose elements are
% xp1(1)=d delta(t)/dt
% xp1(1)=d omega(t)/dt
%
xp1(1)=x(2);
xp1(2)=1-0.976*sin(x(1));
% ODE45 expects the function to return a column vector.
xp1=x;
function xp2 = go2(t,x)
%
% In: t (time) is a scalar.
% x is a 2-element vector whose elements are:
% x(1)= delta(t)
% x(2)= omega(t)
%
% Out: xp2 is a 2-element vector whose elements are
% xp2(1)=d delta(t)/dt
% xp2(1)=d omega(t)/dt
%
xp2(1)=x(2);
xp2(2)=1-1.56*sin(x(1));
% ODE45 expects the function to return a column vector.
xp2=xp2';
 | 
Pages: 1
Prev: working with structures
Next: mupad procedures