将下程序保存为TwoPLL2.m
function xdot=TwoPLL2(t,x,dummy,tau1,tau2,K,Deltaomegao)
xdot=zeros(2,1);
xdot(2)=(-(1/tau1+K*tau2*(1-x(1)^2/2+x(1)^4/24)/tau1)*x(2)-K*(x(1)-x(1)^3/6+x(1)^5/120)/tau1+Deltaomegao/tau1);
xdot(1)=x(2);
程序ode4.m 四阶龙格库塔法
function Y = ode4(odefun,tspan,y0,varargin)
if ~isnumeric(tspan)
error('TSPAN should be a vector of integration steps.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
error(msg);
end
y0 = y0(:); % Make a column vector.
if ~isequal(size(y0),size(f0))
error('Inconsistent sizes of Y0 and f(t0,y0).');
end
neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
F = zeros(neq,4);
Y(:,1) = y0;
for i = 2:N
ti = tspan(i-1);
hi = h(i-1);
yi = Y(:,i-1);
F(:,1) = feval(odefun,ti,yi,varargin{:});
F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:});
F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});
Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));
end
Y = Y.';
主程序:
tau1=5;tau2=1;K=12;Deltaomegao=[1,2,3,4,5,6,7,8,9,10];
tspan=linspace(0,15,1500);
options=odeset('Abstol',[1e-6,1e-6]);
for i=1:10
[t,x]=ode4('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));
figure(1);plot(t,x(:,1)); hold on; title('Phase difference as a function of Time');ylabel('θe(t)/rad');xlabel('t/s');
figure(2);plot(t,x(:,2)); hold on; title('Frequency difference as a function of Time');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('t/s');
figure(3);plot(x(:,1),x(:,2)); hold on; title('Phase portrait');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('/rad');
end
本人找到一个四阶龙格库塔法的程序,但是在主程序中调用时出现如下错误,:@L 请教一下怎么修改?:@L
??? Error using ==> ode4
Too many output arguments.
Error in ==> TwoPLLshzhjiyituduxi at 5
[t,x]=ode4('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i)); |